I wish to convert this Python code using a gdal command line tool to using the gdal python library. I am using it to convert some PNG satellite images I am downloading to georeferenced TIFFs. I manually specify the bounding box of the images when I download them (in EPSG:3857), so I know all the info about the images, I just wish to store all the info in a single georeferenced TIFF file.
My current solution is to use gdal_translate called by subprocess but apparently this can all be done natively in python, so I would prefer to do that rather than rely on this hacked-together workaround.
from osgeo import gdal
import subprocess
import shlex
# ll_x etc... are all corner coords I calculate before hand
subprocess.run(shlex.split("gdal_translate -of GTiff -a_ullr {} {} {} {} -a_srs {} {} {}".format(
ll_x, ur_y, ur_x, ll_y,
"EPSG:3857",
"input_sat_image.png",
"referenced_sat.tif")))
I have seen individual questions that cover pieces of this, but I am struggling to put them all together.
1 Answer 1
I might be wrong, but maybe you can try this way :
from osgeo import gdal
import osr
inputImage = YOURNAME
outputImage = YOURNAME
dataset = gdal.Open(inputImage)
I = dataset.ReadAsArray(0,0,dataset.RasterXSize,dataset.RasterYSize)
outdataset = gdal.GetDriverByName('GTiff')
output_SRS = osr.SpatialReference()
output_SRS.ImportFromEPSG(3857)
outdataset = outdataset.Create(outputImage,dataset.RasterXSize,dataset.RasterYSize,I.shape[0],gdal_Float32)
for nb_band in range(I.shape[0]):
outdataset.GetRasterBand(nb_band+1).WriteArray(I[nb_band,:,:])
gcp_list = []
gcp_list.append(gdal.GCP(lon,lat,alt,col,row)) #Repeat this step for as many GCP as you want (there is a limit but I don't know the number of GCPs)
outdataset.SetProjection(srs.ExportToWkt())
wkt = outdataset.GetProjection()
outdataset.SetGCPs(gcp_list,wkt)
outdataset = None
It worked for me, but it was in a specific case (and I did not need a very precise result...). I think you might have a problem with the sizes, but I actually don't know how to solve it.