Skip to content

Instantly share code, notes, and snippets.

@rukku
Last active May 13, 2021 08:39
Show Gist options
  • Save rukku/41ce00a1cc7653b6670cbb6bbfb6d7f8 to your computer and use it in GitHub Desktop.
Save rukku/41ce00a1cc7653b6670cbb6bbfb6d7f8 to your computer and use it in GitHub Desktop.
HPT georeferencer
from osgeo import gdal, osr
import csv
import os
##
gdal.UseExceptions()
testPoints = "D2_HPT_2021-05-01T061835085_p1.points"
fileName = os.path.basename(testPoints)[0:27]
if str(os.path.basename(testPoints)).find("p2") >= 1:
polyOrder = "2"
if str(os.path.basename(testPoints)).find("p1") >= 1:
polyOrder = "1"
#polyOrder = os.path.basename(pointsFiles[i])[-8:-7]
srcDN = 'D2_HPT_2021-05-01T061835085_DN.tif'
gcpFile = testPoints#pointsFiles[i]
gcpL = []
# for r in range(row):
# gcpL.append(gdal.GCP(gcps.iloc[r,0], gcps.iloc[r,1], 0, gcps.iloc[r,2], gcps.iloc[r,3]*-1)) # multiplied to -1 so that it will not flip the image
with open(testPoints) as csv_file:
csv_reader = csv.reader(csv_file, delimiter=',')
next(csv_reader)
for row in csv_reader:
offset = 1000
mapx, mapy, pixelx, pixely, *dummy = map(float, row)
pixelx = pixelx-offset
pixely = -1*(pixely-offset)
gcpL.append(gdal.GCP(mapx, mapy,0, pixelx, pixely))
# Get spatial reference
testFile = gdal.Open(srcDN)
output_srs = osr.SpatialReference()
output_srs.ImportFromEPSG(4326)
# Translate
tOptions = gdal.TranslateOptions(format="GTiff", resampleAlg= "near", GCPs=gcpL, outputSRS=output_srs, srcWin = [0, 0, 634, 466])
ds = gdal.Translate(destName = 'DN_translate.tif', srcDS = gdal.Open(srcDN), options = tOptions)
ds = None
wOptions = gdal.WarpOptions(format="GTiff", resampleAlg= "near", polynomialOrder=int(polyOrder), srcSRS=output_srs, dstSRS=output_srs)
ds = gdal.Warp(destNameOrDestDS = 'DN_p1_csv_minus1_map_offset.tif', options=wOptions, srcDSOrSrcDSTab = 'DN_translate.tif')
ds = None
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment