Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Generate better cut lines for less-rectangular maps #10

Open
kylebarron opened this issue May 9, 2020 · 3 comments
Open

Generate better cut lines for less-rectangular maps #10

kylebarron opened this issue May 9, 2020 · 3 comments

Comments

@kylebarron
Copy link
Owner

kylebarron commented May 9, 2020

E.g.
image

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
Copy link
Owner Author

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
Copy link
Owner Author

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
Copy link
Owner Author

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())

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant