试图使用Geojson，Python和Rasterio

Question

sissue

：

[{'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.

---

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

#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

Answer 1

以以下方式打开.tif文件：raster_file= rioxarray.open_rasterio（“ cb04.tif”，masked = true）.squeeze（）

1. 检查任何不匹配的CRS

   Print（raster_file.rio.crs） ＃以一种格式来打印CRS：EPSG：32636

打开Geojson文件并检查CRS

2. data= gpd.read_file（“ orcops.geojson”） 打印（data.crs） #prints crs：epsg：xxxxxx

如果要匹配它们，将是一个不匹配的

3. #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

关注此方法，您可以使用

json

软件包用Geojson层执行剪辑。

象：掩码层和图像的坐标参考系统（CRS）必须相同。

import rasterio
from rasterio.mask import mask
import json

# Open geojson file (a file with .geojson extension is valid too)
with open('mask.json') as src:
    geojson_dict = json.load(src)
    # Store the features' geometry in a list
    features = [f["geometry"] for f in geojson_dict["features"]]

# Open the image to start cropping process
with rasterio.open("image.tif") as src:
    out_image, out_transform = mask(src, features, crop=True)
    out_meta = src.meta.copy()

Answer 2