¿Cómo transformar una imagen de geotiff para usar la proyección de otra imagen?

Question

Tengo dos juegos de glib2 de datos de diferentes fuentes en las que me gustaría proceso juntos.

Cada imagen tiene una proyección diferentes y resolución (no sé si la resolución es la palabra correcta aquí.)

por ejemplo, aquí están mis dos imágenes:

y

estas dos imágenes generadas por el análisis de la glib2 de datos y la escritura de la salida a geotiff. Aquí está mi código:

data = gdal.Open(file)
driver = gdal.GetDriverByName('GTiff')

image = driver.Create(output_file, sizex, sizey, 4, gdal.GDT_Byte)

image.SetGeoTransform(data.GetGeoTransform())
image.SetProjection(data.GetProjection())

band = data.GetRasterBand(1).ReadAsArray()
shape = numpy.shape(band)

r = image.GetRasterBand(1)
g = image.GetRasterBand(2)
b = image.GetRasterBand(3)
a = image.GetRasterBand(4)

aarray = a.ReadAsArray()

#.. scale arrays

r.WriteArray(band)
r.FlushCache()
g.WriteArray(band)
g.FlushCache()
b.WriteArray(band)
b.FlushCache()
a.WriteArray(aarray)
a.FlushCache()

image = None

He tratado de analizar en ambos conjuntos de glib2 de datos y, a continuación, en el punto donde estoy ajuste de la proyección y geotransform yo uso el otro conjunto de transformaciones de los datos.

Como este:

data_1 = gdal.Open(file_1)
data_2 = gdal.Open(file_2)

image.SetGeoTransform(data_2.GetGeoTransform())
image.SetProjection(data_2.GetProjection())

band = data_1.GetRasterBand(1).ReadAsArray()

Así que los datos que utiliza es diferente de la transformación y de proyección utilizado.

Cuando hago esto, no hace ninguna diferencia visible. Parece que el cambio en la transformación y proyección no tienen ningún impacto.

Sé que estoy muy lejos de la base aquí, porque yo creo que es necesario utilizar la misma transformación y proyección de datos como puedo hacer que los datos reales, pero no sé cómo transformar la imagen a otra proyección.

¿Cómo puedo hacer esto?

Answer 1

Aquí tiene un ejemplo de lo que está tratando de lograr, creo: http://jgomezdans.github.io/gdal_notes/reprojection.html

Abajo hay una excepción del enlace de arriba.

     g = gdal.Open ( dataset )
    # Get the Geotransform vector
    geo_t = g.GetGeoTransform ()
    x_size = g.RasterXSize # Raster xsize
    y_size = g.RasterYSize # Raster ysize
    # Work out the boundaries of the new dataset in the target projection
    (ulx, uly, ulz ) = tx.TransformPoint( geo_t[0], geo_t[3])
    (lrx, lry, lrz ) = tx.TransformPoint( geo_t[0] + geo_t[1]*x_size, \
                                          geo_t[3] + geo_t[5]*y_size )
    # See how using 27700 and WGS84 introduces a z-value!
    # Now, we create an in-memory raster
    mem_drv = gdal.GetDriverByName( 'MEM' )
    # The size of the raster is given the new projection and pixel spacing
    # Using the values we calculated above. Also, setting it to store one band
    # and to use Float32 data type.
    dest = mem_drv.Create('', int((lrx - ulx)/pixel_spacing), \
            int((uly - lry)/pixel_spacing), 1, gdal.GDT_Float32)
    # Calculate the new geotransform
    new_geo = ( ulx, pixel_spacing, geo_t[2], \
                uly, geo_t[4], -pixel_spacing )
    # Set the geotransform
    dest.SetGeoTransform( new_geo )
    dest.SetProjection ( osng.ExportToWkt() )
    # Perform the projection/resampling 
    res = gdal.ReprojectImage( g, dest, \
                wgs84.ExportToWkt(), osng.ExportToWkt(), \
                gdal.GRA_Bilinear )
 

¡Buena suerte!

Answer 2

Descubrí que es necesario usar la función ReprojectImage y proporcionar tanto la fuente como la proyección de destino como argumentos. La fuente es el archivo que contiene los datos, la coincidencia es la proyección del archivo que quiero copiar.

 from osgeo import gdal,gdalconst

# Open files
data_src = gdal.Open(file_1)
data_match = gdal.Open(file_2)

# Create the result dataset
data_result = gdal.GetDriverByName('MEM').Create('', data_match.RasterXSize, data_match.RasterYSize, 1, gdalconst.GDT_Float32)

# Set the result's projection and transform to be the matching one
data_result.SetGeoTransform(data_match.GetGeoTransform())
data_result.SetProjection(data_match.GetProjection())

# reproject the data
gdal.ReprojectImage(data_src,data_result,data_src.GetProjection(),data_match.GetProjection(), gdalconst.GRA_Bilinear)
 

Ahora el resultado proyectado correctamente se almacena en data_result .