5

I am trying to create a stacked tiff image out of 6 tiff images. Each one of them represents a different band.

When I try to do that, I get an error that the photos are not geo-referenced, which is true, because the corner coordinates of the photo are only known after calculations.

My question is, how can I geo-reference these photos in order to create a stacked tiff image, given that I know the corner coordinates along with the rest of metadata (image resolution, altitude etc.)?

This is the code I (plan to) use to create my stacked tiff image:

my_bands = "my_tif_photo_band*[1-6]*.tif"
stack_band_paths = glob(my_bands)
stack_band_paths.sort()

arr_st, meta = earthpy.spatial.stack(stack_band_paths, nodata=-9999)

I guess it can be done with gdal but didn't find an answer yet.

2
  • Read gdal.org/programs/gdal_translate.html#gdal-translate and gdal.org/programs/gdal_edit.html first. Do you need to give four corner coordinates or is two corners enough for your images?
    – user30184
    Nov 30, 2020 at 16:43
  • If 2 corner coordinates are enough to place the image on a map, then I only need 2. I've looked the links you provided but there are so much options that I'm not sure how I should construct the command that applies the geo-reference to the tif photo. Could you provide an example? Nov 30, 2020 at 16:59

1 Answer 1

6

Using GDAL utilities

As per the comments you may be able to use gdal_translate assuming you know the upper left and bottom right coordinates:

gdal_translate -of GTiff -a_srs 'EPSG:<4326>' -a_ullr <upper left x> <upper left y> <lower right x> <lower right y> <img_path> <dest_img>

The -a_srs allows you to set the spatial reference of the image using either the epsg code or a proj4 string and the -a_ullr dictates the bounds of your image and it will get made into a GeoTiff format (-of option)

Using Python

To georeference an image, you need to know the geotransform and the projection of the image.

The [geotransform][1] will be a tuple of length 6 consisting of the upper left x coordinate (ulx), resolution in the x direction (xres), rotation in the x direction (xrot), upper left y coordinate (uly), rotation in y axis (yrot), resolution in the y direction (yres).

If you know the epsg code of the final image you can also determine the projection information.

You then need to open the dataset using gdal and update this information like so:

from osgeo import gdal, osr
gdal.UseExceptions()
osr.UseExceptions()

img_path = '/path/to/tif'

#input values for geotransform
ulx = 25
uly = 60
xres = 0.5 
yres = -0.5
xrot = 0
yrot = 0
geotransform = (ulx, xres, xrot, uly, yrot, yres)


#find projection
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)

#update image georeference
ds = gdal.Open(img_path, gdal.GA_Update)
ds.SetGeoTransform(geotransform)
ds.SetProjection(srs.ExportToWkt())
ds.FlushCache()
del ds

#do this for each image and then stack

I've just assumed an WGS84 reference but replace with the relevant details for your images. [1]: https://gdal.org/tutorials/geotransforms_tut.html

3
  • Thank you a lot I appreciate it. I've followd the python approach and it works. I have one problem though. When I stack the images and then plot the RGB values, I get an image with an incorrect rotation. It needs to be rotated 180 degrees to the right. I have tried setting the xrot or yrot to 180 but nothing changed. Do you have any idea what might be the issue? Dec 1, 2020 at 12:36
  • Just to clarify do you mean, the image is flipped horizontally? If so it might be due to the y resolution in the geotransform. GDAL reads an image from the top left but as the positive y direction is towards the north , your yres should be negative? If this isn't the case, would you mind posting the geotransform and your result please? Dec 1, 2020 at 18:55
  • Sorry, my bad. It was correct all along. Thank you again :) Dec 2, 2020 at 10:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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