using CSV using Rasters using NighttimeLights using Plots using DataFrames using Dates capex = CSV.read("./DATA/29_thermal_2015_2018_manually_checked_construction_lat_long.csv", DataFrame) radiance_path = "/home/nighttimelights/VIIRS_ANNUAL/v21_median/" dates = range(start = Date(2012), step = Year(1), length = length(readdir(radiance_path))) |> collect timestamps = year.(dates) function agg_radiance_annual(lat, long, R = 10, res = 15) raster1 = Raster(radiance_path * readdir(radiance_path)[1]) bounds, m = annular_ring(raster1, R, lat, long, res) radiances = [Raster(i, lazy = true)[bounds...] for i in radiance_path .* readdir(radiance_path)] series = RasterSeries(radiances, Ti(timestamps)) radiance_datacube = Rasters.combine(series, Ti) return [sum(skipmissing(Array(radiance_datacube[:,:,1,i] .* m))) for i in 1:length(readdir(radiance_path))] end ## 1 km aggs = [] for i in capex.lat_long lat, long = tuple(parse.(Float64, split(i, ","))...) push!(aggs, agg_radiance_annual(lat, long, 1)) end agg_matrix = cat(aggs..., dims = 2)' agg_df = DataFrame(Array(agg_matrix), :auto) rename!(agg_df, Symbol.(dates)) agg_df.lat_long = capex.lat_long capexnl1 = leftjoin(capex, agg_df, on = :lat_long) CSV.write("./DATA/29_thermal_2015_2018_annual_radiance_1km.csv", capexnl1) ## 2 km aggs = [] for i in capex.lat_long lat, long = tuple(parse.(Float64, split(i, ","))...) push!(aggs, agg_radiance_annual(lat, long, 2)) end agg_matrix = cat(aggs..., dims = 2)' agg_df = DataFrame(Array(agg_matrix), :auto) rename!(agg_df, Symbol.(dates)) agg_df.lat_long = capex.lat_long capexnl2 = leftjoin(capex, agg_df, on = :lat_long) CSV.write("./DATA/29_thermal_2015_2018_annual_radiance_2km.csv", capexnl2) ## 4 km aggs = [] for i in capex.lat_long lat, long = tuple(parse.(Float64, split(i, ","))...) push!(aggs, agg_radiance_annual(lat, long, 4)) end agg_matrix = cat(aggs..., dims = 2)' agg_df = DataFrame(Array(agg_matrix), :auto) rename!(agg_df, Symbol.(dates)) agg_df.lat_long = capex.lat_long capexnl = leftjoin(capex, agg_df, on = :lat_long) CSV.write("./DATA/29_thermal_2015_2018_annual_radiance_4km.csv", capexnl) ## 8 km aggs = [] for i in capex.lat_long lat, long = tuple(parse.(Float64, split(i, ","))...) push!(aggs, agg_radiance_annual(lat, long, 8)) end agg_matrix = cat(aggs..., dims = 2)' agg_df = DataFrame(Array(agg_matrix), :auto) rename!(agg_df, Symbol.(dates)) agg_df.lat_long = capex.lat_long capexnl = leftjoin(capex, agg_df, on = :lat_long) CSV.write("./DATA/29_thermal_2015_2018_annual_radiance_8km.csv", capexnl) ## 16 km aggs = [] for i in capex.lat_long lat, long = tuple(parse.(Float64, split(i, ","))...) push!(aggs, agg_radiance_annual(lat, long, 16)) end agg_matrix = cat(aggs..., dims = 2)' agg_df = DataFrame(Array(agg_matrix), :auto) rename!(agg_df, Symbol.(dates)) agg_df.lat_long = capex.lat_long capexnl = leftjoin(capex, agg_df, on = :lat_long) CSV.write("./DATA/29_thermal_2015_2018_annual_radiance_16km.csv", capexnl) ## 20 km aggs = [] for i in capex.lat_long lat, long = tuple(parse.(Float64, split(i, ","))...) push!(aggs, agg_radiance_annual(lat, long, 20)) end agg_matrix = cat(aggs..., dims = 2)' agg_df = DataFrame(Array(agg_matrix), :auto) rename!(agg_df, Symbol.(dates)) agg_df.lat_long = capex.lat_long capexnl = leftjoin(capex, agg_df, on = :lat_long) CSV.write("./DATA/29_thermal_2015_2018_annual_radiance_20km.csv", capexnl) ## 25 km aggs = [] for i in capex.lat_long lat, long = tuple(parse.(Float64, split(i, ","))...) push!(aggs, agg_radiance_annual(lat, long, 25)) end agg_matrix = cat(aggs..., dims = 2) agg_df = DataFrame(Array(agg_matrix), :auto) rename!(agg_df, Symbol.(dates)) agg_df.lat_long = capex.lat_long capexnl = leftjoin(capex, agg_df, on = :lat_long) CSV.write("./DATA/29_thermal_2015_2018_annual_radiance_25km.csv", capexnl)