Skip to content

Commit

Permalink
DEM and layerdatamask
Browse files Browse the repository at this point in the history
  • Loading branch information
@gr4n0t4
gr4n0t4 committed Jul 28, 2022
1 parent 719be24 commit 2450b77
Show file tree
Hide file tree
Showing 6 changed files with 227 additions and 138 deletions.
21 changes: 16 additions & 5 deletions ERS_NRB/ancillary.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@
from scipy.interpolate import griddata
import spatialist
from spatialist import gdalbuildvrt, Raster, bbox
from rsgislib.imageutils import create_mosaic_images_vrt
import pyroSAR
from pyroSAR import identify, finder, examine
from ERS_NRB.metadata.extract import get_uid_sid, etree_from_sid, find_in_annotation
Expand Down Expand Up @@ -376,13 +377,15 @@ def create_data_mask(outname, valid_mask_list, src_files, extent, epsg, driver,

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

vrt_snap_ls = '/vsimem/' + os.path.dirname(outname) + 'snap_ls.vrt'
vrt_snap_valid = '/vsimem/' + os.path.dirname(outname) + 'snap_valid.vrt'
vrt_snap_gamma0 = '/vsimem/' + os.path.dirname(outname) + 'snap_gamma0.vrt'
vrt_snap_ls = os.path.dirname(outname) + '/snap_ls.vrt'
vrt_snap_valid = os.path.dirname(outname) + '/snap_valid.vrt'
vrt_snap_gamma0 = os.path.dirname(outname) + '/snap_gamma0.vrt'
gdalbuildvrt(snap_ls_mask, vrt_snap_ls, options={'outputBounds': tile_bounds}, void=False)
gdalbuildvrt(valid_mask_list, vrt_snap_valid, options={'outputBounds': tile_bounds}, void=False)
gdalbuildvrt(snap_gamma0, vrt_snap_gamma0, options={'outputBounds': tile_bounds}, void=False)

print(f"vrt_snap_ls: {vrt_snap_ls}")
print(f"snap_ls_mask: {snap_ls_mask}")

with Raster(vrt_snap_ls) as ras_snap_ls:
with bbox(extent, crs=epsg) as tile_vec:
rows = ras_snap_ls.rows
Expand All @@ -394,6 +397,9 @@ def create_data_mask(outname, valid_mask_list, src_files, extent, epsg, driver,
# Add Water Body Mask
if wbm is not None:
with Raster(wbm) as ras_wbm:
print(f"ras_wbm: {ras_wbm}")
print(f"ras_snap_ls: {ras_snap_ls}")

arr_wbm = ras_wbm.array()
out_arr = np.where((arr_wbm == 1), 4, arr_snap_dm)
del arr_wbm
Expand All @@ -418,14 +424,16 @@ def create_data_mask(outname, valid_mask_list, src_files, extent, epsg, driver,
if out_format == 'multi-layer':
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,
print(f"len(dm_bands.keys()): {len(dm_bands.keys())}")
ds_tmp = gdriver.Create(outname_tmp, cols, rows, len(dm_bands.keys()), 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)
print(f"band: {band.XSize}x{band.YSize}")
arr_val = v['arr_val']
b_name = v['name']

Expand All @@ -439,6 +447,7 @@ def create_data_mask(outname, valid_mask_list, src_files, extent, epsg, driver,
arr[out_arr == 4] = 1

arr = arr.astype('uint8')
print(f"arr: {arr.shape}x{arr.dtype}")
band.WriteArray(arr)
band.SetNoDataValue(out_nodata)
band.SetDescription(b_name)
Expand Down Expand Up @@ -529,6 +538,8 @@ def create_acq_id_image(ref_tif, valid_mask_list, src_scenes, extent, epsg, driv
creation_opt.append('TIFFTAG_IMAGEDESCRIPTION={}'.format(tag))

with Raster(ref_tif) as ref_ras:
print(f"ref_ras: {ref_ras}")
print(f"out_arr: {out_arr.shape}")
ref_ras.write(outname, format=driver, array=out_arr.astype('uint8'), nodata=out_nodata, overwrite=True,
overviews=overviews, options=creation_opt)

Expand Down
4 changes: 2 additions & 2 deletions ERS_NRB/config.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,7 +71,7 @@ def get_config(config_file, section_name='GENERAL'):
v = int(v)
if k == 'dem_type':
allowed = ['Copernicus 10m EEA DEM', 'Copernicus 30m Global DEM II',
'Copernicus 30m Global DEM', 'GETASSE30']
'Copernicus 30m Global DEM', 'GETASSE30', "Copernicus 90m Global DEM II"]
assert v in allowed, "Parameter '{}': expected to be one of {}; got '{}' instead".format(k, allowed, v)
out_dict[k] = v

Expand Down Expand Up @@ -145,7 +145,7 @@ def geocode_conf(config):
'clean_edges': True,
'clean_edges_npixels': 3,
'test': False,
'cleanup': True,
'cleanup': False,
'rlks': {'IMM': 5, # TODO Completly guess
'IMP': 6, # TODO Completly guess
'APP': 3}[config['acq_mode']], # TODO Completly guess
Expand Down
16 changes: 12 additions & 4 deletions ERS_NRB/metadata/extract.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,11 +40,14 @@ def get_prod_meta(product_id, tif, src_scenes, src_dir):
if tif:
with vec_from_srccoords(coord_list=coord_list) as srcvec:
with Raster(tif) as ras:
image = ras.raster.GetGeoTransform()
vec = ras.bbox()
srs = vec.srs
out['extent'] = vec.extent
out['wkt'] = srs.ExportToWkt()
out['epsg'] = vec.getProjection(type='epsg')
out['rowSpacing'] = image[5]
out['columnSpacing'] = image[1]
out['rows'] = ras.rows
out['cols'] = ras.cols
out['res'] = ras.res
Expand Down Expand Up @@ -479,7 +482,6 @@ def meta_dict(config, target, src_scenes, src_files, proc_time):

# Product metadata (sorted alphabetically)
meta['prod']['access'] = None
meta['prod']['ancillaryData1'] = None
meta['prod']['acquisitionType'] = 'NOMINAL'
# meta['prod']['ascendingNodeDate'] = manifest0.find('.//s1:ascendingNodeTime', nsmap0).text
# meta['prod']['azimuthNumberOfLooks'] = prod_meta['ML_nAzLooks']
Expand All @@ -500,6 +502,7 @@ def meta_dict(config, target, src_scenes, src_files, proc_time):
meta['prod']['demType'] = dem_type
meta['prod']['demURL'] = dem_url
meta['prod']['doi'] = None
meta['prod']['ancillaryData1'] = meta['prod']['demReference']
meta['prod']['ellipsoidalHeight'] = None
meta['prod']['fileBitsPerSample'] = '32'
meta['prod']['fileByteOrder'] = 'little-endian'
Expand Down Expand Up @@ -527,18 +530,21 @@ def meta_dict(config, target, src_scenes, src_files, proc_time):
meta['prod']['griddingConvention'] = 'Military Grid Reference System (MGRS)'
meta['prod']['licence'] = None
meta['prod']['majorCycleID'] = str(sid0.meta['cycleNumber'])
meta['prod']['noiseRemovalApplied'] = True
meta['prod']['noiseRemovalApplied'] = False
meta['prod']['noiseRemovalAlgorithm'] = 'https://doi.org/10.1109/tgrs.2018.2889381'
meta['prod']['numberLines'] = str(prod_meta['rows'])
meta['prod']['numberOfAcquisitions'] = str(len(src_scenes))
meta['prod']['numBorderPixels'] = prod_meta['nodata_borderpx']
meta['prod']['numPixelsPerLine'] = str(prod_meta['cols'])
meta['prod']['columnSpacing'] = str(prod_meta['columnSpacing'])
meta['prod']['rowSpacing'] = str(prod_meta['rowSpacing'])
meta['prod']['pixelCoordinateConvention'] = 'pixel ULC'
meta['prod']['processingCenter'] = 'FSU'
meta['prod']['processingCenter'] = 'TBD'
meta['prod']['location'] = 'TBD'
meta['prod']['processingLevel'] = 'Level 2'
meta['prod']['processingMode'] = 'PROTOTYPE'
meta['prod']['processorName'] = 'ERS_NRB'
meta['prod']['processorVersion'] = '0.1'
meta['prod']['processorVersion'] = 'TBD'
meta['prod']['productName'] = 'NORMALISED RADAR BACKSCATTER'
meta['prod']['pxSpacingColumn'] = str(prod_meta['res'][0])
meta['prod']['pxSpacingRow'] = str(prod_meta['res'][1])
Expand Down Expand Up @@ -678,6 +684,8 @@ def meta_dict(config, target, src_scenes, src_files, proc_time):
# meta['source'][uid]['swaths'] = swaths
meta['source'][uid]['incidenceAngleMax'] = src_sid[uid].meta['incidenceAngleMax']
meta['source'][uid]['incidenceAngleMin'] = src_sid[uid].meta['incidenceAngleMin']
meta['source'][uid]['neszNear'] = src_sid[uid].meta['neszNear']
meta['source'][uid]['neszFar'] = src_sid[uid].meta['neszFar']
meta['source'][uid]['incidenceAngleMidSwath'] = src_sid[uid].meta['incidence']
# return meta, m_sid, m_src
return meta
30 changes: 21 additions & 9 deletions ERS_NRB/metadata/stacparser.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,9 @@ def product_json(meta, target, tifs):
bbox=meta['prod']['geom_stac_bbox_native'],
shape=[int(meta['prod']['numPixelsPerLine']), int(meta['prod']['numberLines'])],
transform=meta['prod']['transform'])


item.properties['processing:location'] = meta['prod']['location']

item.properties['processing:facility'] = meta['prod']['processingCenter']
item.properties['processing:software'] = {meta['prod']['processorName']: meta['prod']['processorVersion']}
item.properties['processing:level'] = meta['prod']['processingLevel']
Expand All @@ -95,7 +97,7 @@ def product_json(meta, target, tifs):
'window_size_col': meta['prod']['filterWindowSizeCol'],
'window_size_line': meta['prod']['filterWindowSizeLine']}
else:
item.properties['card4l:speckle_filtering'] = None
item.properties['card4l:speckle_filtering'] = "No"
item.properties['card4l:noise_removal_applied'] = meta['prod']['noiseRemovalApplied']
item.properties['card4l:conversion_eq'] = meta['prod']['backscatterConversionEq']
item.properties['card4l:relative_radiometric_accuracy'] = meta['prod']['radiometricAccuracyRelative']
Expand All @@ -104,6 +106,8 @@ def product_json(meta, target, tifs):
item.properties['card4l:dem_resampling_method'] = meta['prod']['demResamplingMethod']
item.properties['card4l:egm_resampling_method'] = meta['prod']['demEgmResamplingMethod']
item.properties['card4l:geometric_accuracy_type'] = meta['prod']['geoCorrAccuracyType']
item.properties['card4l:column_spacing'] = meta['prod']['columnSpacing']
item.properties['card4l:row_spacing'] = meta['prod']['rowSpacing']
for x in ['Northern', 'Eastern']:
key = ['geoCorrAccuracy{}{}'.format(x, y) for y in ['STDev', 'Bias']]
stddev = float(meta['prod'][key[0]]) if meta['prod'][key[0]] is not None else None
Expand Down Expand Up @@ -150,6 +154,9 @@ def product_json(meta, target, tifs):
item.add_link(link=pystac.Link(rel='radiometric-terrain-correction',
target=meta['prod']['RTCAlgorithm'],
title='Reference to the Radiometric Terrain Correction algorithm details.'))
item.add_link(link=pystac.Link(rel='rtc-dem',
target=meta['prod']['demReference'],
title='The Digital Elevation Model used during Radiometric Terrain Correction'))
item.add_link(link=pystac.Link(rel='radiometric-accuracy',
target=meta['prod']['radiometricAccuracyReference'],
title='Reference describing the radiometric uncertainty of the product.'))
Expand Down Expand Up @@ -183,7 +190,10 @@ def product_json(meta, target, tifs):
pol = re.search('[vh]{2}', tif).group().lower()
created = datetime.fromtimestamp(os.path.getctime(tif)).isoformat()
extra_fields = {'created': created,
'raster:bands': [{'nodata': 'NaN',
'measurement_type': meta['prod']['backscatterMeasurement'],
'backscatter_convention': meta['prod']['backscatterConvention'],
'raster:bands': [{'unit':'natural',
'nodata': 'NaN',
'data_type': '{}{}'.format(meta['prod']['fileDataType'],
meta['prod']['fileBitsPerSample']),
'bits_per_sample': int(meta['prod']['fileBitsPerSample'])}],
Expand Down Expand Up @@ -365,6 +375,8 @@ def source_json(meta, target):
item.properties['card4l:source_geometry'] = meta['source'][uid]['dataGeometry']
item.properties['card4l:incidence_angle_near_range'] = meta['source'][uid]['incidenceAngleMin']
item.properties['card4l:incidence_angle_far_range'] = meta['source'][uid]['incidenceAngleMax']
item.properties['card4l:nesz_near'] = meta['source'][uid]['neszNear']
item.properties['card4l:nesz_far'] = meta['source'][uid]['neszFar']
item.properties['card4l:noise_equivalent_intensity'] = meta['source'][uid]['perfEstimates']
item.properties['card4l:noise_equivalent_intensity_type'] = meta['source'][uid]['perfNoiseEquivalentIntensityType']
item.properties['card4l:mean_faraday_rotation_angle'] = meta['source'][uid]['faradayMeanRotationAngle']
Expand Down Expand Up @@ -406,12 +418,12 @@ def source_json(meta, target):
title='Reference describing the method used to derive the estimate for the mean'
' Faraday rotation angle.'))

xml_relpath = './' + os.path.relpath(outname.replace('.json', '.xml'), target).replace('\\', '/')
item.add_asset(key='card4l',
asset=pystac.Asset(href=xml_relpath,
title='CARD4L XML Metadata File',
media_type=pystac.MediaType.XML,
roles=['metadata', 'card4l']))
# xml_relpath = './' + os.path.relpath(outname.replace('.json', '.xml'), target).replace('\\', '/')
# item.add_asset(key='card4l',
# asset=pystac.Asset(href=xml_relpath,
# title='CARD4L XML Metadata File',
# media_type=pystac.MediaType.XML,
# roles=['metadata', 'card4l']))
print(outname)
item.save_object(dest_href=outname)

Expand Down
Loading

0 comments on commit 2450b77

Please sign in to comment.