Gridding tutorial boundary gotchas

[joshuacortez]

In the gridding tutorial, the bounding box was directly computed from the epsg:4326 admin bounds

aoi_total_bounds = region3_gdf.total_bounds

These bounds were then passed to SquareGridGenerator

grid_generator1k = grids.SquareGridGenerator(1_000, boundary=aoi_total_bounds)

In my own use case, I tried to make 100x100m grid cells across the entire Indonesian coastline. In the same way as the tutorial, I computed the bounding box on the epsg:4326 geodataframe, and passed them to SquareGridGenerator.

However I found that there was a significant number of missing grid cells on the northwestern side of Indonesia.

There are supposed to be tiles along Aceh's coast. The image below shows what should have been there. (tiles are smaller since they're from a separate run)

After some investigation, I found that it had to do with the bounding box. I should have computed the bounding box on the projected admin boundaries, not the epsg:4326 admin boundaries. Even if the proj_crs is specified in SquareGridGenerator, we’ll still have this error because the bounding box (computed from epsg:4326) was already distorted to begin with.

I resolved the issue by generating the bounding box from the projected admin bounds instead of the epsg:4326 admin bounds, and used SquareGridBoundary as follows

boundary = grids.SquareGridBoundary(*proj_admin_bounds_bbox)

After passing this boundary to SquareGridGenerator, tiles now show up along Aceh's coast.

So technically this isn't a bug with geowrangler itself but it's definitely something to watch out for. Hopefully other people won't make the same mistake. I was thinking if we should be explicit about this in the tutorial. We could:

• Use the projected bbox by default in the tutorial, so people who follow along will use it as best practice. (leaning on this one)
• Or we could keep the same flow but with a warning to use the projected bbox especially for bigger areas.

Let me know what you think! :)

[butchtm]

@joshuacortez -- just clarifying:

in the generate_grid method:

• if you don't pass in a boundary arg when the SquareGridGenerator was created -- the boundary is created from the reprojected total bounds of the gdf you pass in the generate_grid method. -- this is not the use case for your example above

• if a SquareGridBoundary as the boundary arg when SquareGridGenerator was created, it is assumed the user already did the reprojecting and computation of the reprojected boundaries. -- this is what you did to fix the problem, right?

• if an array ([minx,miny,maxx,maxy]) (e.g. aoi.total_bounds) is passed, the boundaries are reprojected and recomputed -- is this not the behavior you want? (If it didn't then this might be a bug?)

   ...
    reprojected_gdf = gdf.to_crs(self.grid_projection)
    if self.boundary is None:
        boundary = SquareGridBoundary(*reprojected_gdf.total_bounds)
    elif isinstance(self.boundary, SquareGridBoundary):
        boundary = self.boundary
    else:
        transformer = Transformer.from_crs(gdf.crs, reprojected_gdf.crs, always_xy=True)
        x_min, y_min = transformer.transform(self.boundary[0], self.boundary[1])
        x_max, y_max = transformer.transform(self.boundary[2], self.boundary[3])
        boundary = SquareGridBoundary(x_min, y_min, x_max, y_max)

[jtmiclat]

what i imagined

>>> region3_gdf = load_gdf()
>>> grid_generator1k = grids.SquareGridGenerator(1_000, boundary=region3.total_bounds)
>>> grid_generator1k.generate_grid(region3_gdf.subset(1))
Warning: The bounds of aoi after reprojection is outside the boundary after reprojection. There might be missing grid cells due to this. Trying adding a buffer to the boundary or reproject aoi before getting the boundary. See geowrangler.thinkingmachin.es/grids#repojectioning-before-gridding 



... default projection is epsg:3857
>>> reprojected_gdf = region3_gdf.to_crs("EPSG:3857")
>>> grid_generator1k = grids.SquareGridGenerator(1_000, boundary=reprojected_gdf.total_bounds)
>>> grid_generator1k.generate_grid(reprojected_gdf.subset(1))



... if with custom crs
>>> reprojected_gdf = region3_gdf.to_crs(custom_crs)
>>> grid_generator1k = grids.SquareGridGenerator(1_000, boundary=reprojected_gdf.total_bounds, grid_projection=custom_crs)
>>> grid_generator1k.generate_grid(reprojected_gdf.subset(1))

[joshuacortez]

Ok got it, so no need to pass in custom_crs to generate_grid

[joshuacortez]

Just to ask, did we get to implement this? :D

[jtmiclat]

We haven't! But we created an issue to track it: #147

[joshuacortez]

cool cool cool

[joshuacortez]

Is addressed by #253

1. SquareGridGenerator and FastSquareGridGenerator should provide a warning if the boundary does not fully enclose the AOI
2. Added markdown text documentation in tutorial.grids.ipynb that explains using boundary of reprojected AOI instead of original AOI