Intentando recortar un raster en QGIS, "La línea de corte no es válida"

Question

Estoy usando QGIS 2.14.5-Essen y estoy tratando de recortar una imagen simplemente usando un polígono cuadrado. Me encuentro con un error que nunca he visto antes. Esto es lo que estoy tratando de hacer:

"clipper" es un shapefile que contiene una característica de polígono rectangular sobre mi área de estudio deseada.

La imagen que estoy tratando de recortar es una imagen de Reflectancia Superficial LANDSAT-8. Cuando la cargo en QGIS, ya está en el sistema de referencia espacial

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

y he comprobado que aparece en la zona correcta cargando el mapa base de OSM y confirmando que está orientado correctamente.

Mi shapefile "clipper" también está en

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

Sin embargo, cuando intento ejecutar la operación, que corresponde al comando gdal

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

Me sale este error:

No entiendo por qué no se recorta normalmente, porque ya lo he hecho antes sin problemas. No entiendo el error en absoluto y no puedo encontrar nada útil con una búsqueda en Google

Intentando recortar la imagen manualmente usando la "Extensión" en el cuadro de diálogo y arrastrando un cuadro se ejecuta la herramienta con éxito, pero obtengo una imagen que va de NaN a NaN

Answer 1

(Opción 1) Podrías probar la opción en gdalwarp de ignorar los polígonos no válidos. Esto se invoca usando "--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

(Opción 2) Dado que su polígono es un cuadrado, podría leer el shapefile en Python/OGR para obtener el área de delimitación de la capa/cuadrado, y luego utilizar las coordenadas de la capa para recortar su 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.")

Answer 2

Normalmente este mensaje de error indica que el anillo exterior del polígono se trata como anillo interior (la dirección de los puntos en el anillo no es la de las agujas del reloj).

Answer 3

Yo también me encontré con este problema y lo solucioné en cuatro pasos.

Paso 1. Si está trabajando con datos Landsat, deberá definir la ausencia de datos para cada banda utilizando la herramienta "Fill nodata" antes de continuar con el paso 2-5. Haga esto para cada banda por separado y luego fusione para el recorte. Asegúrese de seleccionar para poner cada capa en una banda separada cuando fusione utilizando la herramienta GDAL Merge.

Paso 2. Cree un ráster virtual con la capa fusionada como se muestra a continuación:

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

Paso 3. He evitado el shapefile rectangular "clipper" como en su ejemplo. Podría funcionar usando estos procedimientos, pero yo recorté mi raster usando la herramienta GDAL Clip Raster by Extent (la herramienta Clipper regular también funciona). Tenía varios rásteres - así que copié la extensión al portapapeles para repetir el proceso. Tiene que ser muy cuidadoso con sus proyecciones para que esto funcione - de hecho, la causa principal detrás del error es un problema con la proyección (ver Paso 4).

Paso 4. Verás que ocurren cosas raras si guardas el archivo como GeoTiff e intentas añadirlo al espacio de trabajo. Las capas recortadas necesitan ser reproyectadas. Yo usé el GDAL Warp(reproject) para hacer esto. Creo que puedes hacer este paso antes de recortar - el orden de los pasos 3 y 4 no debería importar.

Esto recortó la trama y mantuvo los datos. Estos pasos también resolvieron otro problema que tuve con la fusión, que funcionó a las mil maravillas.

Otra opción es trabajar con el nuevo plugin SCP, que es extremadamente fácil de usar y que probablemente querrá utilizar si está trabajando con datos Landsat. Puedes recortar los rásteres dibujando un rectángulo o utilizando un shapefile. He probado ambos usando QGIS 2.18.21 y funcionó muy bien.