Más

Replicando el resultado de gdalwarp usando enlaces gdal Python

Replicando el resultado de gdalwarp usando enlaces gdal Python


Estoy tratando de volver a proyectar / volver a muestrear con los enlaces de Python GDAL, pero obtengo resultados ligeramente diferentes en comparación con los de la utilidad de línea de comandogdalwarp.

Consulte la actualización a continuación para ver un ejemplo más corto.

Este script ilustra el enfoque de Python:

from osgeo import osr, gdal import numpy def reproject_point (point, srs, target_srs): "Reproyecta un par de coordenadas de un sistema de referencia espacial a otro." transform = osr.CoordinateTransformation (srs, target_srs) (x, y, z) = transform.TransformPoint (* point) return (x, y) def reproject_bbox (top_left, bottom_right, srs, dest_srs): x_min, y_max = top_left x_max, y_min = bottom_right esquinas = [(x_min, y_max), (x_max, y_max) , (x_max, y_min), (x_min, y_min)] projected_corners = [reproject_point (crnr, srs, dest_srs) para crnr en las esquinas] dest_top_left = (min ([crnr [0] para crnr en projected_corners]), max ([crnr [1] para crnr en esquinas_proyectadas])) dest_bottom_right = (max ([crnr [0] para crnr en esquinas_proyectadas]), min ([crnr [1] para crnr en esquinas_proyectadas])) return dest_top_left, dest_bottom_right ##### ########################################################################################################################################################################################################################################################################################## ######################### # Crear datos sintéticos gtiff_drv = gdal.GetDriverByName ('GTiff') w, h = 512, 512 raster = numpy. ceros ((w, h), dtype = numpy.uint8) raster [:: w / 10,:] = 255 raster [:, :: h / 10] = 255 top_left = (-109764, 215677) pixel_size = 45 src_srs = osr.SpatialReference () src_srs. ImportFromEPSG (3413) src_geotran = [top_left [0], pixel_size, 0, top_left [1], 0, -pixel_size] filas, cols = raster.shape src_ds = gtiff_drv.Create ('test_epsg3413.tif', columnas, filas, 1 , gdal.GDT_Byte) src_ds.SetGeoTransform (src_geotran) src_ds.SetProjection (src_srs.ExportToWkt ()) src_ds.GetRasterBand (1) .WriteArray (raster) ########################################################################################################################################################################################################################### ############################################################################################################################################################################################################################################################################################ ############ # Reproyectar a EPSG: 3573 y aumentar la muestra a 7 m dest_pixel_size = 7 dest_srs = osr.SpatialReference () dest_srs.ImportFromEPSG (3573) # Calcular nuevos límites al volver a proyectar esquinas antiguas x_min y_max = top_left bottom_right = (x_min + cols * pixel_size, y_max - rows * pixel_size) dest_top_left, dest_bottom_right = reproject_bbox (top_left, bottom_right, src_srs, dest_srs) # Hacer el conjunto de datos x_min, y_max = dest_top_left_left x tom_right new_rows = int ((x_max - x_min) / float (dest_pixel_size)) new_cols = int ((y_max - y_min) / float (dest_pixel_size)) dest_ds = gtiff_drv.Create ('test_epsg3573.tif', new_rows, g_dalcols, 1 .GDT_Byte) dest_geotran = (dest_top_left [0], dest_pixel_size, 0, dest_top_left [1], 0, -dest_pixel_size) dest_ds.SetGeoTransform (dest_geotran) dest_ds. (src_ds, dest_ds, src_srs.ExportToWkt (), dest_srs.ExportToWkt (), gdal.GRA_NearestNeighbour) dest_data = dest_ds.GetRasterBand (1) .ReadAsArray () # Cerrar conjuntos de datos src_ds = Ninguno dest_ds = Ninguno

Compare con la salida de:

gdalwarp -s_srs EPSG: 3413 -t_srs EPSG: 3573 -tr 7 7 -r cerca de GTiff test_epsg3413.tif test_epsg3573_gdalwarp.tif

Se diferencian en tamaño (en 2 filas y 1 columna), así como con algunos valores de píxeles diferentes cerca de los bordes.

Vea la superposición transparente de test_epsg3573.tif y test_epsg3573_gdalwarp.tif a continuación. Si las imágenes fueran idénticas, solo habría píxeles en blanco y negro, no grises.

Probado con Python 2.7.8, GDAL 1.11.1, Numpy 1.9.1

Actualizar:

Aquí hay un ejemplo mucho más breve. Esto parece no ser causado por el muestreo superior, ya que lo siguiente también produce resultados inconsistentes congdalwarp

from osgeo import osr, gdal import numpy # Crear datos sintéticos gtiff_drv = gdal.GetDriverByName ('GTiff') w, h = 512, 512 raster = numpy.zeros ((w, h), dtype = numpy.uint8) raster [: : w / 10,:] = 255 ráster [:, :: h / 10] = 255 top_left = (-109764, 215677) pixel_size = 45 src_srs = osr.SpatialReference () src_srs.ImportFromEPSG (3413) src_geotran = [top_left [ 0], pixel_size, 0, top_left [1], 0, -pixel_size] filas, cols = raster.shape src_ds = gtiff_drv.Create ('test_epsg3413.tif', cols, rows, 1, gdal.GDT_Byte) src_ds.SetGeoTransform ( src_geotran) src_ds.SetProjection (src_srs.ExportToWkt ()) src_ds. src_srs.ExportToWkt (), dest_srs.ExportToWkt ()) # Hacer el conjunto de datos dest dest_ds = gtiff_drv.Create ('test_epsg3573_avrt.tif', int_ds.RasterXSize, int_ds.RasterYSize, 1, gdal.GDT_dsGeoform. dest. )) dest_ds.Se tProjection (int_ds.GetProjection ()) dest_ds.GetRasterBand (1) .WriteArray (int_ds.GetRasterBand (1) .ReadAsArray ()) # Cerrar conjuntos de datos src_ds = Ninguno dest_ds = Ninguno

Y esta es la llamada gdalwarp que espero que sea la misma, pero no lo es:

gdalwarp -s_srs EPSG: 3413 -t_srs EPSG: 3573 -of GTiff test_epsg3413.tif test_epsg3573_gdalwarp.tif

La siguiente imagen muestra cada imagen binaria resultante superpuesta al 50% de transparencia. Los píxeles de color gris claro son inconsistencias entre los dos resultados.


Obtengo los mismos resultados quegdalwarpdegdal.AutoCreateWarpedVRTsi configuro el umbral de error en 0.125 para que coincida con el predeterminado (-et) en gdalwarp. Alternativamente, puede configurar-et 0.0en tu llamada agdalwarppara que coincida con el predeterminado engdal.AutoCreateWarpedVRT.

Ejemplo

Cree una referencia para comparar con:

gdalwarp -t_srs EPSG: 4326 byte.tif warp_ref.tif

Ejecute la proyección en Python (según el código de la función "warp_27 () en el conjunto de pruebas automáticas GDAL):

# Conjunto de datos de código abierto src_ds = gdal.Open ('byte.tif') # Definir SRS de destino dst_srs = osr.SpatialReference () dst_srs.ImportFromEPSG (4326) dst_wkt = dst_srs.ExportToWkt () error_threshold = 0.125 # umbral de error -> usar mismo valor que en gdalwarp resampling = gdal.GRA_NearestNeighbour # Llame a AutoCreateWarpedVRT () para obtener los valores predeterminados para las dimensiones del ráster objetivo y geotransform tmp_ds = gdal.AutoCreateWarpedVRT (src_ds, None, # src_wkt: de izquierda a un valor predeterminado -> usará el valor predeterminado de source dst_wkt, resampling, error_threshold) # Crea el ráster deformado final dst_ds = gdal.GetDriverByName ('GTiff'). CreateCopy ('warp_test.tif', tmp_ds) dst_ds = None # Comprueba que tenemos el mismo resultado que el producido por 'gdalwarp -rb -t_srs EPSG: 4326… 'ref_ds = gdal.Open (' warp_ref.tif ') ref_cs = ref_ds.GetRasterBand (1) .Checksum () ds = gdal.Open (' warp_test.tif ') cs = ds1.GetRasterBand (1) .Checksum () si cs == ref_cs: print 'éxito, coinciden' else: print "falla, no coinciden"

Ver el vídeo: Reproject, resample and clip raster data with GDAL in Python