7

I am using QGIS 2.14.5-Essen and I am trying to simply clip an image using a square polygon. I am encountering an error that I have never seen before. This is what I am attempting to do:

enter image description here

"clipper" is a shapefile containing a rectangular polygon feature over my desired study area.

The image I am trying to clip is a LANDSAT-8 Surface Reflectance image. When I load it into QGIS, it is already in spatial reference system

+proj=utm +zone=35 +datum=WGS84 +units=m +no_defs

and I verified that it is showing up in the correct area by loading in the OSM basemap and confirming that it is oriented properly.

My "clipper" shapefile is also in

+proj=utm +zone=35 +datum=WGS84 +units=m +no_defs

However, when I try to run the operation,which corresponds to the gdal command

gdalwarp -q -cutline D:/african_towns_identifying/data/original_image/clipper.shp -tr 30.0 30.0 -of GTiff D:/african_towns_identifying/data/original_image/LC81720712016212LGN00_sr_band2.tif D:/african_towns_identifying/data/clipped_image/band2.tif

I get this error:

enter image description here

I don't understand why it just won't clip normally, because I have done this before without any problems. I don't understand the error at all and cannot find anything helpful with a Google search

Trying to clip the image manually using the "Extent" in the dialog box and dragging a box will execute the tool successfully, but I get an image ranging from NaN to NaN

1
  • 4
    Your cutline polygon seems to be invalid based an the error message about duplicate rings. Use topology check or digitize a new polygon. Capturing the geometry as WKT and pasting it into your question would be nice.
    – user30184
    Aug 7, 2016 at 8:10

3 Answers 3

7

(Option 1) You could try the option in gdalwarp to ignore invalid polygons. This is invoked using "--config GDALWARP_IGNORE_BAD_CUTLINE YES":

gdalwarp -of GTiff -co "COMPRESS=DEFLATE" -tr 30 30 -cutline VectorName -crop_to_cutline InputRaster OutputRaster --config GDALWARP_IGNORE_BAD_CUTLINE YES

(Option 2) Since your polygon is a square, you could read the shapefile in Python/OGR to get the bounding area of the layer/square, and then use the layer coordinates to clip your raster:

from osgeo import ogr, gdal
import subprocess

# Define I/O files
InputImage = "Input.tif"
OutImage = "Out.tif"
Shapefile = "SomeName.shp"`

# Open the shapefile in read-only mode and retrieve the spatial reference
Driver = ogr.GetDriverByName("ESRI Shapefile")
VectorDataset = Driver.Open(Shapefile, 0)
layer = VectorDataset.GetLayer()
feature = layer.GetNextFeature()
geom = feature.GetGeometryRef()
spatialRef = geom.GetSpatialReference()`

# Get bounding area of the layer
minX, maxX, minY, maxY = geom.GetEnvelope()`

# Clip the InputImage to the bounding area of the input feature (using Nearest Neighbour interpolation).
OutTile = gdal.Warp(OutImage, InputImage, format="GTiff", outputBounds=[minX, minY, maxX, maxY], xRes=30, yRes=30, dstSRS=spatialRef, resampleAlg=gdal.GRA_NearestNeighbour)`

# Build overviews for the OutImage
subprocess.call("gdaladdo --config COMPRESS_OVERVIEW DEFLATE "+OutImage+" 2 4 8 16 32 64", shell=True)
print("Done.")
0

Usually such error message indicate what you outer ring of polygon treats as inner ring (the direction of points in ring is not clockwise).

0

I ran into this issue as well and solved the problem in four steps.

Step 1. If you are working with Landsat data, then you should define no data for each band using the "Fill nodata" tool before proceeding through step 2-5.Fill nodata Do this for each separate band and then merge for clipping. Make sure you select to put each layer into a separate band when you merge using the GDAL Merge tool. GDAL merge multiple layers

Step 2. Create a virtual raster with the merged layer as shown in the following:

https://www.youtube.com/watch?v=YofH7eFTnNA

Step 3. I avoided the rectangular shapefile "clipper" as in your example. It might work using these procedures, but I clipped my raster using the GDAL Clip Raster by Extent tool (the regular Clipper tool works as well). I had multiple rasters - so I copied the extent to the clipboard to repeat the process. You have to be very careful with your projections to make this work - in fact, the main cause behind the error is a problem with the projection (see Step 4).

Step 4. You will find weird things happen if you save the file as a GeoTiff and try to add these to the workspace. The clipped layers need to be re-projected. I used the GDAL Warp(reproject) to do this. I think you can do this step before you clip - the order for step 3 and 4 shouldn't matter.

This clipped the raster and maintained the data. These steps also solved another issue I faced with merging, which worked like a charm.

Another option is to work with the new SCP plugin - which is extremely user friendly and you will probably want to use it if you are working with Landsat data. You can clip rasters by drawing a rectangle or by using a shapefile. I've tested both using QGIS 2.18.21 and it worked great. SCP raster clipping - panel with highlighted areas showing where to select

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.