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:
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!):
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