EOX GitLab Instance

Commit 442dcba8 authored by Lubomir Doležal's avatar Lubomir Doležal
Browse files

attempt to at least approximately georeference using a footprint from gsc when...

attempt to at least approximately georeference using a footprint from gsc when everything fails for core09
parent 2d62e4c7
from lxml import etree
from .util import pairwise
def evaluate_xpath(root, xpath):
"""
"""
result = root.xpath(xpath, namespaces=root.nsmap)
print(xpath, result)
if result:
if isinstance(result, list):
return result[0]
......@@ -44,4 +43,32 @@ def extract_product_type_and_level(metadata_files, config):
if product_type and product_level:
break
return product_type, product_level
\ No newline at end of file
return product_type, product_level
def parse_ring(string):
raw_coords = string.split()
return [(lat, lon) for lat, lon in pairwise(raw_coords)]
def parse_gcps_gsc(elem):
interior = None
if elem is not None:
interior = parse_ring(
elem.xpath(
"gml:exterior/gml:LinearRing/gml:posList", namespaces=elem.nsmap
)[0].text.strip()
)
return interior
def extract_footprint(metadata_files, footprint_extractor: str = '//gml:target/eop:Footprint/gml:multiExtentOf/gml:MultiSurface/gml:surfaceMembers/gml:Polygon'):
footprint = None
for metadata_file in metadata_files:
with open(metadata_file) as f:
tree = etree.parse(f)
root = tree.getroot()
footprint = evaluate_xpath(root, footprint_extractor)
if footprint is not None:
break
return parse_gcps_gsc(footprint)
......@@ -6,13 +6,12 @@ from lxml import etree
import textwrap
from .util import get_all_data_files, gdal, replace_ext
from .steps.georeference import set_gcps_from_gsc_footprint
logger = logging.getLogger(__name__)
PROJECTIONS = {
"EUROPEAN": "EPSG:3035",
"CZ": "EPSG:32633",
"SK": "EPSG:32633",
"ES_29": "EPSG:32629",
"ES_30": "EPSG:32630",
}
......@@ -125,10 +124,10 @@ def attempt_to_set_projection(source_dir, target_dir, preprocess_config):
if not epsg_code:
logger.warning("Projection %s not in the mapping of projections" % projection)
# use GSC bbox as a fallback if everything failed, risking bad orientation
pass
set_gcps_from_gsc_footprint(ds, preprocess_config)
else:
ds.SetProjection(epsg_code)
del ds
del preprocess_config['stack_bands']
else:
logger.warn("metadata.xml not found. No way to fix .bil without geotransform")
\ No newline at end of file
logger.warn("metadata.xml not found. No way to fix .bil without geotransform")
......@@ -5,6 +5,7 @@ from glob import glob
import shutil
from typing import List, Tuple
from ..metadata import extract_footprint
from ..util import gdal, osr, replace_ext, get_all_data_files
......@@ -287,3 +288,28 @@ def gcps_from_borders(size: Tuple[float, float], coords: List[Tuple[float, float
])
return gcps
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]
]
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)
else:
logger.warning("Failed to set corner GCPs from GSC file, coordinates were: " % coords)
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