尝试使用 GeoJson、Python 和 rasterio 裁剪图像时出现问题答案

【问题标题】：Problem trying to crop image using GeoJson, Python and rasterio尝试使用 GeoJson、Python 和 rasterio 裁剪图像时出现问题
【发布时间】：2020-11-09 17:10:45
【问题描述】：

我正在尝试使用 rasterio、json、GeoJson 和 Python 从卫星图片中裁剪出一张图像。但我面临以下问题。有谁知道如何解决它？

代码

import rasterio
from rasterio.mask import mask
import json

#opens json file
myGeoms_file = open('floripa.json')
myGeoms_json = json.load(myGeoms_file)
myGeoms_file.close()

#parsing my json to a format appropriate for rasterio
coords = myGeoms_json['coordinates'][0]
parsedJson = [{
    'type': 'Polygon', 
    'coordinates': [[]]
}]
for i in coords:
  newCoord = (i[0],i[1])
  parsedJson[0]['coordinates'][0].append(newCoord)

#starting to image cropping process
with rasterio.open("CB04.tif") as src:
    out_image, out_transform = mask(src, parsedJson, crop=True)
    out_meta = src.meta.copy()

问题：

[{'type': 'Polygon', 'coordinates': [[(-48.52309226989746, -27.596847566108902), (-48.52386474609375, -27.59707576282468), (-48.526482582092285, -27.60110715979902), (-48.52652549743652, -27.601943846249203), (-48.52382183074951, -27.604453867271626), (-48.52249145507812, -27.604263716176668), (-48.52184772491455, -27.604453867271626), (-48.51961612701416, -27.607001860107363), (-48.516783714294434, -27.60441583707903), (-48.51738452911376, -27.602248094283674), (-48.51674079895019, -27.599890149923038), (-48.51515293121338, -27.594717707186042), (-48.515281677246094, -27.594033100886566), (-48.51575374603271, -27.593995067077895), (-48.51644039154053, -27.59418523598927), (-48.52038860321045, -27.59791248005127), (-48.52163314819336, -27.597532154830994), (-48.52309226989746, -27.596847566108902)]]}]

WindowError                               Traceback (most recent call last)
/usr/local/lib/python3.6/dist-packages/rasterio/mask.py in raster_geometry_mask(dataset, shapes, all_touched, invert, crop, pad, pad_width)
     83         window = geometry_window(dataset, shapes, north_up=north_up, rotated=rotated,
---> 84                                  pad_x=pad_x, pad_y=pad_y)
     85 

6 frames
WindowError: windows do not intersect

在处理上述异常的过程中，又发生了一个异常：

ValueError                                Traceback (most recent call last)
/usr/local/lib/python3.6/dist-packages/rasterio/mask.py in raster_geometry_mask(dataset, shapes, all_touched, invert, crop, pad, pad_width)
     88         # depending on value of crop
     89         if crop:
---> 90             raise ValueError('Input shapes do not overlap raster.')
     91         else:
     92             warnings.warn('shapes are outside bounds of raster. '

ValueError: Input shapes do not overlap raster.

【问题讨论】：

栅格和多边形是否在同一个投影中？

标签： python geojson crop gdal rasterio

【解决方案1】：

您基本上是在创建光栅文件的掩码版本，为您的解决方案执行以下操作：

#import rioxarray module
import rioxarray
from rasterio.crs import CRS

#open the geojson file with fiona
    with fiona.open("floripa.geojson", "r") as geojson:
      #create features
      features = [feature["geometry"] for feature in geojson]

    #open raster file with rasterio
    with rasterio.open("CB04.tif") as src:
        #clip the raster with polygon
        out_image, out_transform = rasterio.mask.mask(src, features, crop=True)
        #copy meta data of the src
        out_meta = src.meta.copy()

如果需要匹配 crs，请查看有关重新投影光栅文件的更多信息

按以下方式打开您的 .tif 文件：

raster_file = rioxarray.open_rasterio("CB04.tif", masked=True).squeeze()
检查您的 crs 是否有任何不匹配

打印（raster_file.rio.crs） #以以下格式打印 crs：EPSG:32636
打开geojson文件并检查crs

data = gpd.read_file("crops.geojson") 打印（数据.crs） #打印 crs：EPSG:XXXXXX

会有不匹配，如果你想匹配他们重新投影光栅文件CB04.tif

#open raster file ("CB04.tif") using rioxaarray

    raster_file = rioxarray.open_rasterio("CB04.tif", masked=True).squeeze()

#创建crs对象，使dst_crs = EPSG你要改＃例如。 "EPSG:4326" 如果 data.crs = EPSG:4326

    crs_wgs84 = CRS.from_string("EPSG:4326")
    raster_file = raster_file.rio.set_crs(crs_wgs84) 
#check the crs of your raster file:
print(raster_file.crs)
#prints EPSG:4326

【讨论】：