无法使用剪辑图像光栅.mas

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)

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

^{pr2}$

但我得到了一个错误-

File "tryWithNewEPSG.py", line 24, in out_image, out_transform = mask(src, geoms, crop=True) File "/home/ubuntu/.local/lib/python2.7/site-packages/rasterio/mask.py", line 181, in mask pad=pad) File "/home/ubuntu/.local/lib/python2.7/site-packages/rasterio/mask.py", line 87, in raster_geometry_mask raise ValueError('Input shapes do not overlap raster.') ValueError: Input shapes do not overlap raster.

我正在用代码检查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)

在坐标中提及坐标；确保要剪裁的tiff具有爱普生：4326您可以使用转换代码-

^{pr2}$

在哪里NDVIData.tiff文件是我最初的口角蒂芙是我和新epsg 4326的新口角。在

现在尝试运行剪辑代码