Skip to content

Commit

Permalink
Merge pull request #253 from joshuacortez/feat/fast-square-gridding
Browse files Browse the repository at this point in the history 
Add `FastSquareGridGenerator`
  • Loading branch information
@joshuacortez
joshuacortez authored Aug 29, 2024
2 parents 54553e9 + 59664a3 commit e6e18fb
Show file tree
Hide file tree
Showing 8 changed files with 1,935 additions and 419 deletions.
33 changes: 29 additions & 4 deletions geowrangler/_modidx.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -109,8 +109,12 @@
'geowrangler/distance_zonal_stats.py'),
'geowrangler.distance_zonal_stats.create_distance_zonal_stats': ( 'distance_zonal_stats.html#create_distance_zonal_stats',
'geowrangler/distance_zonal_stats.py')},
'geowrangler.gridding_utils.polygon_fill': { 'geowrangler.gridding_utils.polygon_fill.interpolate_x': ( 'polygon_fill.html#interpolate_x',
'geowrangler.gridding_utils.polygon_fill': { 'geowrangler.gridding_utils.polygon_fill.fast_polygon_fill': ( 'polygon_fill.html#fast_polygon_fill',
'geowrangler/gridding_utils/polygon_fill.py'),
'geowrangler.gridding_utils.polygon_fill.interpolate_x': ( 'polygon_fill.html#interpolate_x',
'geowrangler/gridding_utils/polygon_fill.py'),
'geowrangler.gridding_utils.polygon_fill.polygons_to_vertices': ( 'polygon_fill.html#polygons_to_vertices',
'geowrangler/gridding_utils/polygon_fill.py'),
'geowrangler.gridding_utils.polygon_fill.scanline_fill': ( 'polygon_fill.html#scanline_fill',
'geowrangler/gridding_utils/polygon_fill.py'),
'geowrangler.gridding_utils.polygon_fill.voxel_traversal_2d': ( 'polygon_fill.html#voxel_traversal_2d',
Expand Down Expand Up @@ -140,8 +144,6 @@
'geowrangler/grids.py'),
'geowrangler.grids.FastBingTileGridGenerator._lng_to_xtile': ( 'grids.html#fastbingtilegridgenerator._lng_to_xtile',
'geowrangler/grids.py'),
'geowrangler.grids.FastBingTileGridGenerator._polygons_to_vertices': ( 'grids.html#fastbingtilegridgenerator._polygons_to_vertices',
'geowrangler/grids.py'),
'geowrangler.grids.FastBingTileGridGenerator._xtile_to_lng': ( 'grids.html#fastbingtilegridgenerator._xtile_to_lng',
'geowrangler/grids.py'),
'geowrangler.grids.FastBingTileGridGenerator._xy_to_bbox': ( 'grids.html#fastbingtilegridgenerator._xy_to_bbox',
Expand All @@ -152,6 +154,26 @@
'geowrangler/grids.py'),
'geowrangler.grids.FastBingTileGridGenerator.generate_grid': ( 'grids.html#fastbingtilegridgenerator.generate_grid',
'geowrangler/grids.py'),
'geowrangler.grids.FastSquareGridGenerator': ( 'grids.html#fastsquaregridgenerator',
'geowrangler/grids.py'),
'geowrangler.grids.FastSquareGridGenerator.__init__': ( 'grids.html#fastsquaregridgenerator.__init__',
'geowrangler/grids.py'),
'geowrangler.grids.FastSquareGridGenerator._easting_to_xtile': ( 'grids.html#fastsquaregridgenerator._easting_to_xtile',
'geowrangler/grids.py'),
'geowrangler.grids.FastSquareGridGenerator._northing_to_ytile': ( 'grids.html#fastsquaregridgenerator._northing_to_ytile',
'geowrangler/grids.py'),
'geowrangler.grids.FastSquareGridGenerator._northingeasting_to_xy': ( 'grids.html#fastsquaregridgenerator._northingeasting_to_xy',
'geowrangler/grids.py'),
'geowrangler.grids.FastSquareGridGenerator._remove_out_of_bounds_polygons': ( 'grids.html#fastsquaregridgenerator._remove_out_of_bounds_polygons',
'geowrangler/grids.py'),
'geowrangler.grids.FastSquareGridGenerator._xtile_to_easting': ( 'grids.html#fastsquaregridgenerator._xtile_to_easting',
'geowrangler/grids.py'),
'geowrangler.grids.FastSquareGridGenerator._xy_to_bbox': ( 'grids.html#fastsquaregridgenerator._xy_to_bbox',
'geowrangler/grids.py'),
'geowrangler.grids.FastSquareGridGenerator._ytile_to_northing': ( 'grids.html#fastsquaregridgenerator._ytile_to_northing',
'geowrangler/grids.py'),
'geowrangler.grids.FastSquareGridGenerator.generate_grid': ( 'grids.html#fastsquaregridgenerator.generate_grid',
'geowrangler/grids.py'),
'geowrangler.grids.H3GridGenerator': ('grids.html#h3gridgenerator', 'geowrangler/grids.py'),
'geowrangler.grids.H3GridGenerator.__init__': ( 'grids.html#h3gridgenerator.__init__',
'geowrangler/grids.py'),
Expand All @@ -176,7 +198,10 @@
'geowrangler.grids.get_intersect_partition': ( 'grids.html#get_intersect_partition',
'geowrangler/grids.py'),
'geowrangler.grids.get_parallel_intersects': ( 'grids.html#get_parallel_intersects',
'geowrangler/grids.py')},
'geowrangler/grids.py'),
'geowrangler.grids.is_aoi_within_boundary': ( 'grids.html#is_aoi_within_boundary',
'geowrangler/grids.py'),
'geowrangler.grids.setup_boundary': ('grids.html#setup_boundary', 'geowrangler/grids.py')},
'geowrangler.raster_process': { 'geowrangler.raster_process.query_window_by_gdf': ( 'raster_process.html#query_window_by_gdf',
'geowrangler/raster_process.py'),
'geowrangler.raster_process.query_window_by_polygon': ( 'raster_process.html#query_window_by_polygon',
Expand Down
107 changes: 103 additions & 4 deletions geowrangler/gridding_utils/polygon_fill.py
View file
Original file line number Diff line number Diff line change
@@ -1,13 +1,14 @@
# AUTOGENERATED! DO NOT EDIT! File to edit: ../../notebooks/15_polygon_fill.ipynb.

# %% auto 0
__all__ = ['voxel_traversal_2d', 'scanline_fill', 'voxel_traversal_scanline_fill']
__all__ = ['voxel_traversal_2d', 'scanline_fill', 'voxel_traversal_scanline_fill', 'polygons_to_vertices', 'fast_polygon_fill']

# %% ../../notebooks/15_polygon_fill.ipynb 5
from typing import List, Tuple, Set, Optional, Dict, Union

import numpy as np
import pandas as pd
import geopandas as gpd
import polars as pl

# %% ../../notebooks/15_polygon_fill.ipynb 11
Expand All @@ -18,7 +19,8 @@ def voxel_traversal_2d(
) -> List[Tuple[int, int]]:
"""
Returns all pixels between two points as inspired by Amanatides & Woo's “A Fast Voxel Traversal Algorithm For Ray Tracing”
Implementation adapted from https://www.redblobgames.com/grids/line-drawing/
Implementation adapted from https://www.redblobgames.com/grids/line-drawing/ in the supercover lines section
"""

# Setup initial conditions
Expand Down Expand Up @@ -190,8 +192,9 @@ def voxel_traversal_scanline_fill(
debug: bool = False, # if true, prints diagnostic info for both voxel traversal and scanline fill algorithms
) -> Set[Tuple[int, int]]:
"""
Returns pixels that intersect a polygon
This uses voxel traversal to fill the boundary, and scanline fill for the interior. All coordinates are assumed to be nonnegative integers
Returns pixels that intersect a polygon.
This uses voxel traversal to fill the boundary, and scanline fill for the interior. All coordinates are assumed to be integers.
"""

vertices = list(zip(vertices_df[x_col].to_list(), vertices_df[y_col].to_list()))
Expand All @@ -205,3 +208,99 @@ def voxel_traversal_scanline_fill(
polygon_pixels.update(scanline_fill(vertices, debug))

return polygon_pixels

# %% ../../notebooks/15_polygon_fill.ipynb 27
SUBPOLYGON_ID_COL = "__subpolygon_id__"
PIXEL_DTYPE = pl.Int32

# %% ../../notebooks/15_polygon_fill.ipynb 28
def polygons_to_vertices(
polys_gdf: gpd.GeoDataFrame,
unique_id_col: Optional[
str
] = None, # the ids under this column will be preserved in the output tiles
) -> pl.DataFrame:

if unique_id_col is not None:
duplicates_bool = polys_gdf[unique_id_col].duplicated()
if duplicates_bool.any():
raise ValueError(
f"""{unique_id_col} is not unique!
Found {duplicates_bool.sum():,} duplicates"""
)
polys_gdf = polys_gdf.set_index(unique_id_col)
else:
# reset index if it is not unique
if polys_gdf.index.nunique() != len(polys_gdf.index):
polys_gdf = polys_gdf.reset_index(drop=True)
unique_id_col = polys_gdf.index.name

polys_gdf = polys_gdf.explode(index_parts=True)

is_poly_bool = polys_gdf.type == "Polygon"
if not is_poly_bool.all():
raise ValueError(
f"""
All geometries should be polygons or multipolygons but found
{is_poly_bool.sum():,} after exploding the GeoDataFrame"""
)

polys_gdf.index.names = [unique_id_col, SUBPOLYGON_ID_COL]
vertices_df = polys_gdf.get_coordinates().reset_index()
vertices_df = pl.from_pandas(vertices_df)

return vertices_df

# %% ../../notebooks/15_polygon_fill.ipynb 32
def fast_polygon_fill(
vertices_df: pl.DataFrame, # integer vertices of all polygons in the AOI
unique_id_col: Optional[
str
] = None, # the ids under this column will be preserved in the output tiles
) -> pl.DataFrame:

if unique_id_col is not None:
id_cols = [SUBPOLYGON_ID_COL, unique_id_col]
has_unique_id_col = True
else:
complement_cols = ["x", "y", SUBPOLYGON_ID_COL]
unique_id_col = list(set(vertices_df.columns) - set(complement_cols))
assert len(unique_id_col) == 1
unique_id_col = unique_id_col[0]
id_cols = [SUBPOLYGON_ID_COL, unique_id_col]
has_unique_id_col = False

for col in id_cols:
assert col in vertices_df, f"{col} should be column in vertices_df"

polygon_ids = vertices_df.select(id_cols).unique(maintain_order=True).rows()

tiles_in_geom = set()
for polygon_id in polygon_ids:
subpolygon_id, unique_id = polygon_id
filter_expr = (pl.col(SUBPOLYGON_ID_COL) == subpolygon_id) & (
pl.col(unique_id_col) == unique_id
)
poly_vertices = vertices_df.filter(filter_expr)

poly_vertices = poly_vertices.unique(maintain_order=True)
_tiles_in_geom = voxel_traversal_scanline_fill(
poly_vertices, x_col="x", y_col="y"
)

if has_unique_id_col:
_tiles_in_geom = [(x, y, unique_id) for (x, y) in _tiles_in_geom]

tiles_in_geom.update(_tiles_in_geom)

schema = {"x": PIXEL_DTYPE, "y": PIXEL_DTYPE}
if has_unique_id_col:
schema[unique_id_col] = vertices_df[unique_id_col].dtype

tiles_in_geom = pl.from_records(
data=list(tiles_in_geom),
orient="row",
schema=schema,
)

return tiles_in_geom
Loading

0 comments on commit e6e18fb

Please sign in to comment.