Merge pull request #3 from JuliaPlanet/dev

Dev

src/NighttimeLights.jl

@@ -5,7 +5,7 @@ using Dates, ArchGDAL, GeoArrays, DataFrames, Shapefile, GeoInterface, ProgressM
 export load_img, load_datacube, save_img, save_datacube, 
 lat_to_row, row_to_lat, lat_to_row, long_to_column, column_to_long, translate_geometry, Coordinate, CoordinateSystem, image_to_coordinate,
 polygon_mask, load_shapefile, mask_area, LAT, EW, NS, MODEL_EW, MODEL_NS, long_apply, cross_apply, apply_mask, view_img, aggregate, aggregate_timeseries, aggregate_dataframe, bias_correction, bias_correction_datacube, outlier_mask, outlier_ts, linear_interpolation,
-conventional_cleaning, PatnaikSTT2021, threshold_datacube, background_noise_mask, sparse_cube, mark_nan, make_datacube, load_example, radiance_datacube, clouds_datacube, mumbai_map, TILE1_COORDINATE_SYSTEM, TILE2_COORDINATE_SYSTEM, TILE3_COORDINATE_SYSTEM, TILE4_COORDINATE_SYSTEM, TILE5_COORDINATE_SYSTEM, TILE6_COORDINATE_SYSTEM, FULL_COORDINATE_SYSTEM_COORDINATE_SYSTEM, INDIA_COORDINATE_SYSTEM, MUMBAI_COORDINATE_SYSTEM, coordinate_to_image
+conventional_cleaning, PatnaikSTT2021, threshold_datacube, background_noise_mask, sparse_cube, mark_nan, make_datacube, load_example, radiance_datacube, clouds_datacube, mumbai_map, TILE1_COORDINATE_SYSTEM, TILE2_COORDINATE_SYSTEM, TILE3_COORDINATE_SYSTEM, TILE4_COORDINATE_SYSTEM, TILE5_COORDINATE_SYSTEM, TILE6_COORDINATE_SYSTEM, FULL_COORDINATE_SYSTEM_COORDINATE_SYSTEM, INDIA_COORDINATE_SYSTEM, MUMBAI_COORDINATE_SYSTEM, coordinate_to_image, aggregate_per_area_dataframe
 
 
 include("data_io.jl")

src/aggregate.jl

@@ -21,29 +21,60 @@ aggregate_timeseries(rand_datacube, rand_mask)
 ```
 """
 function aggregate_timeseries(datacube, mask)
- datacube = convert(Array{Float32,3}, datacube)
- datacube = sparse_cube(copy(datacube))
- datacube = apply_mask(datacube, mask)
- lights = cross_apply(sum, datacube)
- return lights
+ if typeof(datacube) == Array{Float16, 3}
+ @warn "Float 16 may not be sufficient to aggregate."
+ end
+ tmp = apply_mask(copy(datacube), mask)
+ lights = []
+ for i in 1:length(datacube)
+ push!(lights, sum(tmp[i]))
+ end
+ return Float32.(lights)
 end
 
 """
 Computes the time series of aggregate value of datacube for each row of a shapefile
 # Example:
 ```
-julia> rand_datacube = rand(10, 10, 10)
-julia> shapefile_df = load_shapefile("assets/mumbai_map/mumbai_districts.shp")
-julia> aggregate_dataframe(MUMBAI_COORDINATE_SYSTEM, rand_datacube, shapefile_df, "District")
+rand_datacube = rand(10, 10, 10)
+shapefile_df = load_shapefile("assets/mumbai_map/mumbai_districts.shp")
+aggregate_dataframe(MUMBAI_COORDINATE_SYSTEM, rand_datacube, shapefile_df, "District")
 ```
 """
 function aggregate_dataframe(geometry::CoordinateSystem, datacube, shapefile_df, attribute)
+ if typeof(datacube) == Array{Float16, 3}
+ @warn "Float 16 may not be sufficient to aggregate."
+ end
+ datacube = sparse_cube(datacube)
+ df = DataFrame()
+ @showprogress for i in 1:length(shapefile_df[:, 1])
+ shapefile_row = shapefile_df[i, :]
+ geom_polygon = polygon_mask(geometry, shapefile_row)
+ df[!, shapefile_df[!, attribute][i]] = aggregate_timeseries(datacube, geom_polygon)
+ i = i + 1
+ end
+ return df
+end
+
+"""
+Computes the time series of aggregate value per area of datacube for each row of a shapefile
+# Example:
+```
+rand_datacube = rand(10, 10, 10)
+shapefile_df = load_shapefile("assets/mumbai_map/mumbai_districts.shp")
+aggregate_per_area_dataframe(MUMBAI_COORDINATE_SYSTEM, rand_datacube, shapefile_df, "District")
+```
+"""
+function aggregate_per_area_dataframe(geometry::CoordinateSystem, datacube, shapefile_df, attribute)
+ if typeof(datacube) == Array{Float16, 3}
+ @warn "Float 16 may not be sufficient to aggregate."
+ end
+ datacube = sparse_cube(datacube)
  df = DataFrame()
  @showprogress for i in 1:length(shapefile_df[:, 1])
  shapefile_row = shapefile_df[i, :]
  geom_polygon = polygon_mask(geometry, shapefile_row)
- data = aggregate_timeseries(datacube,geom_polygon)
- df[!, shapefile_df[!, attribute][i]] = data
+ df[!, shapefile_df[!, attribute][i]] = aggregate_timeseries(datacube, geom_polygon) ./mask_area(geometry, geom_polygon)
  i = i + 1
  end
  return df

src/other/masks.jl

@@ -25,11 +25,10 @@ function apply_mask(data, mask = ones((size(data)[1], size(data)[2])))
  return data .* mask
  end
  if typeof(data) == VectorOfArray{Any,3,Array{Any,1}}
- masked_data = []
  for i in 1:length(data)
- push!(masked_data,data[i] .* mask)
+ data[i] = data[i] .* mask
  end
- return VectorOfArray(masked_data)
+ return data
  end
  masked_data = copy(data)
  for i in 1:size(data)[3]