Rasterio Error – Fixing Input Shapes Do Not Overlap Raster Issue

gdalmaskingpythonrasterio

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?

Best Answer

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://www.rfc-editor.org/rfc/rfc7946 so older GeoJSON might not be in WGS84 which rasterio might assume to be the case.

Related Question