diff --git a/process/1_process_add.py b/process/1_process_add.py
index 6fe33e8c09bf9fb2d4a92b915de90ab3ce1462a5..722c17dc797928faf1ff99f2a10de5cab755c07f 100644
--- a/process/1_process_add.py
+++ b/process/1_process_add.py
@@ -22,8 +22,8 @@ and all geometries are fixed (with zero buffer).
#=============================================================================
__author__ = 'Simon Dreutter'
-__version__ = '0.1'
-__date__ = '2020-12-17'
+__version__ = '0.3'
+__date__ = '2021-11-05'
__email__ = 'simon.dreutter@awi.de'
__status__ = 'Developement'
@@ -55,24 +55,55 @@ add_coastline_in_name = os.path.splitext(os.path.basename(ADD_COASTLINE_IN))[0]
add_ice_sheet_name = os.path.splitext(os.path.basename(ADD_ICE_SHEET))[0]
add_shelf_ice_name = os.path.splitext(os.path.basename(ADD_SHELF_ICE))[0]
add_rock_outcrop_in_name = os.path.splitext(os.path.basename(ADD_ROCK_OUTCROP_IN))[0]
+# DEBUG
# inconsistancy with add 7.3 layer names, variable needs to be overwritten
+# remove the DEBUG block once this is fixed in the ADD layer
add_rock_outcrop_in_name = 'add_rock_outcrop_medium_res_polygon_v7'
+#/DEBUG
add_rock_outcrop_name = os.path.splitext(os.path.basename(ADD_ROCK_OUTCROP))[0]
# heading
print()
-print('============================================================')
-print(' ANTARCTIC DIGITAL DATANASE / ADD ')
+print('====================================================================')
+print(f'{"ANTARCTIC DIGITAL DATANASE / ADD":^68}')
msg = 'Reprojecting and extracting land mask from ADD coastlines...'
-cmd = f'ogr2ogr -s_srs {EPSG_ADD} -t_srs {EPSG_WORLD} -wrapdateline -nln "{add_ice_sheet_name}" -nlt MULTIPOLYGON -f GPKG -overwrite {ADD_ICE_SHEET} {ADD_COASTLINE_IN} -dialect SQLite -sql "SELECT ST_Buffer(ST_Multi(geom),0.0) from \'{add_coastline_in_name}\' WHERE surface=\'land\'"'
+cmd = (f'ogr2ogr'
+ f' -s_srs {EPSG_ADD}'
+ f' -t_srs {EPSG_WORLD}'
+ f' -wrapdateline -nln "{add_ice_sheet_name}"'
+ f' -nlt MULTIPOLYGON'
+ f' -f GPKG'
+ f' -overwrite'
+ f' {ADD_ICE_SHEET} {ADD_COASTLINE_IN}'
+ f' -dialect SQLite'
+ f' -sql "SELECT ST_Buffer(ST_Multi(geom),0.0) from \'{add_coastline_in_name}\' WHERE surface=\'land\'"')
run(msg,cmd)
msg = 'Reprojecting and extracting ice mask from ADD coastlines...'
-cmd = f'ogr2ogr -s_srs {EPSG_ADD} -t_srs {EPSG_WORLD} -wrapdateline -nln "{add_shelf_ice_name}" -nlt MULTIPOLYGON -f GPKG -overwrite {ADD_SHELF_ICE} {ADD_COASTLINE_IN} -dialect SQLite -sql "SELECT ST_Buffer(ST_Multi(geom),0.0) from \'{add_coastline_in_name}\' WHERE surface IN (\'ice shelf\',\'rumple\',\'ocean\',\'ice tongue\')"'
+cmd = (f'ogr2ogr'
+ f' -s_srs {EPSG_ADD}'
+ f' -t_srs {EPSG_WORLD}'
+ f' -wrapdateline'
+ f' -nln "{add_shelf_ice_name}"'
+ f' -nlt MULTIPOLYGON'
+ f' -f GPKG -overwrite'
+ f' {ADD_SHELF_ICE} {ADD_COASTLINE_IN}'
+ f' -dialect SQLite'
+ f' -sql "SELECT ST_Buffer(ST_Multi(geom),0.0) from \'{add_coastline_in_name}\' WHERE surface IN (\'ice shelf\',\'rumple\',\'ocean\',\'ice tongue\')"')
run(msg,cmd)
msg = 'Reprojecting ADD rock outcrops...'
-cmd = f'ogr2ogr -s_srs {EPSG_ADD} -t_srs {EPSG_WORLD} -wrapdateline -nln "{add_rock_outcrop_name}" -nlt MULTIPOLYGON -f GPKG -overwrite {ADD_ROCK_OUTCROP} {ADD_ROCK_OUTCROP_IN} -dialect SQLite -sql "SELECT ST_Buffer(ST_Multi(geom),0.0) from \'{add_rock_outcrop_in_name}\'"'
+cmd = (f'ogr2ogr'
+ f' -s_srs {EPSG_ADD}'
+ f' -t_srs {EPSG_WORLD}'
+ f' -wrapdateline'
+ f' -nln "{add_rock_outcrop_name}"'
+ f' -nlt MULTIPOLYGON'
+ f' -f GPKG'
+ f' -overwrite'
+ f' {ADD_ROCK_OUTCROP} {ADD_ROCK_OUTCROP_IN}'
+ f' -dialect SQLite'
+ f' -sql "SELECT ST_Buffer(ST_Multi(geom),0.0) from \'{add_rock_outcrop_in_name}\'"')
run(msg,cmd)
diff --git a/process/2_process_gimp.py b/process/2_process_gimp.py
index 5e8a2f2182671a1ce55cc960c84920a993fe909f..d258d64236b902d714b87413ea32da1f6ccc4816 100644
--- a/process/2_process_gimp.py
+++ b/process/2_process_gimp.py
@@ -20,8 +20,8 @@ handling and then polygonized to get a vector dataset.
#=============================================================================
__author__ = 'Simon Dreutter'
-__version__ = '0.1'
-__date__ = '2020-12-17'
+__version__ = '0.3'
+__date__ = '2021-11-05'
__email__ = 'simon.dreutter@awi.de'
__status__ = 'Developement'
@@ -45,15 +45,26 @@ gimp_name = os.path.splitext(os.path.basename(GIMP))[0]
# heading
print()
-print('============================================================')
-print(' GREENLAND ICE MAPPING PROJECT / GIMP ')
+print('====================================================================')
+print(f'{"GREENLAND ICE MAPPING PROJECT / GIMP":^68}')
msg = 'Reprojecting and downsampling GIMP raster...'
-cmd = f'gdalwarp -t_srs {EPSG_WORLD} -dstnodata 0.0 -tr {GIMP_RESOLUTION} {GIMP_RESOLUTION} -r near -of GTiff {GIMP_IN} {GIMP_RASTER}'
+cmd = (f'gdalwarp'
+ f' -t_srs {EPSG_WORLD}'
+ f' -dstnodata 0.0'
+ f' -tr {GIMP_RESOLUTION} {GIMP_RESOLUTION}'
+ f' -r near'
+ f' -of GTiff'
+ f' {GIMP_IN} {GIMP_RASTER}')
run(msg,cmd)
msg = 'Polygonizing GIMP raster...'
-cmd = f'{PYTHON} -m gdal_polygonize {GIMP_RASTER} {GIMP} -b 1 -f "GPKG" {gimp_name} DN'
+cmd = (f'gdal_polygonize'
+ f' {GIMP_RASTER} {GIMP}'
+ f' -b 1'
+ f' -f "GPKG"'
+ f' {gimp_name} DN')
run(msg,cmd)
+# cleanup
cleanup(GIMP_RASTER)
diff --git a/process/3_process_glims.py b/process/3_process_glims.py
index 191207b6f618608ecb2d3aa38c690e2e2cb41251..b36249b8fef7ff469cec8d16431e9921854ac4cd 100644
--- a/process/3_process_glims.py
+++ b/process/3_process_glims.py
@@ -21,8 +21,8 @@ Arctic and the Antarctic.
#=============================================================================
__author__ = 'Simon Dreutter'
-__version__ = '0.1'
-__date__ = '2020-12-17'
+__version__ = '0.3'
+__date__ = '2021-11-05'
__email__ = 'simon.dreutter@awi.de'
__status__ = 'Developement'
@@ -49,17 +49,30 @@ glims_antarctic_name = os.path.splitext(os.path.basename(GLIMS_ANTARCTIC))[0]
# heading
print()
-print('============================================================')
-print(' GLOBAL LAND ICE MEASUREMENTS FROM SPACE / GLIMS ')
+print('====================================================================')
+print(f'{"GLOBAL LAND ICE MEASUREMENTS FROM SPACE / GLIMS":^68}')
msg = 'Fixing GLIMS geometries...'
-cmd = f'ogr2ogr {GLIMS_WORLD} {GLIMS_IN} -dialect sqlite -sql "SELECT ST_Buffer(geometry, 0.0) AS geometry,* FROM \"{glims_in_name}\"" -nln "{glims_world_name}" -f "GPKG"'
+cmd = (f'ogr2ogr'
+ f' {GLIMS_WORLD} {GLIMS_IN}'
+ f' -dialect sqlite'
+ f' -sql "SELECT ST_Union(ST_Buffer(geometry, 0.0)) AS geometry,* FROM \"{glims_in_name}\"" -nln "{glims_world_name}" -f "GPKG"')
run(msg,cmd)
msg = 'Clipping GLIMS for the Arctic...'
-cmd = f'ogr2ogr -spat -180.0 90.0 180.0 {ARCTIC_EXTENT_LAT} -clipsrc spat_extent {GLIMS_ARCTIC} {GLIMS_WORLD} {glims_world_name} -nln "{glims_arctic_name}" -f "GPKG"'
+cmd = (f'ogr2ogr'
+ f' -spat -180.0 90.0 180.0 {ARCTIC_EXTENT_LAT}'
+ f' -clipsrc spat_extent'
+ f' {GLIMS_ARCTIC} {GLIMS_WORLD} {glims_world_name}'
+ f' -nln "{glims_arctic_name}"'
+ f' -f "GPKG"')
run(msg,cmd)
msg = 'Clipping GLIMS for the Antarctic...'
-cmd = f'ogr2ogr -spat -180.0 {ANTARCTIC_EXTENT_LAT} 180.0 -90.0 -clipsrc spat_extent {GLIMS_ANTARCTIC} {GLIMS_WORLD} {glims_world_name} -nln "{glims_antarctic_name}" -f "GPKG"'
+cmd = (f'ogr2ogr'
+ f' -spat -180.0 {ANTARCTIC_EXTENT_LAT} 180.0 -90.0'
+ f' -clipsrc spat_extent'
+ f' {GLIMS_ANTARCTIC} {GLIMS_WORLD} {glims_world_name}'
+ f' -nln "{glims_antarctic_name}"'
+ f' -f "GPKG"')
run(msg,cmd)
diff --git a/process/4_process_gebco.py b/process/4_process_gebco.py
index 57d90b45756580a55bc14e5be7fe2936e16fb09b..167b6dcffa539758b54993e63f8f95afba7e92b5 100644
--- a/process/4_process_gebco.py
+++ b/process/4_process_gebco.py
@@ -24,8 +24,8 @@ target resolution (RESOLUTION_M).
#=============================================================================
__author__ = 'Simon Dreutter'
-__version__ = '0.1'
-__date__ = '2020-12-17'
+__version__ = '0.3'
+__date__ = '2021-11-05'
__email__ = 'simon.dreutter@awi.de'
__status__ = 'Developement'
@@ -34,6 +34,8 @@ __status__ = 'Developement'
# IMPORT
#=============================================================================
+import shutil
+
# import AWI Basemap configuration
from config import *
@@ -44,38 +46,198 @@ from config import *
# heading
print()
-print('============================================================')
-print(' GENERAL BATHYMETRIC CHART OF THE OCEANS / GEBCO ')
+print('====================================================================')
+print(f'{"GENERAL BATHYMETRIC CHART OF THE OCEANS / GEBCO":^68}')
+
+# World
+# ice surface
+msg = 'Resampling GEBCO ice surface...'
+cmd = (f'gdalwarp'
+ f' -s_srs {EPSG_WORLD}'
+ f' -t_srs {EPSG_WORLD}'
+ f' -tr {RESOLUTION_DEG} {RESOLUTION_DEG}'
+ f' -r cubic'
+ f' -of GTiff'
+ f' {WARP_OPTIONS}'
+ f' {GDAL_CREATE_OPTIONS}'
+ f' {GEBCO_IN} {GEBCO_WORLD_ICESURFACE_TMP1}')
+run(msg,cmd)
+
+msg = 'Copying World DEM...'
+print('====================================================================')
+print()
+print(msg)
+print()
+shutil.copy(GEBCO_WORLD_ICESURFACE_TMP1, GEBCO_WORLD)
+
+msg = 'Extracting World DEM above sealevel...'
+cmd = (f'gdal_calc'
+ f' --calc "(A >= 0) * A"'
+ f' --format GTiff'
+ f' -A {GEBCO_WORLD_ICESURFACE_TMP1}'
+ f' --A_band 1'
+ f' --overwrite'
+ f' --outfile {GEBCO_WORLD_ICESURFACE_TMP2}'
+ f' --NoDataValue=0')
+run(msg,cmd)
+
+msg = 'Fixing Antarctic DEM above sealevel NoData...'
+cmd = (f'gdalwarp'
+ f' -overwrite'
+ f' -srcnodata 0'
+ f' {GEBCO_WORLD_ICESURFACE_TMP2}'
+ f' {GEBCO_WORLD_ICESURFACE}')
+run(msg,cmd)
-msg = 'Resampling GEBCO...'
-cmd = f'gdalwarp -s_srs {EPSG_WORLD} -t_srs {EPSG_WORLD} -tr {RESOLUTION_DEG} {RESOLUTION_DEG} -r cubic -of GTiff {WARP_OPTIONS} {GDAL_CREATE_OPTIONS} {GEBCO_IN} {GEBCO_WORLD}'
+# sub ice
+msg = 'Resampling GEBCO sub ice...'
+cmd = (f'gdalwarp'
+ f' -s_srs {EPSG_WORLD}'
+ f' -t_srs {EPSG_WORLD}'
+ f' -tr {RESOLUTION_DEG} {RESOLUTION_DEG}'
+ f' -r cubic'
+ f' -of GTiff'
+ f' {WARP_OPTIONS}'
+ f' {GDAL_CREATE_OPTIONS}'
+ f' {GEBCO_SUBICE_IN} {GEBCO_WORLD_SUBICE}')
run(msg,cmd)
+# Arctic
msg = 'Clipping GEBCO for the Arctic...'
-cmd = f'gdal_translate -projwin -180.0 90.0 180.0 {ARCTIC_EXTENT_LAT} -a_srs {EPSG_WORLD} -of VRT {GEBCO_IN} {GEBCO_ARCTIC_CLIP_VRT}'
+cmd = (f'gdal_translate'
+ f' -projwin -180.0 90.0 180.0 {ARCTIC_EXTENT_LAT}'
+ f' -a_srs {EPSG_WORLD}'
+ f' -of VRT'
+ f' {GEBCO_IN} {GEBCO_ARCTIC_CLIP_VRT}')
run(msg,cmd)
msg = 'Reprojecting GEBCO for the Arctic...'
-cmd = f'gdalwarp -s_srs {EPSG_WORLD} -t_srs {EPSG_ARCTIC} -r near -of VRT {WARP_OPTIONS} {GEBCO_ARCTIC_CLIP_VRT} {GEBCO_ARCTIC_FULL_VRT}'
+cmd = (f'gdalwarp'
+ f' -s_srs {EPSG_WORLD}'
+ f' -t_srs {EPSG_ARCTIC}'
+ f' -r near'
+ f' -of VRT'
+ f' {WARP_OPTIONS}'
+ f' {GEBCO_ARCTIC_CLIP_VRT} {GEBCO_ARCTIC_FULL_VRT}')
run(msg,cmd)
msg = 'Resampling Arctic to target resolution...'
-cmd = f'gdalwarp -s_srs {EPSG_ARCTIC} -t_srs {EPSG_ARCTIC} -tr {RESOLUTION_M} {RESOLUTION_M} -r cubic -of GTiff {WARP_OPTIONS} {GDAL_CREATE_OPTIONS} {GEBCO_ARCTIC_FULL_VRT} {GEBCO_ARCTIC}'
+cmd = (f'gdalwarp'
+ f' -s_srs {EPSG_ARCTIC}'
+ f' -t_srs {EPSG_ARCTIC}'
+ f' -tr {RESOLUTION_M} {RESOLUTION_M}'
+ f' -r cubic'
+ f' -of GTiff'
+ f' {WARP_OPTIONS}'
+ f' {GDAL_CREATE_OPTIONS}'
+ f' {GEBCO_ARCTIC_FULL_VRT} {GEBCO_ARCTIC}')
+run(msg,cmd)
+
+# Antarctic
+# Due to ice shelfs there is a very steep drop in the DEM from shelf ice edge
+# down to the seabed, which creates hillshade artifacts. To avoid these, the
+# Antarctic GEBCO grid will processed twice, once the ice surface DEM and once
+# the sub ice topography DEM. The main hillshade is created from the sub ice
+# DEM. Then all data above sea level is extracted from the ice surface DEM, a
+# hillshade is created for these areas and the sub ice hillshade is replaced
+# in these areas.
+
+# ice surface
+msg = 'Clipping GEBCO ice surface for the Antarctic...'
+cmd = (f'gdal_translate'
+ f' -projwin -180.0 {ANTARCTIC_EXTENT_LAT} 180.0 -90.0'
+ f' -a_srs {EPSG_WORLD}'
+ f' -of VRT'
+ f' {GEBCO_IN} {GEBCO_ANTARCTIC_ICESURFACE_CLIP_VRT}')
+run(msg,cmd)
+
+msg = 'Reprojecting GEBCO ice surface for the Antarctic...'
+cmd = (f'gdalwarp'
+ f' -s_srs {EPSG_WORLD}'
+ f' -t_srs {EPSG_ANTARCTIC}'
+ f' -r near'
+ f' -of VRT'
+ f' {WARP_OPTIONS}'
+ f' {GEBCO_ANTARCTIC_ICESURFACE_CLIP_VRT} {GEBCO_ANTARCTIC_ICESURFACE_FULL_VRT}')
+run(msg,cmd)
+
+msg = 'Resampling Antarctic ice surface to target resolution...'
+cmd = (f'gdalwarp'
+ f' -s_srs {EPSG_ANTARCTIC}'
+ f' -t_srs {EPSG_ANTARCTIC}'
+ f' -tr {RESOLUTION_M} {RESOLUTION_M}'
+ f' -r cubic'
+ f' -of GTiff'
+ f' {WARP_OPTIONS}'
+ f' {GDAL_CREATE_OPTIONS}'
+ f' {GEBCO_ANTARCTIC_ICESURFACE_FULL_VRT} {GEBCO_ANTARCTIC_ICESURFACE_TMP1}')
+run(msg,cmd)
+
+msg = 'Copying Antarctic DEM...'
+print('====================================================================')
+print()
+print(msg)
+print()
+shutil.copy(GEBCO_ANTARCTIC_ICESURFACE_TMP1, GEBCO_ANTARCTIC)
+
+msg = 'Extracting Antarctic DEM above sealevel...'
+cmd = (f'gdal_calc'
+ f' --calc "(A >= 0) * A"'
+ f' --format GTiff'
+ f' -A {GEBCO_ANTARCTIC_ICESURFACE_TMP1}'
+ f' --A_band 1'
+ f' --overwrite'
+ f' --outfile {GEBCO_ANTARCTIC_ICESURFACE_TMP2}'
+ f' --NoDataValue=0')
+run(msg,cmd)
+
+msg = 'Fixing Antarctic DEM above sealevel NoData...'
+cmd = (f'gdalwarp'
+ f' -overwrite'
+ f' -srcnodata 0'
+ f' {GEBCO_ANTARCTIC_ICESURFACE_TMP2}'
+ f' {GEBCO_ANTARCTIC_ICESURFACE}')
run(msg,cmd)
-msg = 'Clipping GEBCO for the Antarctic...'
-cmd = f'gdal_translate -projwin -180.0 {ANTARCTIC_EXTENT_LAT} 180.0 -90.0 -a_srs {EPSG_WORLD} -of VRT {GEBCO_IN} {GEBCO_ANTARCTIC_CLIP_VRT}'
+# sub ice
+msg = 'Clipping GEBCO sub ice for the Antarctic...'
+cmd = (f'gdal_translate'
+ f' -projwin -180.0 {ANTARCTIC_EXTENT_LAT} 180.0 -90.0'
+ f' -a_srs {EPSG_WORLD}'
+ f' -of VRT'
+ f' {GEBCO_SUBICE_IN} {GEBCO_ANTARCTIC_SUBICE_CLIP_VRT}')
run(msg,cmd)
-msg = 'Reprojecting GEBCO for the Antarctic...'
-cmd = f'gdalwarp -s_srs {EPSG_WORLD} -t_srs {EPSG_ANTARCTIC} -r near -of VRT {WARP_OPTIONS} {GEBCO_ANTARCTIC_CLIP_VRT} {GEBCO_ANTARCTIC_FULL_VRT}'
+msg = 'Reprojecting GEBCO sub ice for the Antarctic...'
+cmd = (f'gdalwarp'
+ f' -s_srs {EPSG_WORLD}'
+ f' -t_srs {EPSG_ANTARCTIC}'
+ f' -r near'
+ f' -of VRT'
+ f' {WARP_OPTIONS}'
+ f' {GEBCO_ANTARCTIC_SUBICE_CLIP_VRT} {GEBCO_ANTARCTIC_SUBICE_FULL_VRT}')
run(msg,cmd)
-msg = 'Resampling Antarctic to target resolution...'
-cmd = f'gdalwarp -s_srs {EPSG_ANTARCTIC} -t_srs {EPSG_ANTARCTIC} -tr {RESOLUTION_M} {RESOLUTION_M} -r cubic -of GTiff {WARP_OPTIONS} {GDAL_CREATE_OPTIONS} {GEBCO_ANTARCTIC_FULL_VRT} {GEBCO_ANTARCTIC}'
+msg = 'Resampling Antarctic sub ice to target resolution...'
+cmd = (f'gdalwarp'
+ f' -s_srs {EPSG_ANTARCTIC}'
+ f' -t_srs {EPSG_ANTARCTIC}'
+ f' -tr {RESOLUTION_M} {RESOLUTION_M}'
+ f' -r cubic'
+ f' -of GTiff'
+ f' {WARP_OPTIONS}'
+ f' {GDAL_CREATE_OPTIONS}'
+ f' {GEBCO_ANTARCTIC_SUBICE_FULL_VRT} {GEBCO_ANTARCTIC_SUBICE}')
run(msg,cmd)
+# cleanup
+cleanup(GEBCO_WORLD_ICESURFACE_TMP1)
+cleanup(GEBCO_WORLD_ICESURFACE_TMP2)
cleanup(GEBCO_ARCTIC_CLIP_VRT)
cleanup(GEBCO_ARCTIC_FULL_VRT)
-cleanup(GEBCO_ANTARCTIC_CLIP_VRT)
-cleanup(GEBCO_ANTARCTIC_FULL_VRT)
+cleanup(GEBCO_ANTARCTIC_ICESURFACE_CLIP_VRT)
+cleanup(GEBCO_ANTARCTIC_ICESURFACE_FULL_VRT)
+cleanup(GEBCO_ANTARCTIC_ICESURFACE_TMP1)
+cleanup(GEBCO_ANTARCTIC_ICESURFACE_TMP2)
+cleanup(GEBCO_ANTARCTIC_SUBICE_CLIP_VRT)
+cleanup(GEBCO_ANTARCTIC_SUBICE_FULL_VRT)
diff --git a/process/5_create_shading_layers.py b/process/5_create_shading_layers.py
index 1c26dd6c07f4763d13835f90d10a29471d7d5b04..8fd6fe2e00177dcffe1015d3cdc4b9e742428530 100644
--- a/process/5_create_shading_layers.py
+++ b/process/5_create_shading_layers.py
@@ -11,9 +11,10 @@
This scipt computes the hillshade layers for the World, the Arctic
and the Antarctic basemap. The hillshades are computed with a
synthetic light source from 315° Azimuth and 45° Altitude and with
-and exaggeration / z factor of 10. The combined shading is used to
-get Simon's magical relief visualization, which combines shading
-by synthetic illumination and shading by slope inclination.
+and exaggeration / z factor of 10 (Z_FACTOR in config). The combined
+shading is used to get Simon's magical relief visualization, which
+combines shading by synthetic illumination and shading by slope
+inclination.
"""
@@ -23,8 +24,8 @@ by synthetic illumination and shading by slope inclination.
#=============================================================================
__author__ = 'Simon Dreutter'
-__version__ = '0.1'
-__date__ = '2020-12-17'
+__version__ = '0.3'
+__date__ = '2021-11-05'
__email__ = 'simon.dreutter@awi.de'
__status__ = 'Developement'
@@ -43,17 +44,131 @@ from config import *
# heading
print()
-print('============================================================')
-print(' COMPUTE SHADING LAYERS ')
+print('====================================================================')
+print(f'{"COMPUTE SHADING LAYERS":^68}')
+# World
msg = 'Creating World Hillshade...'
-cmd = f'gdaldem hillshade -b 1 -z {Z_FACTOR} -s 111120.0 -az 315.0 -alt 45.0 -compute_edges -combined -of GTiff {GDAL_CREATE_OPTIONS} {GEBCO_WORLD} {GEBCO_WORLD_HILLSHADE}'
+cmd = (f'gdaldem'
+ f' hillshade'
+ f' -b 1'
+ f' -z {Z_FACTOR}'
+ f' -s 111120.0'
+ f' -az 315.0'
+ f' -alt 45.0'
+ f' -compute_edges'
+ f' -combined'
+ f' -of GTiff'
+ f' {GDAL_CREATE_OPTIONS}'
+ f' {GEBCO_WORLD} {GEBCO_WORLD_HILLSHADE}')
run(msg,cmd)
+msg = 'Creating World ice surface Hillshade...'
+cmd = (f'gdaldem'
+ f' hillshade'
+ f' -b 1'
+ f' -z {Z_FACTOR}'
+ f' -s 111120.0'
+ f' -az 315.0'
+ f' -alt 45.0'
+ f' -compute_edges'
+ f' -combined'
+ f' -of GTiff'
+ f' {GDAL_CREATE_OPTIONS}'
+ f' {GEBCO_WORLD_ICESURFACE} {GEBCO_WORLD_ICESURFACE_HILLSHADE}')
+run(msg,cmd)
+
+msg = 'Creating World sub ice Hillshade...'
+cmd = (f'gdaldem'
+ f' hillshade'
+ f' -b 1'
+ f' -z {Z_FACTOR}'
+ f' -s 111120.0'
+ f' -az 315.0'
+ f' -alt 45.0'
+ f' -compute_edges'
+ f' -combined'
+ f' -of GTiff'
+ f' {GDAL_CREATE_OPTIONS}'
+ f' {GEBCO_WORLD_SUBICE} {GEBCO_WORLD_SUBICE_HILLSHADE}')
+run(msg,cmd)
+
+msg = 'Merging World Hillshades...'
+cmd = (f'gdal_merge'
+ f' -o {GEBCO_WORLD_HILLSHADE}'
+ f' -n 0'
+ f' -a_nodata 0'
+ f' {GEBCO_WORLD_SUBICE_HILLSHADE} {GEBCO_WORLD_ICESURFACE_HILLSHADE}')
+run(msg,cmd)
+
+# Arctic
msg = 'Creating Arctic Hillshade...'
-cmd = f'gdaldem hillshade -b 1 -z {Z_FACTOR} -s 1.0 -az 315.0 -alt 45.0 -compute_edges -combined -of GTiff {GDAL_CREATE_OPTIONS} {GEBCO_ARCTIC} {GEBCO_ARCTIC_HILLSHADE}'
+cmd = (f'gdaldem'
+ f' hillshade'
+ f' -b 1'
+ f' -z {Z_FACTOR}'
+ f' -s 1.0'
+ f' -az 315.0'
+ f' -alt 45.0'
+ f' -compute_edges'
+ f' -combined'
+ f' -of GTiff'
+ f' {GDAL_CREATE_OPTIONS}'
+ f' {GEBCO_ARCTIC} {GEBCO_ARCTIC_HILLSHADE}')
run(msg,cmd)
-msg = 'Creating Antarctic Hillshade...'
-cmd = f'gdaldem hillshade -b 1 -z {Z_FACTOR} -s 1.0 -az 315.0 -alt 45.0 -compute_edges -combined -of GTiff {GDAL_CREATE_OPTIONS} {GEBCO_ANTARCTIC} {GEBCO_ANTARCTIC_HILLSHADE}'
+# Antarctic
+# Due to ice shelfs there is a very steep drop in the DEM from shelf ice edge
+# down to the seabed, which creates hillshade artifacts. To avoid these, the
+# Antarctic GEBCO grid will processed twice, once the ice surface DEM and once
+# the sub ice topography DEM. The main hillshade is created from the sub ice
+# DEM. Then all data above sea level is extracted from the ice surface DEM, a
+# hillshade is created for these areas and the sub ice hillshade is replaced
+# in these areas.
+msg = 'Creating Antarctic ice surface Hillshade...'
+cmd = (f'gdaldem'
+ f' hillshade'
+ f' -b 1'
+ f' -z {Z_FACTOR}'
+ f' -s 1.0'
+ f' -az 315.0'
+ f' -alt 45.0'
+ f' -compute_edges'
+ f' -combined'
+ f' -of GTiff'
+ f' {GDAL_CREATE_OPTIONS}'
+ f' {GEBCO_ANTARCTIC_ICESURFACE} {GEBCO_ANTARCTIC_ICESURFACE_HILLSHADE}')
run(msg,cmd)
+
+msg = 'Creating Antarctic sub ice Hillshade...'
+cmd = (f'gdaldem'
+ f' hillshade'
+ f' -b 1'
+ f' -z {Z_FACTOR}'
+ f' -s 1.0'
+ f' -az 315.0'
+ f' -alt 45.0'
+ f' -compute_edges'
+ f' -combined'
+ f' -of GTiff'
+ f' {GDAL_CREATE_OPTIONS}'
+ f' {GEBCO_ANTARCTIC_SUBICE} {GEBCO_ANTARCTIC_SUBICE_HILLSHADE}')
+run(msg,cmd)
+
+msg = 'Merging Antarctic Hillshades...'
+cmd = (f'gdal_merge'
+ f' -o {GEBCO_ANTARCTIC_HILLSHADE}'
+ f' -n 0'
+ f' -a_nodata 0'
+ f' {GEBCO_ANTARCTIC_SUBICE_HILLSHADE} {GEBCO_ANTARCTIC_ICESURFACE_HILLSHADE}')
+run(msg,cmd)
+
+# cleanup
+cleanup(GEBCO_WORLD_ICESURFACE)
+cleanup(GEBCO_WORLD_ICESURFACE_HILLSHADE)
+cleanup(GEBCO_WORLD_SUBICE)
+cleanup(GEBCO_WORLD_SUBICE_HILLSHADE)
+cleanup(GEBCO_ANTARCTIC_ICESURFACE)
+cleanup(GEBCO_ANTARCTIC_ICESURFACE_HILLSHADE)
+cleanup(GEBCO_ANTARCTIC_SUBICE)
+cleanup(GEBCO_ANTARCTIC_SUBICE_HILLSHADE)
diff --git a/process/6_burn_vector_layers.py b/process/6_burn_vector_layers.py
index 3b508e65495889f3eeeba0dc42a15aaa91786413..f21133e426b24ae5fd03e565ec870bef1f382307 100644
--- a/process/6_burn_vector_layers.py
+++ b/process/6_burn_vector_layers.py
@@ -11,7 +11,7 @@
This script burns ths vector layers into the raster DEM layers. That
means, that raster cells are overwritten with the values 3000 (rock
outcrops) 20000 (shelf ice) and 30000 (ice sheets and glaciers)
-whereever they intersect with the vector layers.
+wherever they intersect with the vector layers.
"""
@@ -21,8 +21,8 @@ whereever they intersect with the vector layers.
#=============================================================================
__author__ = 'Simon Dreutter'
-__version__ = '0.1'
-__date__ = '2020-12-17'
+__version__ = '0.3'
+__date__ = '2021-11-05'
__email__ = 'simon.dreutter@awi.de'
__status__ = 'Developement'
@@ -57,53 +57,89 @@ V_ICE_SHEET = 30000.0
# heading
print()
-print('============================================================')
-print(' BURN VECTOR LAYERS ')
+print('====================================================================')
+print(f'{"BURN VECTOR LAYERS":^68}')
msg = 'Burning Antarctic ice shelfs in World raster...'
-cmd = f'gdal_rasterize -l {add_shelf_ice_name} -burn {V_SHELF_ICE} {ADD_SHELF_ICE} {GEBCO_WORLD}'
+cmd = (f'gdal_rasterize'
+ f' -l {add_shelf_ice_name}'
+ f' -burn {V_SHELF_ICE}'
+ f' {ADD_SHELF_ICE} {GEBCO_WORLD}')
run(msg,cmd)
msg = 'Burning Antarctic ice sheet in World raster...'
-cmd = f'gdal_rasterize -l {add_ice_sheet_name} -burn {V_ICE_SHEET} {ADD_ICE_SHEET} {GEBCO_WORLD}'
+cmd = (f'gdal_rasterize'
+ f' -l {add_ice_sheet_name}'
+ f' -burn {V_ICE_SHEET}'
+ f' {ADD_ICE_SHEET} {GEBCO_WORLD}')
run(msg,cmd)
msg = 'Burning Antarctic rock outcrops in World raster...'
-cmd = f'gdal_rasterize -l {add_rock_outcrop_name} -burn {V_ROCK_OUTCROPS} -at {ADD_ROCK_OUTCROP} {GEBCO_WORLD}'
+cmd = (f'gdal_rasterize'
+ f' -l {add_rock_outcrop_name}'
+ f' -burn {V_ROCK_OUTCROPS}'
+ f' -at'
+ f' {ADD_ROCK_OUTCROP} {GEBCO_WORLD}')
run(msg,cmd)
msg = 'Burning Greenland ice sheet in World raster...'
-cmd = f'gdal_rasterize -l {gimp_name} -burn {V_ICE_SHEET} {GIMP} {GEBCO_WORLD}'
+cmd = (f'gdal_rasterize'
+ f' -l {gimp_name}'
+ f' -burn {V_ICE_SHEET}'
+ f' {GIMP} {GEBCO_WORLD}')
run(msg,cmd)
msg = 'Burning glaciers in World raster...'
-cmd = f'gdal_rasterize -l {glims_world_name} -burn {V_ICE_SHEET} {GLIMS_WORLD} {GEBCO_WORLD}'
+cmd = (f'gdal_rasterize'
+ f' -l {glims_world_name}'
+ f' -burn {V_ICE_SHEET}'
+ f' {GLIMS_WORLD} {GEBCO_WORLD}')
run(msg,cmd)
msg = 'Burning Greenland ice sheet in Arctic raster...'
-cmd = f'gdal_rasterize -l {gimp_name} -burn {V_ICE_SHEET} {GIMP} {GEBCO_ARCTIC}'
+cmd = (f'gdal_rasterize'
+ f' -l {gimp_name}'
+ f' -burn {V_ICE_SHEET}'
+ f' {GIMP} {GEBCO_ARCTIC}')
run(msg,cmd)
msg = 'Burning glaciers in Arctic raster...'
-cmd = f'gdal_rasterize -l {glims_arctic_name} -burn {V_ICE_SHEET} {GLIMS_ARCTIC} {GEBCO_ARCTIC}'
+cmd = (f'gdal_rasterize'
+ f' -l {glims_arctic_name}'
+ f' -burn {V_ICE_SHEET}'
+ f' {GLIMS_ARCTIC} {GEBCO_ARCTIC}')
run(msg,cmd)
msg = 'Burning Antarctic ice shelfs in Antarctic raster...'
-cmd = f'gdal_rasterize -l {add_shelf_ice_name} -burn {V_SHELF_ICE} {ADD_SHELF_ICE} {GEBCO_ANTARCTIC}'
+cmd = (f'gdal_rasterize'
+ f' -l {add_shelf_ice_name}'
+ f' -burn {V_SHELF_ICE}'
+ f' {ADD_SHELF_ICE} {GEBCO_ANTARCTIC}')
run(msg,cmd)
msg = 'Burning Antarctic ice sheet in Antarctic raster...'
-cmd = f'gdal_rasterize -l {add_ice_sheet_name} -burn {V_ICE_SHEET} {ADD_ICE_SHEET} {GEBCO_ANTARCTIC}'
+cmd = (f'gdal_rasterize'
+ f' -l {add_ice_sheet_name}'
+ f' -burn {V_ICE_SHEET}'
+ f' {ADD_ICE_SHEET} {GEBCO_ANTARCTIC}')
run(msg,cmd)
msg = 'Burning Antarctic rock outcrops in Antarctic raster...'
-cmd = f'gdal_rasterize -l {add_rock_outcrop_name} -burn {V_ROCK_OUTCROPS} -at {ADD_ROCK_OUTCROP} {GEBCO_ANTARCTIC}'
+cmd = (f'gdal_rasterize'
+ f' -l {add_rock_outcrop_name}'
+ f' -burn {V_ROCK_OUTCROPS}'
+ f' -at'
+ f' {ADD_ROCK_OUTCROP} {GEBCO_ANTARCTIC}')
run(msg,cmd)
msg = 'Burning glaciers in Antarctic raster...'
-cmd = f'gdal_rasterize -l {glims_antarctic_name} -burn {V_ICE_SHEET} {GLIMS_ANTARCTIC} {GEBCO_ANTARCTIC}'
+cmd = (f'gdal_rasterize'
+ f' -l {glims_antarctic_name}'
+ f' -burn {V_ICE_SHEET}'
+ f' {GLIMS_ANTARCTIC} {GEBCO_ANTARCTIC}')
run(msg,cmd)
+# cleanup
cleanup(ADD_SHELF_ICE)
cleanup(ADD_ICE_SHEET)
cleanup(ADD_ROCK_OUTCROP)
diff --git a/process/7_create_basemaps.py b/process/7_create_basemaps.py
index a1205cc9baa975d274f47d017515b68199174713..f7b90142957f36bb81d7c5ced242c8e058a85a36 100644
--- a/process/7_create_basemaps.py
+++ b/process/7_create_basemaps.py
@@ -14,7 +14,7 @@ the DEM is rendered with the color palette into a RGB color raster.
Then the shaded map is calculated by combining the colored raster
with the shading raster.
At the end, the raster files are finalized by creating tiles and
-pyramids.
+pyramids and by adding metadata tags.
"""
@@ -24,8 +24,8 @@ pyramids.
#=============================================================================
__author__ = 'Simon Dreutter'
-__version__ = '0.1'
-__date__ = '2020-12-17'
+__version__ = '0.3'
+__date__ = '2021-11-05'
__email__ = 'simon.dreutter@awi.de'
__status__ = 'Developement'
@@ -34,6 +34,8 @@ __status__ = 'Developement'
# IMPORT
#=============================================================================
+import rasterio
+
# import AWI Basemap configuration
from config import *
@@ -53,29 +55,93 @@ OVERVIEW_METHOD = 'cubic'
# FUNCTIONS
#=============================================================================
-def create_basemap(map_params):
- name, gebco, hs, cp, epsg, output = map_params
+def create_basemap(name, grid, hillshade, color_palette, epsg, output):
+ """Create the AWI Basemap RGB grid from various inputs
+
+ Parameters
+ ----------
+ name : str
+ Name of grid version
+ grid : str
+ Path to input DEM grid
+ hillshade : str
+ Path to input hillshade grid
+ color_palette : str
+ Path to input color palette file
+ epsg : str
+ EPSG code
+ output : str
+ Path to output RGB grid
+
+ """
+ # create RGB version with grid and color palette
msg = f'Creating {name} RGB...'
- cmd = f'gdaldem color-relief {gebco} {cp} {TMP_RGB} -of GTiff -b 1 -compute_edges {GDAL_CREATE_OPTIONS}'
+ cmd = (f'gdaldem'
+ f' color-relief'
+ f' {grid}'
+ f' {color_palette}'
+ f' {TMP_RGB}'
+ f' -of GTiff'
+ f' -b 1'
+ f' -compute_edges'
+ f' {GDAL_CREATE_OPTIONS}')
run(msg,cmd)
+ # compile shaded and colored RGB grid
msg = f'Creating shaded {name} map...'
- cmd = f'{PYTHON} -m gdal_calc --calc "uint8(((A/255.)*(((B*0.75)+(255.*0.25))/255.))*255)" --format GTiff --type Byte -A {TMP_RGB} --A_band 1 -B {hs} --B_band 1 --allBands=A --overwrite --outfile {TMP_SHADED}'
+ cmd = (f'gdal_calc'
+ f' --calc "uint8(((A/255.)*(((B*0.75)+(255.*0.25))/255.))*254)"'
+ f' --format GTiff'
+ f' --type Byte'
+ f' -A {TMP_RGB}'
+ f' --A_band 1'
+ f' -B {hillshade}'
+ f' --B_band 1'
+ f' --allBands=A'
+ f' --overwrite'
+ f' --outfile {TMP_SHADED}')
# {GDAL_CREATE_OPTIONS_PY}
run(msg,cmd)
+ # cleanup
cleanup(TMP_RGB)
- msg = f'Finalizing AWI Basemap {name}...'
- cmd = f'gdal_translate -a_srs {epsg} -stats -of GTiff {GDAL_CREATE_OPTIONS_FINAL} {TMP_SHADED} {output}'
+ # finalize grid with GDAL create options
+ msg = f'Finalizing {name}...'
+ cmd = (f'gdal_translate'
+ f' -a_srs {epsg}'
+ f' -stats'
+ f' -of GTiff'
+ f' {GDAL_CREATE_OPTIONS_FINAL}'
+ f' {TMP_SHADED} {output}')
run(msg,cmd)
+ # cleanup
cleanup(TMP_SHADED)
cleanup(f'{TMP_SHADED}.aux.xml')
- msg = f'Creating pyramids for AWI Basemap {name}...'
- cmd = f'gdaladdo -r {OVERVIEW_METHOD} --config COMPRESS_OVERVIEW NONE --config INTERLEAVE_OVERVIEW PIXEL --config BIGTIFF_OVERVIEW IF_NEEDED --config GDAL_TIFF_OVR_BLOCKSIZE 512 {output} 2 4 8 16 32 64'
+ # create internal pyramid layers
+ msg = f'Creating pyramids for {name}...'
+ cmd = (f'gdaladdo'
+ f' -r {OVERVIEW_METHOD}'
+ f' --config COMPRESS_OVERVIEW NONE'
+ f' --config INTERLEAVE_OVERVIEW PIXEL'
+ f' --config BIGTIFF_OVERVIEW IF_NEEDED'
+ f' --config GDAL_TIFF_OVR_BLOCKSIZE 512'
+ f' {output}'
+ f' 2 4 8 16 32 64')
run(msg,cmd)
+
+ # set metadata in GeoTIFF tags
+ print('====================================================================')
+ print()
+ print(f'Updating Metadata for {name}...')
+ METADATA['NAME'] = name
+ METADATA['EPSG'] = epsg
+ with rasterio.open(output, 'r+') as ds:
+ # update metadata
+ ds.update_tags(**METADATA)
+ print()
#=============================================================================
@@ -84,24 +150,14 @@ def create_basemap(map_params):
# heading
print()
-print('============================================================')
-print(' CREATE BASEMAPS ')
-
-maps = [
- ('World [default]', GEBCO_WORLD, GEBCO_WORLD_HILLSHADE, CP_DEFAULT, EPSG_WORLD, BASEMAP_WORLD ),
- ('Arctic [default]', GEBCO_ARCTIC, GEBCO_ARCTIC_HILLSHADE, CP_DEFAULT, EPSG_ARCTIC, BASEMAP_ARCTIC ),
- ('Antarctic [default]', GEBCO_ANTARCTIC, GEBCO_ANTARCTIC_HILLSHADE, CP_DEFAULT, EPSG_ANTARCTIC, BASEMAP_ANTARCTIC ),
- ('World [bright]', GEBCO_WORLD, GEBCO_WORLD_HILLSHADE, CP_BRIGHT, EPSG_WORLD, BASEMAP_WORLD_BRIGHT ),
- ('Arctic [bright]', GEBCO_ARCTIC, GEBCO_ARCTIC_HILLSHADE, CP_BRIGHT, EPSG_ARCTIC, BASEMAP_ARCTIC_BRIGHT ),
- ('Antarctic [bright]', GEBCO_ANTARCTIC, GEBCO_ANTARCTIC_HILLSHADE, CP_BRIGHT, EPSG_ANTARCTIC, BASEMAP_ANTARCTIC_BRIGHT ),
- ('World [dark]', GEBCO_WORLD, GEBCO_WORLD_HILLSHADE, CP_DARK, EPSG_WORLD, BASEMAP_WORLD_DARK ),
- ('Arctic [dark]', GEBCO_ARCTIC, GEBCO_ARCTIC_HILLSHADE, CP_DARK, EPSG_ARCTIC, BASEMAP_ARCTIC_DARK ),
- ('Antarctic [dark]', GEBCO_ANTARCTIC, GEBCO_ANTARCTIC_HILLSHADE, CP_DARK, EPSG_ANTARCTIC, BASEMAP_ANTARCTIC_DARK )
-]
-
-for map_params in maps:
- create_basemap(map_params)
+print('====================================================================')
+print(f'{"CREATE BASEMAPS":^68}')
+
+# create the different basemaps with map parameters from config.py
+for kwargs in MAPS:
+ create_basemap(**kwargs)
+# cleanup
cleanup(GEBCO_WORLD)
cleanup(GEBCO_WORLD_HILLSHADE)
cleanup(GEBCO_ARCTIC)
@@ -111,9 +167,14 @@ cleanup(GEBCO_ANTARCTIC_HILLSHADE)
#=============================================================================
-# TODO
+# METADATA
#=============================================================================
-"""
-Metadata
-"""
+# write metadata to text file
+with open(METADATA_FILE, 'w') as f:
+ for key, val in METADATA.items():
+ if not key in ['NAME','EPSG']:
+ f.write(f'====================================================================\n')
+ f.write(f'{key.replace("_"," "):^68}\n')
+ f.write(f'====================================================================\n')
+ f.write(f'{val}\n\n\n')
diff --git a/process/RUN.py b/process/RUN.py
index 2eb531d7eb58493355de7074b05b00c215371ea5..da48da6a8e58c160fab234bace7df784f190c1e7 100644
--- a/process/RUN.py
+++ b/process/RUN.py
@@ -4,7 +4,6 @@
####################################################################
# #
# AWI Basemap #
-# RUN THE ENTIRE PROCESS #
# #
####################################################################
@@ -16,8 +15,8 @@
#=============================================================================
__author__ = 'Simon Dreutter'
-__version__ = '0.1'
-__date__ = '2021-05-19'
+__version__ = '0.3'
+__date__ = '2021-11-05'
__email__ = 'simon.dreutter@awi.de'
__status__ = 'Developement'
@@ -37,29 +36,29 @@ from config import PYTHON
# RUN
#=============================================================================
-cmd = [PYTHON, '1_process_add.py']
-subprocess.call(cmd)
-
-cmd = [PYTHON, '2_process_gimp.py']
-subprocess.call(cmd)
-
-cmd = [PYTHON, '3_process_glims.py']
-subprocess.call(cmd)
-
-cmd = [PYTHON, '4_process_gebco.py']
-subprocess.call(cmd)
-
-cmd = [PYTHON, '5_create_shading_layers.py']
-subprocess.call(cmd)
+print(__doc__)
+print('--------------------------------------------------------------------')
+print(f'version: {__version__}')
+print(f'author: {__author__}')
+print(f'email: {__email__}')
+print(f'status: {__status__}')
+print('--------------------------------------------------------------------')
+print()
-cmd = [PYTHON, '6_burn_vector_layers.py']
-subprocess.call(cmd)
+#=============================================================================
-cmd = [PYTHON, '7_create_basemaps.py']
-subprocess.call(cmd)
+subprocess.call([PYTHON, '1_process_add.py'])
+subprocess.call([PYTHON, '2_process_gimp.py'])
+subprocess.call([PYTHON, '3_process_glims.py'])
+subprocess.call([PYTHON, '4_process_gebco.py'])
+subprocess.call([PYTHON, '5_create_shading_layers.py'])
+subprocess.call([PYTHON, '6_burn_vector_layers.py'])
+subprocess.call([PYTHON, '7_create_basemaps.py'])
#=============================================================================
+print()
+print('--------------------------------------------------------------------')
print()
print('All done!')
input('\nPress Enter to exit!\n')
diff --git a/process/config.py b/process/config.py
index 7b123d4c109d2ef55f958bd3127d5b7cbaee0168..834d38e05307558804934fac9f50b8649e5db877 100644
--- a/process/config.py
+++ b/process/config.py
@@ -13,11 +13,12 @@ process. Adjust the settings below and the other basemap python
scripts will source these settings. The following can be adjusted
here:
+- Metadata
- Input folder paths
- Input file paths
- General create options (like resolution, extent, etc.)
- Paths to temporary files
-- Output file paths
+- Map versions (output files, names, color palettes)
- Operational parameters for the commands
- A small set of functions used by the other scripts
@@ -29,8 +30,8 @@ here:
#=============================================================================
__author__ = 'Simon Dreutter'
-__version__ = '0.1'
-__date__ = '2020-12-17'
+__version__ = '0.3'
+__date__ = '2021-11-05'
__email__ = 'simon.dreutter@awi.de'
__status__ = 'Developement'
@@ -39,9 +40,60 @@ __status__ = 'Developement'
# IMPORT
#=============================================================================
+import datetime as dt
import os
import sys
+from osgeo import gdal
+import rasterio
+
+
+#=============================================================================
+# METADATA
+#=============================================================================
+
+# metadata fields as convenient text blocks for easy writing
+ABSTRACT = f'''
+The AWI Basemap is a global basemap intended to be used as a background for GIS applications or web map viewers like the DAM marine data portal. The map is a rendered and shaded RGB version of the current GEBCO grid with ice overlays from ADD, GLIMS and GIMP. It is available in three projections, EPSG:4326 for the global map as well as EPSG:3995 and EPSG:3031 for the polar stereographic versions.
+AWI Basemaps are created by {__author__} ({__email__}), Bathymetry Group, Alfred Wegener Institute for Polar and Marine Research, Bremerhaven, Germany.
+'''
+
+SOURCES = '''
+GEBCO Compilation Group (2021) GEBCO 2021 Grid (doi:10.5285/c6612cbe-50b3-0cff-e053-6c86abc09f8f)
+GLIMS and NSIDC (2005, updated 2020): Global Land Ice Measurements from Space glacier database. Compiled and made available by the international GLIMS community and the National Snow and Ice Data Center, Boulder CO, U.S.A. DOI:10.7265/N5V98602
+Howat, I. 2017. MEaSUREs Greenland Ice Mapping Project (GIMP) Land Ice and Ocean Classification Mask, Version 1. [GimpIceMask_90m_v1.1]. Boulder, Colorado USA. NASA National Snow and Ice Data Center Distributed Active Archive Center. doi: https://doi.org/10.5067/B8X58MQBFUPA. [2020-11-19].
+Gerrish, L., Fretwell, P., & Cooper, P. (2021). Medium resolution vector polylines of the Antarctic coastline (7.4) [Data set]. UK Polar Data Centre, Natural Environment Research Council, UK Research & Innovation. https://doi.org/10.5285/824b5350-763e-4933-bb76-09f5d24cb033
+Gerrish, L., Fretwell, P., & Cooper, P. (2020). Medium resolution vector polygons of Antarctic rock outcrop (7.3) [Data set]. UK Polar Data Centre, Natural Environment Research Council, UK Research & Innovation. https://doi.org/10.5285/077e1f04-7068-4327-a4f2-71d863f70064
+'''
+
+LICENSE = '''
+CC BY 4.0
+'''
+
+USAGE_FEES = '''
+none
+'''
+
+CREATION_COMMENTS = f'''
+Python version {sys.version}
+GDAL version {gdal.__version__}
+Rasterio version {rasterio.__version__}
+AWI Basemap script version {__version__} ({__date__})
+grid creation date: {dt.date.today():%Y-%m-%d}
+'''
+
+#=============================================================================
+
+# metadata compiled as dictionary for later script usage
+METADATA = {
+ 'ABSTRACT' : ABSTRACT.strip(),
+ 'SOURCES' : SOURCES.strip(),
+ 'LICENSE' : LICENSE.strip(),
+ 'USAGE_FEES' : USAGE_FEES.strip(),
+ 'CREATION_COMMENTS' : CREATION_COMMENTS.strip()
+}
+
+
#=============================================================================
# INPUT FOLDERS
#=============================================================================
@@ -75,8 +127,9 @@ ADD_COASTLINE_IN = os.path.join(PATH_ADD,'add_coastline_medium_res_polygon_v7_4.
# input filepath to ADD Rock Outcrops
ADD_ROCK_OUTCROP_IN = os.path.join(PATH_ADD,'add_rock_outcrop_medium_res_polygon_v7.3.gpkg')
-# input filepath to GEBCO netCDF grid
-GEBCO_IN = os.path.join(PATH_GEBCO,'GEBCO_2020.nc')
+# input filepath to GEBCO netCDF grids
+GEBCO_IN = os.path.join(PATH_GEBCO,'GEBCO_2021.nc')
+GEBCO_SUBICE_IN = os.path.join(PATH_GEBCO,'GEBCO_2021_sub_ice_topo.nc')
# input filepath to GIMP GeoTIFF
GIMP_IN = os.path.join(PATH_GIMP,'GimpIceMask_90m_v1.1.tif')
@@ -84,11 +137,6 @@ GIMP_IN = os.path.join(PATH_GIMP,'GimpIceMask_90m_v1.1.tif')
# input filepath to GLIMS Esri Shapefile
GLIMS_IN = os.path.join(PATH_GLIMS,'glims_polygons.shp')
-# color palettes
-CP_DEFAULT = os.path.abspath(os.path.join('..','styles','color_palettes','basemap_default.txt'))
-CP_BRIGHT = os.path.abspath(os.path.join('..','styles','color_palettes','basemap_bright.txt'))
-CP_DARK = os.path.abspath(os.path.join('..','styles','color_palettes','basemap_dark.txt'))
-
#=============================================================================
# GENERAL CREATE OPTIONS
@@ -145,35 +193,37 @@ ANTARCTIC_EXTENT_LAT= -30.0
# scripts.
# ADD:
-# ADD land with cut rock outcrops
ADD_ICE_SHEET = os.path.join(PATH_ADD,'add_ice_sheet.gpkg')
-# ADD ice mask
ADD_SHELF_ICE = os.path.join(PATH_ADD,'add_shelf_ice.gpkg')
-# ADD rock outcrops
ADD_ROCK_OUTCROP = os.path.join(PATH_ADD,'add_rock_outcrop.gpkg')
#=============================================================================
# GIMP:
-# GIMP raster
GIMP_RASTER = os.path.join(PATH_GIMP,'GimpIceMask.tif')
-# GIMP vector
GIMP = os.path.join(PATH_GIMP,'GimpIceMask.gpkg')
#=============================================================================
# GLIMS:
-# GLIMS World
GLIMS_WORLD = os.path.join(PATH_GLIMS,'glims_world.gpkg')
-# GLIMS Arctic
GLIMS_ARCTIC = os.path.join(PATH_GLIMS,'glims_arctic.gpkg')
-# GLIMS Antarctic
GLIMS_ANTARCTIC = os.path.join(PATH_GLIMS,'glims_antarctic.gpkg')
#=============================================================================
# GEBCO:
-# GEBCO World
+# GEBCO World ice surface
+GEBCO_WORLD_ICESURFACE_TMP1 = os.path.join(PATH_GEBCO,'GEBCO_world_icesurface_tmp1.tif')
+GEBCO_WORLD_ICESURFACE_TMP2 = os.path.join(PATH_GEBCO,'GEBCO_world_icesurface_tmp2.tif')
+GEBCO_WORLD_ICESURFACE = os.path.join(PATH_GEBCO,'GEBCO_world_icesurface.tif')
+GEBCO_WORLD_ICESURFACE_HILLSHADE = os.path.join(PATH_GEBCO,'GEBCO_world_icesurface_hillshade.tif')
+
+# GEBCO World sub ice
+GEBCO_WORLD_SUBICE = os.path.join(PATH_GEBCO,'GEBCO_world_subice.tif')
+GEBCO_WORLD_SUBICE_HILLSHADE = os.path.join(PATH_GEBCO,'GEBCO_world_subice_hillshade.tif')
+
+# GEBCO World combined
GEBCO_WORLD = os.path.join(PATH_GEBCO,'GEBCO_world.tif')
GEBCO_WORLD_HILLSHADE = os.path.join(PATH_GEBCO,'GEBCO_world_hillshade.tif')
@@ -183,33 +233,29 @@ GEBCO_ARCTIC_FULL_VRT = os.path.join(PATH_GEBCO,'GEBCO_arctic_full.vrt')
GEBCO_ARCTIC = os.path.join(PATH_GEBCO,'GEBCO_arctic.tif')
GEBCO_ARCTIC_HILLSHADE = os.path.join(PATH_GEBCO,'GEBCO_arctic_hillshade.tif')
-# GEBCO Antarctic
-GEBCO_ANTARCTIC_CLIP_VRT = os.path.join(PATH_GEBCO,'GEBCO_antarctic_clip.vrt')
-GEBCO_ANTARCTIC_FULL_VRT = os.path.join(PATH_GEBCO,'GEBCO_antarctic_full.vrt')
+# GEBCO Antarctic ice surface
+GEBCO_ANTARCTIC_ICESURFACE_CLIP_VRT = os.path.join(PATH_GEBCO,'GEBCO_antarctic_icesurface_clip.vrt')
+GEBCO_ANTARCTIC_ICESURFACE_FULL_VRT = os.path.join(PATH_GEBCO,'GEBCO_antarctic_icesurface_full.vrt')
+GEBCO_ANTARCTIC_ICESURFACE_TMP1 = os.path.join(PATH_GEBCO,'GEBCO_antarctic_icesurface_tmp1.tif')
+GEBCO_ANTARCTIC_ICESURFACE_TMP2 = os.path.join(PATH_GEBCO,'GEBCO_antarctic_icesurface_tmp2.tif')
+GEBCO_ANTARCTIC_ICESURFACE = os.path.join(PATH_GEBCO,'GEBCO_antarctic_icesurface.tif')
+GEBCO_ANTARCTIC_ICESURFACE_HILLSHADE = os.path.join(PATH_GEBCO,'GEBCO_antarctic_icesurface_hillshade.tif')
+
+# GEBCO Antarctic sub ice
+GEBCO_ANTARCTIC_SUBICE_CLIP_VRT = os.path.join(PATH_GEBCO,'GEBCO_antarctic_subice_clip.vrt')
+GEBCO_ANTARCTIC_SUBICE_FULL_VRT = os.path.join(PATH_GEBCO,'GEBCO_antarctic_subice_full.vrt')
+GEBCO_ANTARCTIC_SUBICE = os.path.join(PATH_GEBCO,'GEBCO_antarctic_subice.tif')
+GEBCO_ANTARCTIC_SUBICE_HILLSHADE = os.path.join(PATH_GEBCO,'GEBCO_antarctic_subice_hillshade.tif')
+
+# Antarctic combined
GEBCO_ANTARCTIC = os.path.join(PATH_GEBCO,'GEBCO_antarctic.tif')
GEBCO_ANTARCTIC_HILLSHADE = os.path.join(PATH_GEBCO,'GEBCO_antarctic_hillshade.tif')
-# temporary files
-TMP_RGB = os.path.join(PATH_GEBCO,'tmp_rgb.tif')
-TMP_SHADED = os.path.join(PATH_GEBCO,'tmp_shaded.tif')
-
#=============================================================================
-# Final Output:
-# default colors
-BASEMAP_WORLD = os.path.join(PATH_RESULT,'AWI_Basemap_World.tif')
-BASEMAP_ARCTIC = os.path.join(PATH_RESULT,'AWI_Basemap_Arctic.tif')
-BASEMAP_ANTARCTIC = os.path.join(PATH_RESULT,'AWI_Basemap_Antarctic.tif')
-
-# bright colors
-BASEMAP_WORLD_BRIGHT = os.path.join(PATH_RESULT,'AWI_Basemap_World_bright.tif')
-BASEMAP_ARCTIC_BRIGHT = os.path.join(PATH_RESULT,'AWI_Basemap_Arctic_bright.tif')
-BASEMAP_ANTARCTIC_BRIGHT = os.path.join(PATH_RESULT,'AWI_Basemap_Antarctic_bright.tif')
-
-# dark colors
-BASEMAP_WORLD_DARK = os.path.join(PATH_RESULT,'AWI_Basemap_World_dark.tif')
-BASEMAP_ARCTIC_DARK = os.path.join(PATH_RESULT,'AWI_Basemap_Arctic_dark.tif')
-BASEMAP_ANTARCTIC_DARK = os.path.join(PATH_RESULT,'AWI_Basemap_Antarctic_dark.tif')
+# temporary files for final shading
+TMP_RGB = os.path.join(PATH_GEBCO,'tmp_rgb.tif')
+TMP_SHADED = os.path.join(PATH_GEBCO,'tmp_shaded.tif')
#=============================================================================
@@ -229,6 +275,109 @@ EPSG_ARCTIC = 'EPSG:3995' # WGS 84 / Arctic Polar Stereographic
EPSG_ANTARCTIC = 'EPSG:3031' # WGS 84 / Antarctic Polar Stereographic
+#=============================================================================
+# MAP VERSIONS
+#=============================================================================
+
+# In this section the different map versions are defined, based on the World,
+# Arctic and Antarctic files. The version defintion includes output, name and
+# color palette information.
+
+#=============================================================================
+
+# Final Output:
+# blue colors
+BASEMAP_WORLD_BLUE = os.path.join(PATH_RESULT,'AWI_Basemap_World_blue.tif')
+BASEMAP_ARCTIC_BLUE = os.path.join(PATH_RESULT,'AWI_Basemap_Arctic_blue.tif')
+BASEMAP_ANTARCTIC_BLUE = os.path.join(PATH_RESULT,'AWI_Basemap_Antarctic_blue.tif')
+
+# greyblue colors
+BASEMAP_WORLD_GREYBLUE = os.path.join(PATH_RESULT,'AWI_Basemap_World_greyblue.tif')
+BASEMAP_ARCTIC_GREYBLUE = os.path.join(PATH_RESULT,'AWI_Basemap_Arctic_greyblue.tif')
+BASEMAP_ANTARCTIC_GREYBLUE = os.path.join(PATH_RESULT,'AWI_Basemap_Antarctic_greyblue.tif')
+
+# grey colors
+BASEMAP_WORLD_GREY = os.path.join(PATH_RESULT,'AWI_Basemap_World_grey.tif')
+BASEMAP_ARCTIC_GREY = os.path.join(PATH_RESULT,'AWI_Basemap_Arctic_grey.tif')
+BASEMAP_ANTARCTIC_GREY = os.path.join(PATH_RESULT,'AWI_Basemap_Antarctic_grey.tif')
+
+#=============================================================================
+
+# color palettes
+CP_BLUE = os.path.abspath(os.path.join('..','styles','color_palettes','basemap_blue.txt'))
+CP_GREYBLUE = os.path.abspath(os.path.join('..','styles','color_palettes','basemap_greyblue.txt'))
+CP_GREY = os.path.abspath(os.path.join('..','styles','color_palettes','basemap_grey.txt'))
+
+#=============================================================================
+
+MAPS = [
+ {'name' : 'AWI Basemap World blue',
+ 'grid' : GEBCO_WORLD,
+ 'hillshade' : GEBCO_WORLD_HILLSHADE,
+ 'color_palette' : CP_BLUE,
+ 'epsg' : EPSG_WORLD,
+ 'output' : BASEMAP_WORLD_BLUE},
+
+ {'name' : 'AWI Basemap Arctic blue',
+ 'grid' : GEBCO_ARCTIC,
+ 'hillshade' : GEBCO_ARCTIC_HILLSHADE,
+ 'color_palette' : CP_BLUE,
+ 'epsg' : EPSG_ARCTIC,
+ 'output' : BASEMAP_ARCTIC_BLUE},
+
+ {'name' : 'AWI Basemap Antarctic blue',
+ 'grid' : GEBCO_ANTARCTIC,
+ 'hillshade' : GEBCO_ANTARCTIC_HILLSHADE,
+ 'color_palette' : CP_BLUE,
+ 'epsg' : EPSG_ANTARCTIC,
+ 'output' : BASEMAP_ANTARCTIC_BLUE},
+
+ {'name' : 'AWI Basemap World greyblue',
+ 'grid' : GEBCO_WORLD,
+ 'hillshade' : GEBCO_WORLD_HILLSHADE,
+ 'color_palette' : CP_GREYBLUE,
+ 'epsg' : EPSG_WORLD,
+ 'output' : BASEMAP_WORLD_GREYBLUE },
+
+ {'name' : 'AWI Basemap Arctic greyblue',
+ 'grid' : GEBCO_ARCTIC,
+ 'hillshade' : GEBCO_ARCTIC_HILLSHADE,
+ 'color_palette' : CP_GREYBLUE,
+ 'epsg' : EPSG_ARCTIC,
+ 'output' : BASEMAP_ARCTIC_GREYBLUE },
+
+ {'name' : 'AWI Basemap Antarctic greyblue',
+ 'grid' : GEBCO_ANTARCTIC,
+ 'hillshade' : GEBCO_ANTARCTIC_HILLSHADE,
+ 'color_palette' : CP_GREYBLUE,
+ 'epsg' : EPSG_ANTARCTIC,
+ 'output' : BASEMAP_ANTARCTIC_GREYBLUE },
+
+ {'name' : 'AWI Basemap World grey',
+ 'grid' : GEBCO_WORLD,
+ 'hillshade' : GEBCO_WORLD_HILLSHADE,
+ 'color_palette' : CP_GREY,
+ 'epsg' : EPSG_WORLD,
+ 'output' : BASEMAP_WORLD_GREY },
+
+ {'name' : 'AWI Basemap Arctic grey',
+ 'grid' : GEBCO_ARCTIC,
+ 'hillshade' : GEBCO_ARCTIC_HILLSHADE,
+ 'color_palette' : CP_GREY,
+ 'epsg' : EPSG_ARCTIC,
+ 'output' : BASEMAP_ARCTIC_GREY },
+
+ {'name' : 'AWI Basemap Antarctic grey',
+ 'grid' : GEBCO_ANTARCTIC,
+ 'hillshade' : GEBCO_ANTARCTIC_HILLSHADE,
+ 'color_palette' : CP_GREY,
+ 'epsg' : EPSG_ANTARCTIC,
+ 'output' : BASEMAP_ANTARCTIC_GREY },
+]
+
+METADATA_FILE = os.path.join(PATH_RESULT,'AWI_Basemap_Metadata.txt')
+
+
#=============================================================================
# OPERATIONAL PARAMETERS
#=============================================================================
@@ -257,13 +406,14 @@ GDAL_CREATE_OPTIONS_PY = '--co COMPRESS=NONE --co BIGTIFF=IF_NEEDED'
# GDAL warp options
WARP_OPTIONS = '-wo SOURCE_EXTRA=1000 -et 0'
+
#=============================================================================
# RUN FUNCTION
#=============================================================================
# function to run a shell command and print a message to stdout
def run(msg,cmd):
- print('============================================================')
+ print('====================================================================')
print()
print(msg)
print()
@@ -279,12 +429,28 @@ def run(msg,cmd):
# list of temporary files that shall not be deleted
DONT_DELETE = [
-# GEBCO_WORLD,
-# GEBCO_WORLD_HILLSHADE,
-# GEBCO_ARCTIC,
-# GEBCO_ARCTIC_HILLSHADE,
-# GEBCO_ANTARCTIC,
-# GEBCO_ANTARCTIC_HILLSHADE
+ #ADD_SHELF_ICE,
+ #ADD_ICE_SHEET,
+ #ADD_ROCK_OUTCROP,
+ #GIMP_RASTER,
+ #GIMP,
+ #GLIMS_WORLD,
+ #GLIMS_ARCTIC,
+ #GLIMS_ANTARCTIC,
+ #GEBCO_WORLD_ICESURFACE,
+ #GEBCO_WORLD_ICESURFACE_HILLSHADE,
+ #GEBCO_WORLD_SUBICE,
+ #GEBCO_WORLD_SUBICE_HILLSHADE,
+ #GEBCO_WORLD,
+ #GEBCO_WORLD_HILLSHADE,
+ #GEBCO_ARCTIC,
+ #GEBCO_ARCTIC_HILLSHADE,
+ #GEBCO_ANTARCTIC_ICESURFACE,
+ #GEBCO_ANTARCTIC_ICESURFACE_HILLSHADE,
+ #GEBCO_ANTARCTIC_SUBICE,
+ #GEBCO_ANTARCTIC_SUBICE_HILLSHADE,
+ #GEBCO_ANTARCTIC,
+ #GEBCO_ANTARCTIC_HILLSHADE
]
# function to cleanup temporary files if they are not in the DONT_DELETE list