0

I have raster files and shapefiles. I am trying to clip raster image by shapefile using gdal.Warp, however, the result is somehow disoriented.

I tried to use coordinate system parameters, but nothing worked.

Note: both shapefile and raster files are in the same coordinate system of 'EPSG:32635'.

The function for Python I'm using looks like this:

        OutTile = gdal.Warp(fname_out, image_path, format="GTiff",
                        srcSRS = 'EPSG:32635',
                        dstSRS='EPSG:32635',
                        cutlineDSName=shape_file,
                        cutlineWhere=f"id={id:d}",
                        cropToCutline=True,
                        copyMetadata = True,
                        dstNodata = 0)
        OutTile = None 

The result looks like this: enter image description here

nmtoken
  • 13,355
  • 5
  • 38
  • 87
g123456k
  • 608
  • 6
  • 19
  • related? https://gis.stackexchange.com/questions/45053/gdalwarp-cutline-along-with-shapefile – nmtoken Aug 19 '22 at 11:34
  • It is impossible to tell without test data (raster and shape) if there is an issue or not. What cropToCutline does is "Crop the extent of the target dataset to the extent of the cutline." Extent of a raster is always a bounding box filled with full pixels. If the cutline is not a north-up rectangle it is normal that pixels do not follow the boundary of the cutline in every place. However, I believe that the cutline should be covered by pixels with the cropToCutline option. The normal cutline operation selects the pixels whose majority is within the cutline polygon. – user30184 Aug 19 '22 at 11:35
  • My problem was my cutline shape wasn't perfectly rectangular. You can easily verify this by editing in QGIS (edit shape, select "Vertex Tool", right-click on a corner vertex to open table of vertices). – cm1 Feb 29 '24 at 17:58

1 Answers1

1
# Because you need to have a perfect geometry, you can get the extent of the image with the next code:

from osgeo import gdal
from osgeo.gdal import gdalconst
import geopandas as gpd

inputpath = r'input_path_raster' # Original raster

outputpath = r'output_path_raster' # Save new raster

shapefile = r'input_path_geometry' # shapefile or geojson

# We have to get the extent of the raster

bounds = gpd.read_file(shapefile)
bounds.bounds

minx = bounds.bounds['minx'][0]
miny = bounds.bounds['miny'][0]
maxx = bounds.bounds['maxx'][0]
maxy = bounds.bounds['maxy'][0]

print('minx =', minx)
print('miny =', miny)
print('maxx =', maxx)
print('maxy =', maxy)

# First option to cut the raster

cut = gdal.Warp(destNameOrDestDS = outputpath, 
                srcDSOrSrcDSTab  = inputpath,
                outputBounds     = (minx, miny, maxx, maxy),
                cropToCutline    = True,
                format           = "GTiff",
                dstSRS           = "EPSG:32635",
                outputType       = gdalconst.GDT_UInt16) # Choose the data type format

cut = None

# Second option to cut the raster

cut = gdal.Warp(destNameOrDestDS = outputpath, 
                srcDSOrSrcDSTab  = inputpath,
                outputBounds     = (minx, miny, maxx, maxy),
                cropToCutline    = True,
                copyMetadata     = True)

cut = None
Helios
  • 427
  • 2
  • 8
  • Thanks! This makes sense, however, even with the bounding box coordinates raster is not clipped exactly where are the shape borders, strange. All coordinate systems are the same! – g123456k Aug 22 '22 at 11:07
  • Did you make a change of coordinate system? – Helios Aug 26 '22 at 18:09
  • No, actually I resolved the problem by using polygon corner pixel coordinates and then I clipped the raster. I couldn't find a solution with gdal, so instead of using shapefile, I used pixel indexes. The code for that is here https://gis.stackexchange.com/a/439054/178904 – g123456k Sep 06 '22 at 07:22