1

My apologies if this question has been asked before; but I could not find my specific question answered elsewhere.

I define a bounding box with bounds (xmin,xmax, ymin, ymax) in a rotated lon/lat system. If I'd supersample the points along this bounding box and transform the coordinates to a regular WGS84 system, it becomes clear (as you would expect for a rotated grid) that the bounding box ceases to be a box, and instead takes on a curved form in the other CRS:

rotated grid

Now my problem is the following: I want to mask (e.g., with rasterio) a dataset (given in regular lon/lat coordinates) with my given bounding box (given in rotated lon/lat coordinates).

A simple, but wrong, solution is to transform the bounding box coordinates to the regular lon/lat coordinates; as rasterio will then assume a straight line between the points of the polygon, i.e., it will mask following the red lines in the image below. So, the following is not the desired behavior (corner points are preserved correctly, edges are straight but should be curved!):

rotated grid wrong cut

One solution is to reproject my entire dataset into the rotated coordinate system. This is, however, not really the cheapest operation (for something I'll have to do many times over, and want to make reasonably interactive). Another solution is to do as written above, i.e., supersample the points along the bounding box, transform each of those points to the other CRS, and mask along the supersampled bounding box points. This can also get quite expensive, and it's hard to define when the curved cells are appropriately captured by the supersampling.

So I wonder if another clean solution exists.

BTW, the standard cropping/masking code is this

import rasterio
from rasterio.mask import mask

IMAGE_path = '....tif'
POL = ...
with rasterio.open(IMAGE_path) as src:
    cropped_image, _ = mask(src, 
                            POL, 
                            nodata=0,
                            crop=True,
                            all_touched=True)
    return cropped_image

1 Answer 1

0

I believe you are looking for transform bounds. This can be found in rasterio and pyproj. You need to use this to transform the bounds to the target projection.

https://rasterio.readthedocs.io/en/latest/api/rasterio.warp.html#rasterio.warp.transform_bounds

https://pyproj4.github.io/pyproj/stable/api/transformer.html#pyproj.transformer.Transformer.transform_bounds

6
  • It's not quite what I want; running, for example, print( rasterio.warp.transform_bounds('epsg:3035', 'epsg:4326', 4070757, 2964068, 4794123,3142802, densify_pts=21) ) simply gives me back another 4 corner-points for left/bottom/right/top. It, unfortunately, doesn't provide a 'supersampled' curved polygon to use for the mask.
    – Erik
    Oct 13, 2021 at 13:14
  • Right, it works as documented Transforms bounds from src_crs to dst_crs, optionally densifying the edges (to account for nonlinear transformations along these edges) and extracting the outermost bounds. You must take the bbox, densify, re-project the densified polygon, keep it with all vertices, and use that for masking. Ogr2ogr with -segmentize could probably be used for testing the polygon creation part gdal.org/programs/ogr2ogr.html.
    – user30184
    Oct 13, 2021 at 13:32
  • @user30184: the whole approach of a supersampled polygon is actually not quite what I want (as given in the original question), because it will come with a significant overhead (a lot of extra coordinates to transform, probably with a lot of wasted work...). In this question, I was hoping to find a more efficient approach to warp the vertices/edges of a polygon mask for use with rasterio, but perhaps that simply doesn't exist. Thanks for the heads-up about segmentize, I will have a look.
    – Erik
    Oct 13, 2021 at 13:50
  • I would not be afraid of transforming a few more coordinates. Even the densified bbox polygon is super easy as a GIS feature. And rasterio.warp.transform_bounds is doing just that, plus it is throwing away most of the transformed coordinates. Do you experience that it is using more time than a few milliseconds? Actually transform bounds does not select four corner points because for example in your example the max northing is taken from the middle of the line.
    – user30184
    Oct 13, 2021 at 14:31
  • @user30184: I'm actually transforming 1+ million cells (which have at least 4 boundary points each), for which it takes ~20 minutes to compute all masks. The problem is that the fastest Python implementation transforms all coordinates upfront (before masking), which thus has a serious(!) memory cost if I oversample the polygon. This is why I thought "is there a way to get away with my original 4 corner points, and let the GIS system transform the full Polygon shape", but I suppose there isn't! Thanks though, that's also relevant to know! :-)
    – Erik
    Oct 13, 2021 at 17:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.