From f48b3d33ed3c30d11bb520fcbfddce75a425fde5 Mon Sep 17 00:00:00 2001
From: simondreutter <simon.dreutter@awi.de>
Date: Wed, 26 Jan 2022 11:18:41 +0100
Subject: [PATCH] process update 0.3
---
process/1_process_add.py | 45 ++++-
process/2_process_gimp.py | 23 ++-
process/3_process_glims.py | 27 ++-
process/4_process_gebco.py | 196 +++++++++++++++++++--
process/5_create_shading_layers.py | 137 +++++++++++++--
process/6_burn_vector_layers.py | 68 ++++++--
process/7_create_basemaps.py | 125 ++++++++++----
process/RUN.py | 41 +++--
process/config.py | 262 +++++++++++++++++++++++------
9 files changed, 759 insertions(+), 165 deletions(-)
diff --git a/process/1_process_add.py b/process/1_process_add.py
index 6fe33e8..722c17d 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 5e8a2f2..d258d64 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 191207b..b36249b 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 57d90b4..167b6dc 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 1c26dd6..8fd6fe2 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 3b508e6..f21133e 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 a1205cc..f7b9014 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 2eb531d..da48da6 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 7b123d4..834d38e 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
--
GitLab