无法使用rasterio.mask剪辑图像

Question

我正在尝试使用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)

创建上面使用的形状文件的代码是-

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')

但我错了-

"/home/ubuntu/.local/lib/python2.7/site-packages/rasterio/mask.py"，文件"tryWithNewEPSG.py"，第24行，在out_image中，out_transform =掩码(src，geoms，crop=True)文件

第181行，在掩码pad=pad中，在文件out_transform第87行，在raster_geometry_mask raise (‘输入形状不重叠光栅’)中。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],

Answer 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)

在这些和弦中，请提到这些和弦；确保您正在裁剪的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。

现在尝试运行裁剪代码。