1

I'm trying to use rasterio (v1.0.13) and fiona to perform a raster clip on a geotiff using a geojson polygon. For reference, the clip works perfectly from the command line using GDAL:

gdalwarp -cutline tmp_yard.geojson  -crop_to_cutline -of GTiff -r cubic leaf_index.tif tmp.tif

enter image description here

enter image description here

But the following does not (as per rasterio docs: https://rasterio.readthedocs.io/en/latest/topics/masking-by-shapefile.html#)

with fiona.open("tmp_yard.geojson", "r") as geojson:
    features = [feature["geometry"] for feature in geojson]
    pprint.pprint(features)
    with rasterio.open("leaf_index.tif") as src:
        out_image, out_transform = rasterio.mask.mask(src, features,
                                                           crop=True)
        out_meta = src.meta.copy()

Produces error:

[{'coordinates': [[(-11680153.379982425, 4825380.211452511), (-11680101.881582342, 4825359.649286965), (-11680099.47141429, 4825337.825410188), (-11680153.677325355, 4825348.074481825), (-11680153.379982425, 4825380.211452511)]], 'type': 'Polygon'}]

ValueError: Input shapes do not overlap raster.

I tried to apply the solution here: Masking GeoTIFF file after GeoJSON through rasterio - "Input shapes do not overlap raster" to reproject my geoson, using:

import shapely
import pyproj
def project_wsg_shape_to_csr(shape, csr):
    project = lambda x, y: pyproj.transform(
        pyproj.Proj(init='epsg:4326'),
        pyproj.Proj(init=csr),
        x,
        y
      )
    return shapely.ops.transform(project, shape)

yard=project_wsg_shape_to_csr(features, 'epsg:4326')

However, this produces the error:

AttributeError: 'list' object has no attribute 'is_empty'

Anyone got any ideas?

  • Without the full traceback for your error, it's tough to say for sure where the error is happening, BUT, it sounds like your shape might be a list instead of a Shapely geometry. Could you check that your geometry is correct? – James Jan 8 at 0:56
  • Yes, I had that wrong but I couldn't figure out how to fix it. I solved the problem by giving up on reprojecting on the geoJSON side and instead reprojecting the geoTIFF to WSG84 instead before applying the clip. – metalaureate Jan 8 at 14:19
1

Your code tries to convert from EPSG:4326 to EPSG:4326. But your data is actually in EPSG:3857 (wild guess on my end). You need to use EPSG:3857 as "source" CRS.

Try this (untested):

project = lambda x, y: pyproj.transform(
    pyproj.Proj(init='epsg:3857'),
    pyproj.Proj(init='epsg:4326'),
    x,
    y
)
yard = shapely.ops.transform(project, features)

GeoJSON had no specification until https://tools.ietf.org/html/rfc7946 so older GeoJSON might not be in WGS84 which rasterio might assume to be the case.

0

In the end, reprojecting the geoTIFF using QGIS (instead of trying to unsuccessfully reproject geoJSON) to WSG84 did the trick and made the rasterio example work.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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