Skip to content

Commit

Permalink
Merge pull request #90 from wri/CIF-293-Set-tile-size-in-meters-for-f…
Browse files Browse the repository at this point in the history 
…ishnet-functions

selectively backed out conversion to meters
  • Loading branch information
@kcartier-wri
kcartier-wri authored Nov 5, 2024
2 parents 845ec8a + 32dcff7 commit 052d978
Show file tree
Hide file tree
Showing 4 changed files with 95 additions and 79 deletions.
52 changes: 32 additions & 20 deletions city_metrix/layers/layer.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@
import shapely.geometry as geometry
import pandas as pd

MAX_TILE_SIDE_SIZE_METERS = 60000 # Approx. 0.5 degrees latitudinal offset. TODO How selected?
MAX_TILE_SIZE_DEGREES = 0.5 # TODO Why was this value selected?

class Layer:
def __init__(self, aggregate=None, masks=[]):
Expand Down Expand Up @@ -55,22 +55,22 @@ def groupby(self, zones, layer=None):
"""
return LayerGroupBy(self.aggregate, zones, layer, self.masks)

def write(self, bbox, output_path, tile_side_meters=None, tile_buffer_meters=None, **kwargs):
def write(self, bbox, output_path, tile_degrees=None, buffer_size=None, **kwargs):
"""
Write the layer to a path. Does not apply masks.
:param bbox: (min x, min y, max x, max y)
:param output_path: local or s3 path to output to
:param tile_side_meters: optional param to tile the results into multiple files specified as tile length on a side in meters
: param tile_buffer_meters: tile buffer distance in meters
:param tile_degrees: optional param to tile the results into multiple files specified as tile length on a side
:param buffer_size: tile buffer distance
:return:
"""

if tile_side_meters is None:
if tile_degrees is None:
clipped_data = self.aggregate.get_data(bbox)
write_layer(output_path, clipped_data)
else:
tile_grid, unbuffered_tile_grid = _get_tile_boundaries(bbox, tile_side_meters, tile_buffer_meters)
tile_grid, unbuffered_tile_grid = _get_tile_boundaries(bbox, tile_degrees, buffer_size)

# write raster data to files
if not os.path.exists(output_path):
Expand All @@ -92,13 +92,13 @@ def write(self, bbox, output_path, tile_side_meters=None, tile_buffer_meters=Non
_write_tile_grid(unbuffered_tile_grid, output_path, 'tile_grid_unbuffered.geojson')


def _get_tile_boundaries(bbox, tile_side_meters, tile_buffer_meters):
has_buffer = True if tile_buffer_meters is not None and tile_buffer_meters != 0 else False
def _get_tile_boundaries(bbox, tile_degrees, buffer_size):
has_buffer = True if buffer_size is not None and buffer_size != 0 else False
if has_buffer:
tiles = create_fishnet_grid(*bbox, tile_side_meters, tile_buffer_meters)
unbuffered_tiles = create_fishnet_grid(*bbox, tile_side_meters)
tiles = create_fishnet_grid(*bbox, tile_degrees, buffer_size)
unbuffered_tiles = create_fishnet_grid(*bbox, tile_degrees)
else:
tiles = create_fishnet_grid(*bbox, tile_side_meters)
tiles = create_fishnet_grid(*bbox, tile_degrees)
unbuffered_tiles = None

tile_grid = []
Expand Down Expand Up @@ -136,7 +136,7 @@ def count(self):
return self._zonal_stats("count")

def _zonal_stats(self, stats_func):
if box(*self.zones.total_bounds).area <= MAX_TILE_SIDE_SIZE_METERS**2:
if box(*self.zones.total_bounds).area <= MAX_TILE_SIZE_DEGREES**2:
stats = self._zonal_stats_tile(self.zones, [stats_func])
else:
stats = self._zonal_stats_fishnet(stats_func)
Expand All @@ -160,7 +160,7 @@ def group_layer_values(df):

def _zonal_stats_fishnet(self, stats_func):
# fishnet GeoDataFrame into smaller tiles
fishnet = create_fishnet_grid(*self.zones.total_bounds, MAX_TILE_SIDE_SIZE_METERS)
fishnet = create_fishnet_grid(*self.zones.total_bounds, MAX_TILE_SIZE_DEGREES)

# spatial join with fishnet grid and then intersect geometries with fishnet tiles
joined = self.zones.sjoin(fishnet)
Expand Down Expand Up @@ -257,17 +257,29 @@ def get_utm_zone_epsg(bbox) -> str:
return f"EPSG:{epsg}"


def create_fishnet_grid(min_lon, min_lat, max_lon, max_lat, tile_side_meters, tile_buffer_meters=0):
def create_fishnet_grid(min_lon, min_lat, max_lon, max_lat, cell_size, buffer_size=0, tile_units_in_degrees=True):
lon_coord, lat_coord = (min_lon, min_lat)
geom_array = []

center_lat = (min_lat + max_lat) / 2
lon_side_offset, lat_side_offset = offset_meters_to_geographic_degrees(center_lat, tile_side_meters)
if tile_buffer_meters == 0:
lon_buffer_offset = 0
lat_buffer_offset = 0
if tile_units_in_degrees:
if cell_size > 0.5:
raise Exception('Value for cell_size must be < 0.5 degrees.')

lon_side_offset = cell_size
lat_side_offset = cell_size
lon_buffer_offset = buffer_size
lat_buffer_offset = buffer_size
else:
lon_buffer_offset, lat_buffer_offset = offset_meters_to_geographic_degrees(center_lat, tile_buffer_meters)
if cell_size < 10:
raise Exception('Value for cell_size must be >= 10 meters.')

center_lat = (min_lat + max_lat) / 2
lon_side_offset, lat_side_offset = offset_meters_to_geographic_degrees(center_lat, cell_size)
if buffer_size == 0:
lon_buffer_offset = 0
lat_buffer_offset = 0
else:
lon_buffer_offset, lat_buffer_offset = offset_meters_to_geographic_degrees(center_lat, buffer_size)

# Polygon Size
while lat_coord < max_lat:
Expand Down
52 changes: 18 additions & 34 deletions tests/conftest.py
View file
Original file line number Diff line number Diff line change
@@ -1,38 +1,19 @@
from city_metrix.layers import Layer

import shapely.geometry as geometry
import geopandas as gpd
from city_metrix.layers.layer import create_fishnet_grid
from geocube.api.core import make_geocube


def create_fishnet_grid(min_x, min_y, max_x, max_y, cell_size):
x, y = (min_x, min_y)
geom_array = []

# Polygon Size
while y < (max_y - 0.000001):
while x < (max_x - 0.000001):
geom = geometry.Polygon(
[
(x, y),
(x, y + cell_size),
(x + cell_size, y + cell_size),
(x + cell_size, y),
(x, y),
]
)
geom_array.append(geom)
x += cell_size
x = min_x
y += cell_size

fishnet = gpd.GeoDataFrame(geom_array, columns=["geometry"]).set_crs("EPSG:4326")
def test_create_fishnet_grid(min_x, min_y, max_x, max_y, cell_size):
# Slightly reduce aoi to avoid partial cells
reduction = 0.000001
max_y = (max_y - reduction)
max_x = (max_x - reduction)
fishnet = create_fishnet_grid(min_x, min_y, max_x, max_y, cell_size)
fishnet.drop('fishnet_geometry', axis=1, inplace=True)
return fishnet


# Test zones of a regular 0.01x0.01 grid over a 0.1x0.1 extent
ZONES = create_fishnet_grid(106.7, -6.3, 106.8, -6.2, 0.01).reset_index()
LARGE_ZONES = create_fishnet_grid(106, -7, 107, -6, 0.1).reset_index()
# Test zones of a regular 0.01x0.01 grid over a 0.1x0.1 extent by degrees
ZONES = test_create_fishnet_grid(106.7, -6.3, 106.8, -6.2, 0.01).reset_index()
LARGE_ZONES = test_create_fishnet_grid(106, -7, 107, -6, 0.1).reset_index()


class MockLayer(Layer):
Expand All @@ -54,7 +35,8 @@ class MockMaskLayer(Layer):
Simple layer where even indices are masked
"""
def get_data(self, bbox):
mask_gdf = create_fishnet_grid(*bbox, 0.01).reset_index()
cell_size_degrees = 0.01
mask_gdf = test_create_fishnet_grid(*bbox, cell_size_degrees).reset_index()
mask_gdf['index'] = mask_gdf['index'] % 2
mask = make_geocube(
vector_data=mask_gdf,
Expand All @@ -72,7 +54,8 @@ class MockGroupByLayer(Layer):
Simple categorical layer with alternating 1s and 2s
"""
def get_data(self, bbox):
group_by_gdf = create_fishnet_grid(*bbox, 0.001).reset_index()
cell_size_degrees = 0.001
group_by_gdf = test_create_fishnet_grid(*bbox, cell_size_degrees).reset_index()
group_by_gdf['index'] = (group_by_gdf['index'] % 2) + 1
group_by = make_geocube(
vector_data=group_by_gdf,
Expand Down Expand Up @@ -104,7 +87,8 @@ class MockLargeGroupByLayer(Layer):
"""

def get_data(self, bbox):
group_by_gdf = create_fishnet_grid(*bbox, 0.01).reset_index()
cell_size_degrees = 0.01
group_by_gdf = test_create_fishnet_grid(*bbox, cell_size_degrees).reset_index()
group_by_gdf['index'] = (group_by_gdf['index'] % 2) + 1
group_by = make_geocube(
vector_data=group_by_gdf,
Expand All @@ -113,4 +97,4 @@ def get_data(self, bbox):
output_crs=4326,
).index

return group_by
return group_by
Loading

0 comments on commit 052d978

Please sign in to comment.