1

I have a folder with many unprojected images (i.e. LocX_1_LocY_1.tiff to LocX_32_LocY_32.tiff) which have a naming convention that is associated with a text file that has surveyed locations of the images (i.e. LocX_1, LocY_1, UTM_Easting, UTM_Northing to LocX_32, LocY_32, 499999.795400,4480849.220620).

I wanted to ask if Python GDAL could be used to apply a coordinate system and project the images to overlay the images and view them spatially on a map with ArcGIS?

  • There are tools available including Python libraries that will write information to the EXIF header of an image file but this will not georeference your images for plotting in a GIS. Better would be to feed your images into some Structure from Motion software like Agisoft Photoscan or OpenDroneMap. This will get you a mosaic image which you could then georeference and view in ArcGIS. – GBG Mar 12 '19 at 22:00
1

The methods you are looking for are gdal.SetProjection() and gdal.SetGeoTranform(). Besides information about the projection, gdal only needs the coordinates of one corner, the size of the pixel and the number of pixels.

Working code

import gdal

name_of_old_tiff = r"C:\path\some.tif"
name_of_new_tiff = r"C:\path\some_projected.tif"
driver = gdal.GetDriverByName('GTiff')
old_tiff = gdal.Open(name_of_old_tiff)
array = old_tiff.GetRasterBand(1).ReadAsArray() #this is a numpy array

projection = """copy/paste the OGC WKT info
    for your EPSG as you find it at http://spatialreference.org
    Note there are many methods for this, but this way you avoid to import other modules.
    if you need automatic way because of many different projections you can look at
    this answer https://gis.stackexchange.com/a/193475/98816
    """

bottom_left_x_coo = 288217.9 #use units of projections ie meters
pixel_x_size = 30 #use units of projections ie meters
bottom_left_y_coo = 3638738.602#use units of projections ie meters
pixel_y_size = 30 #use units of projections ie meters
geotransform = (bottom_left_x_coo, pixel_x_size, 0.0, bottom_left_y_coo,0.0, -pixel_y_size)

dtype = gdal.GDT_Int16
#         Change data type as needed:
#                gdal.GDT_Byte
#                gdal.GDT_UInt16
#                gdal.GDT_Int16
#                gdal.GDT_UInt32
#                gdal.GDT_Int32
#                gdal.GDT_Float32
#                gdal.GDT_Float64
#                gdal.GDT_CInt16
#                gdal.GDT_CInt32
#                gdal.GDT_CFloat32
#                gdal.GDT_CFloat64
new_tiff = driver.Create(name_of_new_tiff, old_tiff.RasterXSize, old_tiff.RasterYSize, 1, dtype)
new_tiff.SetGeoTransform(geotransform)
new_tiff.SetProjection(projection)

new_tiff.GetRasterBand(1).WriteArray(array) #not saved yet
new_tiff.FlushCache() #save 
new_tiff = None #remove from memory

Now test it

test = gdal.Open(name_of_new_tiff)
projection_t = test.GetProjection()
geotransform_t = test.GetGeoTransform()

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.