Commit e08e16e9 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Fixed issue #75 (Black border around L2B products).

parent 62a9f03c
......@@ -15,9 +15,9 @@ except ImportError:
from geoarray import GeoArray
from py_tools_ds.geo.coord_calc import calc_FullDataset_corner_positions
from py_tools_ds.geo.coord_trafo import pixelToLatLon
from py_tools_ds.geo.coord_trafo import pixelToLatLon, imYX2mapYX
from py_tools_ds.geo.map_info import mapinfo2geotransform
from py_tools_ds.geo.projection import EPSG2WKT
from py_tools_ds.geo.projection import EPSG2WKT, isProjectedOrGeographic
from ..options.config import GMS_config as CFG
from . import geoprocessing as GEOP
......@@ -536,7 +536,7 @@ class L1A_object(GMS_object):
rows, cols = self.shape_fullArr[:2]
CorPosXY = [(0, 0), (cols, 0), (0, rows), (cols, rows)]
gt = self.mask_nodata.geotransform
# using EPSG code ensures that excactly the same WKT code is used in case of in-mem and disk serialization
# using EPSG code ensures that exactly the same WKT code is used in case of in-mem and disk serialization
prj = EPSG2WKT(self.mask_nodata.epsg)
CorLatLon = pixelToLatLon(CorPosXY, geotransform=gt, projection=prj)
self.corner_lonlat = [tuple(reversed(i)) for i in CorLatLon]
......@@ -561,6 +561,12 @@ class L1A_object(GMS_object):
trueCorLatLon = pixelToLatLon(trueCorPosXY, geotransform=gt, projection=prj)
self.trueDataCornerLonLat = [tuple(reversed(i)) for i in trueCorLatLon]
# set true data corner positions (UTM coordinates)
# calculate trueDataCornerUTM
if isProjectedOrGeographic(prj) == 'projected':
data_corners_utmYX = [imYX2mapYX(imYX, gt=gt) for imYX in self.trueDataCornerPos]
self.trueDataCornerUTM = [(YX[1], YX[0]) for YX in data_corners_utmYX]
# set full scene corner positions (image and lonlat coordinates)
if is_dataset_provided_as_fullScene(self.GMS_identifier):
assert len(self.trueDataCornerLonLat) == 4, \
......@@ -569,7 +575,9 @@ class L1A_object(GMS_object):
% (self.scene_ID, self.entity_ID, len(self.trueDataCornerLonLat))
self.fullSceneCornerPos = self.trueDataCornerPos
self.fullSceneCornerLonLat = self.trueDataCornerLonLat
else:
if re.search('AVNIR', self.sensor):
self.fullSceneCornerPos = calc_FullDataset_corner_positions(self.mask_nodata, algorithm='numpy',
assert_four_corners=False)
......@@ -577,14 +585,18 @@ class L1A_object(GMS_object):
trueCorPosXY = [tuple(reversed(i)) for i in self.trueDataCornerPos]
trueCorLatLon = pixelToLatLon(trueCorPosXY, geotransform=gt, projection=prj)
self.fullSceneCornerLonLat = [tuple(reversed(i)) for i in trueCorLatLon]
else: # RapidEye or Sentinel-2 data
else:
# RapidEye or Sentinel-2 data
if re.search('Sentinel-2', self.satellite):
# get fullScene corner coordinates by database query
# -> calculate footprints for all granules of the same S2 datatake
# -> merge them and calculate overall corner positions
self.fullSceneCornerPos = [tuple(reversed(i)) for i in CorPosXY] # outer corner positions
self.fullSceneCornerLonLat = self.corner_lonlat # outer corner positions
else: # RapidEye
else:
# RapidEye
raise NotImplementedError # FIXME
def calc_center_AcqTime(self):
......
......@@ -8,9 +8,9 @@ import logging
import numpy as np
from geoarray import GeoArray, NoDataMask, CloudMask
from py_tools_ds.geo.projection import WKT2EPSG, get_UTMzone
from py_tools_ds.geo.projection import WKT2EPSG, get_UTMzone, isProjectedOrGeographic
from py_tools_ds.geo.coord_calc import calc_FullDataset_corner_positions
from py_tools_ds.geo.coord_trafo import pixelToLatLon, pixelToMapYX, imXY2mapXY
from py_tools_ds.geo.coord_trafo import pixelToLatLon, pixelToMapYX, imXY2mapXY, imYX2mapYX
from py_tools_ds.geo.map_info import geotransform2mapinfo, mapinfo2geotransform
from py_tools_ds.numeric.array import get_array_tilebounds
......@@ -514,10 +514,12 @@ class Dataset(object):
sub_GMS_obj.trueDataCornerLonLat = [(YX[1], YX[0]) for YX in data_corners_LatLon]
# calculate trueDataCornerUTM
data_corners_utmYX = pixelToMapYX([ULxy, URxy, LLxy, LRxy],
geotransform=subArr.gt,
projection=subArr.prj) # FIXME asserts gt in UTM coordinates
sub_GMS_obj.trueDataCornerUTM = [(YX[1], YX[0]) for YX in data_corners_utmYX]
if isProjectedOrGeographic(subArr.prj) == 'projected':
data_corners_utmYX = [imYX2mapYX(imYX, gt=subArr.gt) for imYX in corners_imYX]
sub_GMS_obj.trueDataCornerUTM = [(YX[1], YX[0]) for YX in data_corners_utmYX]
else:
self.logger.error("Error while setting 'trueDataCornerUTM' due to unexpected projection of subArr: %s."
% subArr.prj)
return sub_GMS_obj
......
......@@ -163,10 +163,16 @@ def L2A_map(L1C_objs, block_size=None, return_tiles=True):
# initialize L2A objects
L2A_objs = [L2A_P.L2A_object(L1C_obj) for L1C_obj in L1C_objs]
# get common extent (NOTE: using lon/lat coordinates here would cause black borders around the UTM image
# because get_common_extent() just uses the envelop without regard to the projection
clip_extent = \
get_common_extent([obj.trueDataCornerUTM for obj in L2A_objs]) \
if len(L2A_objs) > 1 else L2A_objs[0].trueDataCornerUTM
# correct geometric displacements and resample to target grid
common_extent = get_common_extent([obj.trueDataCornerLonLat for obj in L2A_objs]) if len(L2A_objs) > 1 else None
for L2A_obj in L2A_objs:
L2A_obj.correct_spatial_shifts(cliptoextent=CFG.clip_to_extent, clipextent=common_extent, clipextent_prj=4326)
L2A_obj.correct_spatial_shifts(cliptoextent=CFG.clip_to_extent,
clipextent=clip_extent, clipextent_prj=L2A_obj.arr.prj)
# merge multiple subsystems belonging to the same scene ID to a single GMS object
L2A_obj = L2A_P.L2A_object.from_sensor_subsystems(L2A_objs) if len(L2A_objs) > 1 else L2A_objs[0]
......
......@@ -203,6 +203,8 @@ class BaseTestCases:
# cls.PC.config.spathomo_estimate_accuracy = True
# cls.PC.config.ac_estimate_accuracy = True # FIXME
# cls.PC.config.spechomo_estimate_accuracy = True # FIXME
# cls.PC.config.exec_L1CP = [1, 1, 0]
# cls.PC.config.exec_2ACP = [1, 1, 0]
def test_run_all_processors(self):
self.PC.run_all_processors()
......
Markdown is supported
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