[GIS] How do a transform a geotiff image to use the projection of another image

osgeopython

I have two sets of glib2 data from different sources which I would like to process together.

Each image has a different projection and resolution (I don't know if resolution is the correct word to use here.)

for example, here are my two images:

image 1

and

image 2

both of these images I generated by parsing the glib2 data and writing the output to geotiff. Here's my code:

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

I have tried parsing in both sets of glib2 data and then at the point where I'm setting the projection and geotransform I use the other set of transforms than the data.

Like this:

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()

So the data used is different from the transform and projection used.

When I do this, it makes no visible difference. It seems like the change in the transform and projection have no impact.

I know I'm way off base here, because I think I need to use the same transform and projection data as I do the actual data, but I don't know how to transform the image to another projection.

How do I do this?

Best Answer

Here you have example of what you are trying to achieve I think: http://jgomezdans.github.io/gdal_notes/reprojection.html

Below is an except from the link above.

    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 )

Good luck!

Related Question