Невозможно обрезать изображение с помощью rasterio.mask - PullRequest
Новая
       104

Невозможно обрезать изображение с помощью rasterio.mask

0 голосов
/ 24 сентября 2019

Я пытаюсь обрезать мой файл TIFF, используя файл формы или геоджон в Python.Код для вырезания изображения: - 

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

Но я получаю сообщение об ошибке -

Файл "tryWithNewEPSG.py", строка 24, в out_image, out_transform = mask (src, geoms, crop = True) Файл "/home/ubuntu/.local/lib/python2.7/site-packages/rasterio/mask.py ", строка 181, в маске pad = pad) Файл" /home/ubuntu/.local/lib/python2.7/site-packages/rasterio/mask.py ", строка 87, в raster_geometry_mask поднять 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 Ответ

0 голосов
/ 24 сентября 2019

Проблема заключается в методе создания файла формы;используйте следующий код для создания файла формы - 


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.

Теперь попробуйте запустить код отсечения

...