EOX GitLab Instance

Commit 6b37bc8e authored by Lubomir Doležal's avatar Lubomir Doležal
Browse files

make georeference by gsc footprint more precise, taking rotated raster in account

will still be imprecise of there is a no-data "ring" around the actual 
raster values
parent 442dcba8
......@@ -48,18 +48,36 @@ def extract_product_type_and_level(metadata_files, config):
def parse_ring(string):
raw_coords = string.split()
return [(lat, lon) for lat, lon in pairwise(raw_coords)]
return [(lon, lat) for lat, lon in pairwise(raw_coords)]
def parse_gcps_gsc(elem):
interior = None
if elem is not None:
interior = parse_ring(
def serialize_coord_list(coords):
return ','.join(
f'{x} {y}' for x, y in coords
)
def parse_polygons_gsc(elem):
interior = serialize_coord_list(
parse_ring(
elem.xpath(
"gml:exterior/gml:LinearRing/gml:posList", namespaces=elem.nsmap
)[0].text.strip()
)
return interior
)
exteriors = [
f'''({
serialize_coord_list(
parse_ring(poslist_elem.text.strip())
)
})'''
for poslist_elem in elem.xpath(
"gml:interior/gml:LinearRing/gml:posList", namespaces=elem.nsmap
)
]
return f"POLYGON(({interior}){',' if exteriors else ''}{','.join(exteriors)})"
def extract_footprint(metadata_files, footprint_extractor: str = '//gml:target/eop:Footprint/gml:multiExtentOf/gml:MultiSurface/gml:surfaceMembers/gml:Polygon'):
......@@ -71,4 +89,4 @@ def extract_footprint(metadata_files, footprint_extractor: str = '//gml:target/e
footprint = evaluate_xpath(root, footprint_extractor)
if footprint is not None:
break
return parse_gcps_gsc(footprint)
return parse_polygons_gsc(footprint)
......@@ -6,7 +6,7 @@ import shutil
from typing import List, Tuple
from ..metadata import extract_footprint
from ..util import gdal, osr, replace_ext, get_all_data_files
from ..util import gdal, osr, replace_ext, get_all_data_files, ogr
logger = logging.getLogger(__name__)
......@@ -293,12 +293,15 @@ def gcps_from_borders(size: Tuple[float, float], coords: List[Tuple[float, float
def set_gcps_from_gsc_footprint(ds, config):
logger.info("Trying to set GCPs from GSC footprint.")
gsc_files = get_all_data_files('extra', config, ['*GSC*.xml'])
coords = extract_footprint(gsc_files)
if coords and len(coords) == 5:
tl, bl, br, tr = [
[float(num) for num in coord_pair]
for coord_pair in coords[0:4]
]
wkt = extract_footprint(gsc_files)
if wkt is not None:
# assume rotated raster, so get bbox of this footprint to match raster file
geom = ogr.CreateGeometryFromWkt(wkt)
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),
......@@ -311,5 +314,3 @@ def set_gcps_from_gsc_footprint(ds, config):
ds.SetGCPs(gcps, sr.ExportToWkt())
geotransform = gdal.GCPsToGeoTransform(ds.GetGCPs())
ds.SetGeoTransform(geotransform)
else:
logger.warning("Failed to set corner GCPs from GSC file, coordinates were: " % coords)
......@@ -21,6 +21,13 @@ except ImportError:
osr.UseExceptions()
try:
from osgeo import ogr
except ImportError:
import ogr
ogr.UseExceptions()
def replace_ext(filename: os.PathLike, new_ext: str, force_dot: bool=True) -> os.PathLike:
return splitext(filename)[0] + ('' if new_ext.startswith('.') or not force_dot else '.') + new_ext
......
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