Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Simplified datamask for ORB product #122

Merged
merged 7 commits into from
Sep 1, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 6 additions & 5 deletions S1_NRB/metadata/mapping.py
Original file line number Diff line number Diff line change
Expand Up @@ -117,11 +117,12 @@
'unit': None,
'role': 'data-mask',
'title': 'Data Mask Image',
'values': {0: 'not layover, nor shadow',
1: 'layover',
2: 'shadow',
3: 'layover and shadow',
4: 'ocean water'}},
'allowed': ['not layover, nor shadow',
'layover',
'shadow',
'layover and shadow',
'ocean water',
'land']},
'-ei.tif': {'type': 'Angle',
'unit': 'deg',
'role': 'ellipsoid-incidence-angle',
Expand Down
4 changes: 2 additions & 2 deletions S1_NRB/metadata/stac.py
Original file line number Diff line number Diff line change
Expand Up @@ -509,9 +509,9 @@ def _asset_handle_raster_ext(stac_asset, nodata, key=None, meta=None, asset=None
with Raster(asset) as dm_ras:
band_descr = [dm_ras.raster.GetRasterBand(band).GetDescription() for band in
range(1, dm_ras.bands + 1)]
samples = {k: v for k, v in ASSET_MAP[key]['values'].items() if v in band_descr}
samples = [x for x in band_descr if x in ASSET_MAP[key]['allowed']]
bands = []
for sample in samples.values():
for sample in samples:
band = RasterBand.create(nodata=nodata,
data_type=DataType.UINT8,
unit=ASSET_MAP[key]['unit'])
Expand Down
15 changes: 6 additions & 9 deletions S1_NRB/metadata/xml.py
Original file line number Diff line number Diff line change
Expand Up @@ -335,15 +335,12 @@ def product_xml(meta, target, assets, nsmap, ard_ns, exist_ok=False):
with Raster(asset) as dm_ras:
band_descr = [dm_ras.raster.GetRasterBand(band).GetDescription() for band in
range(1, dm_ras.bands + 1)]
if 1 < len(band_descr) < len(ASSET_MAP[key]['values']):
samples = {key: val for key, val in ASSET_MAP[key]['values'].items() if val in band_descr}
for i, sample_val in enumerate(samples.values()):
bitValue = etree.SubElement(productInformation, _nsc('_:bitValue', nsmap, ard_ns=ard_ns),
attrib={'band': str(i + 1),
'name': sample_val})
bitValue.text = '1'
else:
raise RuntimeError('{} contains an unexpected number of bands!'.format(asset))
samples = [x for x in band_descr if x in ASSET_MAP[key]['allowed']]
for i, sample in enumerate(samples):
bitValue = etree.SubElement(productInformation, _nsc('_:bitValue', nsmap, ard_ns=ard_ns),
attrib={'band': str(i + 1),
'name': sample})
bitValue.text = '1'
else: # key == '-id.tif'
src_list = list(meta['source'].keys())
src_target = [os.path.basename(meta['source'][src]['filename']).replace('.SAFE',
Expand Down
45 changes: 22 additions & 23 deletions S1_NRB/nrb.py
Original file line number Diff line number Diff line change
Expand Up @@ -231,7 +231,7 @@ def format(config, scenes, datadir, outdir, tile, extent, epsg, wbm=None,
create_data_mask(outname=dm_path, datasets=datasets_sar, extent=extent, epsg=epsg,
driver=driver, creation_opt=write_options['dm'],
overviews=overviews, overview_resampling=ovr_resampling,
wbm=wbm, dst_nodata=dst_nodata_byte)
wbm=wbm, dst_nodata=dst_nodata_byte, orb=orb)
datasets_ard['dm'] = dm_path

# create acquisition ID image raster (-id.tif)
Expand Down Expand Up @@ -748,10 +748,10 @@ def calc_product_start_stop(src_ids, extent, epsg):


def create_data_mask(outname, datasets, extent, epsg, driver, creation_opt,
overviews, overview_resampling, dst_nodata, wbm=None):
overviews, overview_resampling, dst_nodata, wbm=None, orb=False):
"""
maawoo marked this conversation as resolved.
Show resolved Hide resolved
Creation of the Data Mask image.

Parameters
----------
outname: str
Expand All @@ -774,7 +774,7 @@ def create_data_mask(outname, datasets, extent, epsg, driver, creation_opt,
Resampling method for overview levels.
dst_nodata: int or str
Nodata value to write to the output raster.
wbm: str or None
wbm: str or None, optional
Path to a water body mask file with the dimensions of an MGRS tile.
"""
measurement_keys = [x for x in datasets[0].keys() if re.search('[gs]-lin', x)]
Expand All @@ -787,14 +787,13 @@ def create_data_mask(outname, datasets, extent, epsg, driver, creation_opt,
else:
return # do not create a data mask if not all scenes have a layover-shadow mask
print(outname)
dm_bands = {1: {'arr_val': 0,
'name': 'not layover, nor shadow'},
2: {'arr_val': 1,
'name': 'layover'},
3: {'arr_val': 2,
'name': 'shadow'},
4: {'arr_val': 4,
'name': 'ocean water'}}
dm_bands = [{'arr_val': 0, 'name': 'not layover, nor shadow'},
{'arr_val': 1, 'name': 'layover'},
{'arr_val': 2, 'name': 'shadow'},
# {'arr_val': 3, 'name': 'layover and shadow'}, # just for context, not used as an individual band
{'arr_val': 4, 'name': 'ocean water'}] if not orb else \
[{'arr_val': [0, 1, 2], 'name': 'land'},
{'arr_val': 4, 'name': 'ocean water'}]

tile_bounds = [extent['xmin'], extent['ymin'], extent['xmax'], extent['ymax']]

Expand Down Expand Up @@ -831,9 +830,9 @@ def create_data_mask(outname, datasets, extent, epsg, driver, creation_opt,
arr_wbm = ras_wbm.array()
out_arr = np.where((arr_wbm == 1), 4, arr_dm)
del arr_wbm
else:
elif wbm is None and not orb:
out_arr = arr_dm
dm_bands.pop(4)
del dm_bands[3]
del arr_dm

# Extend the shadow class of the data mask with nodata values from backscatter data and create final array
Expand All @@ -851,25 +850,25 @@ def create_data_mask(outname, datasets, extent, epsg, driver, creation_opt,

outname_tmp = '/vsimem/' + os.path.basename(outname) + '.vrt'
gdriver = gdal.GetDriverByName('GTiff')
ds_tmp = gdriver.Create(outname_tmp, rows, cols, len(dm_bands.keys()), gdal.GDT_Byte,
ds_tmp = gdriver.Create(outname_tmp, rows, cols, len(dm_bands), gdal.GDT_Byte,
options=['ALPHA=UNSPECIFIED', 'PHOTOMETRIC=MINISWHITE'])
gdriver = None
ds_tmp.SetGeoTransform(geotrans)
ds_tmp.SetProjection(proj)

for k, v in dm_bands.items():
band = ds_tmp.GetRasterBand(k)
arr_val = v['arr_val']
b_name = v['name']
for i, _dict in enumerate(dm_bands):
band = ds_tmp.GetRasterBand(i+1)
arr_val = _dict['arr_val']
b_name = _dict['name']

arr = np.full((rows, cols), 0)
arr[out_arr == dst_nodata] = dst_nodata
if arr_val == 0:
arr[out_arr == 0] = 1
if arr_val in [0, 4]:
arr[out_arr == arr_val] = 1
elif arr_val in [1, 2]:
arr[(out_arr == arr_val) | (out_arr == 3)] = 1
elif arr_val == 4:
arr[out_arr == 4] = 1
elif arr_val == [0, 1, 2]:
arr[out_arr != 4] = 1

arr = arr.astype('uint8')
band.WriteArray(arr)
Expand Down