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


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)


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


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)


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?


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!


I discovered it's necessary to use the function ReprojectImage and supply both the source and dest projection as arguments. Source is the file containing the data, match is the file projection I want to copy.

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

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

Now the correctly projected result is stored in data_result.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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