Wednesday 26 July 2017

python - Rasterio error: Input shapes do not overlap raster but reprojection not working


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?



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


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...