EOX GitLab Instance

Commit 0b9c2e13 authored by Lubomir Doležal's avatar Lubomir Doležal
Browse files

update georeference from gsc footprint taking no-data areas in account

getting the actual places where footprint should be pointing to, this is 
only a guess though, plus order of points can not be taken as granted
parent dd7599e2
......@@ -4,11 +4,11 @@ import logging
from glob import glob
import shutil
from typing import List, Tuple
import numpy as np
from ..metadata import extract_footprint
from ..util import gdal, osr, replace_ext, get_all_data_files, ogr
logger = logging.getLogger(__name__)
......@@ -310,27 +310,72 @@ def gsc_footprint_georef(input_filename: os.PathLike, target_filename: os.PathLi
set_gcps_from_gsc_footprint(ds, {})
def set_gcps_from_gsc_footprint(ds, config):
def set_gcps_from_gsc_footprint(ds, config, noData=0):
logger.info("Trying to set GCPs from GSC footprint.")
gsc_files = get_all_data_files('extra', config, ['*GSC*.xml'])
wkt_footprint = extract_footprint(gsc_files)
if wkt_footprint is not None:
# assume rotated raster, so get bbox of this footprint to match raster file
# assuming geometry is correct
geom = ogr.CreateGeometryFromWkt(wkt_footprint)
extent = geom.GetEnvelope()
bl = (extent[2], extent[0])
br = (extent[2], extent[1])
tr = (extent[3], extent[1])
tl = (extent[3], extent[0])
gcps = gcps_from_borders(
(ds.RasterXSize, ds.RasterYSize),
(bl, br, tl, tr),
'descending',
False
)
sr = osr.SpatialReference()
sr.ImportFromEPSG(4326)
ds.SetGCPs(gcps, sr.ExportToWkt())
geotransform = gdal.GCPsToGeoTransform(ds.GetGCPs())
ds.SetGeoTransform(geotransform)
polygon = geom.GetGeometryRef(0)
if polygon.GetPointCount() == 5:
gcps = []
tl, bl, br, tr = [
(polygon.GetPoint(i)[0], polygon.GetPoint(i)[1]) for i in range(0, 4)
]
np_array = np.array(ds.GetRasterBand(1).ReadAsArray())
# find first non-nodatavalue pixel in cols/rows
x_size = ds.RasterXSize
y_size = ds.RasterYSize
found = False
# find topleft
for i, row in enumerate(np_array):
if not found:
for j, col in enumerate(row):
if col != noData:
gcps.append(gdal.GCP(tl[0], tl[1], 0, j + 0.5, i + 0.5))
found = True
break
else:
break
# find bottomright
found = False
for i, row in enumerate(np_array[::-1]):
if not found:
for j, col in enumerate(row):
if col != noData:
gcps.append(gdal.GCP(br[0], br[1], 0, j + 0.5, y_size - i - 0.5))
found = True
break
else:
break
# find bottomleft
transpose = np_array.transpose()
found = False
for i, col in enumerate(transpose):
if not found:
for j, row in enumerate(col):
if row != noData:
gcps.append(gdal.GCP(bl[0], bl[1], 0, i + 0.5, j + 0.5))
found = True
break
else:
break
found = False
# find bottomright
for i, col in enumerate(transpose[::-1]):
if not found:
for j, row in enumerate(col):
if row != noData:
gcps.append(gdal.GCP(tr[0], tr[1], 0, x_size - i - 0.5, j + 0.5))
found = True
break
else:
break
sr = osr.SpatialReference()
sr.ImportFromEPSG(4326)
ds.SetGCPs(gcps, sr.ExportToWkt())
geotransform = gdal.GCPsToGeoTransform(ds.GetGCPs())
ds.SetGeoTransform(geotransform)
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment