0

I am trying to clip my tiff file either using a shape or geojson file in python. The code for clipping the image is -

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)

The code to create the above used shape file is -

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

But I am getting the error-

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.

I am currently checking the epsg using the code-

import rasterio

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

and have confirmed it is the same; I have even tried by changing the epsg of both to 4326; but still did not work.

    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],
Smita Ahinave
  • 1,901
  • 7
  • 23
  • 42

1 Answers1

0

The problem is the method in which the shape file was created; use the following code to create the shape file-


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)

In the coords mention the coords; make sure that the tiff you are clipping has epsg:4326 you can do it using the convert code-

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)

where NDVIData.tiff is my original tiff and NDVINew.tiff is my new tiff with new epsg 4326.

Now try running the clipping code

Brian
  • 23
  • 4