HPT georeferencer

georefTest_HPT.py

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