Generate better cut lines for less-rectangular maps

[kylebarron]

E.g.

There was some rasterio function that had a densify_edges argument? You should be able to use that when converting back to image coordinates from geospatial coordinates to have a better curve.

[kylebarron]

You may want to use transform_geom instead of transform_bounds: https://rasterio.readthedocs.io/en/latest/api/rasterio.warp.html#rasterio.warp.transform_geom

[kylebarron]

Concise:

# Create a cutline in image coordinates
proj_bounds = Polygon.from_bounds(*transform_bounds(CRS.from_epsg(4326), profile['crs'], *bounds))
cutline = transform(lambda x,y,z=None: list(inv_geotransform * (x,y)), proj_bounds).wkt

You probably still want to use transform_geom instead of transform_bounds though.

[kylebarron]

Full example:

with rasterio.open('AK_Ruby_361345_1951_250000_geo.tif') as ds:
    profile = ds.profile
    geotransform = profile['transform']
    inv_geotransform = ~geotransform
    # Create a cutline in image coordinates
    proj_bounds = Polygon.from_bounds(*transform_bounds(CRS.from_epsg(4326), profile['crs'], *bounds))
    cutline = transform(lambda x,y,z=None: list(inv_geotransform * (x,y)), proj_bounds).wkt
    # Apply cutline
    with WarpedVRT(ds, cutline=cutline) as vrt:
        out_profile = vrt.profile
        data = vrt.read()
        out_profile['driver'] = 'GTiff'
        out_profile['nodata'] = 0
        with rasterio.open('clipped.tif', 'w', **out_profile) as dst:
            dst.write(vrt.read())