无法使用 rasterio.mask 剪辑图像

Posted 2023-02-25

【中文标题】无法使用 rasterio.mask 剪辑图像

【英文标题】：Unable to clip image using rasterio.mask

【发布时间】：2019-09-24 10:12:00

【问题描述】：

我正在尝试使用 python 中的 shape 或 geojson 文件来剪辑我的 tiff 文件。裁剪图像的代码是 -

from datetime import date
import geopandas as gpd
import rasterio
import rasterio.features
import rasterio.warp 
from shapely.geometry import MultiPolygon, Polygon
import subprocess
import matplotlib.pyplot as plt
import geopandas as gpd
from rasterio.mask import mask


nReserve = gpd.read_file('poly.shp')    
# the polygon GeoJSON geometry

nReserve_proj = nReserve.to_crs('init': 'epsg:32643')
with rasterio.open("RGBNew.tiff") as src:
    out_image, out_transform = rasterio.mask.mask(src, nReserve_proj.geometry,crop=True)
    out_meta = src.meta.copy()
    out_meta.update("driver": "GTiff",
                 "height": out_image.shape[1],
                 "width": out_image.shape[2],
                 "transform": out_transform)

with rasterio.open("RGB_masked.tif", "w", **out_meta) as dest:
    dest.write(out_image)

创建上面用到的shape文件的代码是-

import pandas as pd
from shapely.geometry import Polygon

def bbox(lat,lng, margin):
        #return (Polygon([[19.386508042365318,72.83352971076965],[19.386163944316234,72.83352971076965],[19.387256959134895,72.83352971076965],[19.38746948894201,72.83352971076965],[19.386508042365318,72.83352971076965]]))
        return (Polygon([[72.83352971076965,19.386508042365318],[72.83352971076965,19.386163944316234],[72.83352971076965,19.387256959134895],[72.83352971076965,19.38746948894201],[72.83352971076965,19.386508042365318]]))
gpd.GeoDataFrame(crs = 'init':'epsg:32643',geometry = [bbox(10,10, 0.25)]).to_file('poly.shp')

但我得到了错误-

文件“tryWithNewEPSG.py”，第 24 行，在 out_image，out_transform = mask（src，geoms，crop=True）文件“/home/ubuntu/.local/lib/python2.7/site-packages/rasterio/mask.py”， 第 181 行，在掩码中 垫=垫）文件“/home/ubuntu/.local/lib/python2.7/site-packages/rasterio/mask.py”， 第 87 行，在 raster_geometry_mask raise ValueError('输入形状不重叠光栅。') ValueError: 输入形状不重叠光栅。

我目前正在使用代码检查 epsg-

import rasterio

with rasterio.open('NDVI.tif') as src:
        print (src.crs)

并确认是一样的；我什至尝试将两者的 epsg 更改为 4326；但还是不行。

    out_image, out_transform = rasterio.mask.mask(src, nReserve_proj.geometry,crop=True)
    out_meta = src.meta.copy()
    out_meta.update("driver": "GTiff",
                 "height": out_image.shape[1],

【参考方案1】：

问题在于创建形状文件的方法；使用以下代码创建形状文件-

def create_polygon(coords):          
    ring = ogr.Geometry(ogr.wkbLinearRing)
    for coord in coords:
        ring.AddPoint(coord[0], coord[1])

    # Create polygon
    poly = ogr.Geometry(ogr.wkbPolygon)
    poly.AddGeometry(ring)
    return poly.ExportToWkt()


def write_shapefile(poly, out_shp):
    """
    https://gis.stackexchange.com/a/52708/8104
    """
    # Now convert it to a shapefile with OGR    
    driver = ogr.GetDriverByName('Esri Shapefile')
    ds = driver.CreateDataSource(out_shp)
    layer = ds.CreateLayer('', None, ogr.wkbPolygon)
    # Add one attribute
    layer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))
    defn = layer.GetLayerDefn()

    ## If there are multiple geometries, put the "for" loop here

    # Create a new feature (attribute and geometry)
    feat = ogr.Feature(defn)
    feat.SetField('id', 123)

    # Make a geometry, from Shapely object
    #geom = ogr.CreateGeometryFromWkb(poly.wkb)
    geom = ogr.CreateGeometryFromWkt(poly)
    feat.SetGeometry(geom)

    layer.CreateFeature(feat)
    feat = geom = None  # destroy these

    # Save and close everything
    ds = layer = feat = geom = None

def main(coords, out_shp):
    poly = create_polygon(coords)
    write_shapefile(poly, out_shp)

if __name__ == "__main__":
    #coords = [(73.0367660522461,18.979999927217715), (73.03436279296875,18.96019467106038), (73.05976867675781,18.96214283338193), (73.06011199,18.979999927217715),(73.0367660522461,18.979999927217715)]
    coords = [(73.97164740781432,18.527929607251075),(73.97185125569945,18.528550143773767),(73.97234478215819,18.528305998525394),(73.97215166310912,18.52775667044178),(73.97164740781432,18.527929607251075)]
    out_shp = r'polygon.shp'
    main(coords, out_shp)

在坐标中提到坐标；确保您正在剪辑的 tiff 具有 epsg:4326 您可以使用转换代码来做到这一点-

import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling

dst_crs = 'epsg:4326'

with rasterio.open('NDVIData.tiff') as src:
  transform, width, height = calculate_default_transform(
      src.crs, dst_crs, src.width, src.height, *src.bounds)
  kwargs = src.meta.copy()
  kwargs.update(
      'crs': dst_crs,
      'transform': transform,
      'width': width,
      'height': height
  )

  with rasterio.open('NDVINew.tiff', 'w', **kwargs) as dst:
      for i in range(1, src.count + 1):
          reproject(
              source=rasterio.band(src, i),
              destination=rasterio.band(dst, i),
              src_transform=src.transform,
              src_crs=src.crs,
              dst_transform=transform,
              dst_crs=dst_crs,
              resampling=Resampling.nearest)

其中 NDVIData.tiff 是我的原始 tiff，NDVINew.tiff 是我的新 tiff，带有新的 epsg 4326。

现在尝试运行剪辑代码