Applying coordinate systems and projections to unprojected images

Question

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?

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