Adjustments here and there.

Version v0.1

INSTRUCTIONS.md

+# Instructions to create the AWI Basemap
+
+Detailled instructions on how to create the AWI Basemap.
+
+> Note: All commands and scripts are expected to be executed in `~/process/`. If used otherwise, don't forget to adjust the paths.
+
+`Section 0` quickly describes the necessary software. `Section 1` lists the data that needs to be downloaded. `Section 2` describes all the necesaary adjustments to the scripts. And the rest is just run and wait.
+
+# 0. Software requirements
+
+The process described below is entirely based on open source GIS software. Most of the steps could also be done within the graphical user interface of QGIS. Others rely on commandline usage of the GDAL package (ogr2ogr) which comes with the QGIS installation. To run the scripts, you need the following:
+
+- python 3.x
+- gdal
+- gdal-python
+
+The scripts where designed for and tested with the OSGeo4W installation of the above mentioned packages.
+
+# 1. Download data
+
+Before processing and combining all input layers to the basemap, we need to download some geodata.
+
+##### General Bathymetric Chart of the Oceans (GEBCO)
+GEBCO is the DEM basis for the basemap. Download GEBCO in its global netCDF version to avoid the GeoTIFF tiles. Details on the GEBCO dataset can be found [here](https://gitlab.awi.de/sdreutte/basemap/-/blob/master/data/gebco.md).
+Download GEBCO here: https://www.gebco.net/data_and_products/gridded_bathymetry_data/
+The netCDF grid `GEBCO_2020.nc` should be placed in `~/data/gebco/`.
+
+##### Global Land Ice Measurements from Space (GLIMS)
+To overlay the map with global glaciers, GLIMS stores the necessary data. Download the latest GLIMS dataset as ESRI Shapefile Polygon. Details on the GLIMS dataset can be found [here](https://gitlab.awi.de/sdreutte/basemap/-/blob/master/data/glims.md).
+Download GLIMS here: http://www.glims.org/download/latest
+The SHP file `glims_polygons.shp` (and all auxiliary files) should be placed in `~data/glims/`.
+
+##### Greenland Ice Mapping Project (GIMP)
+To overlay the map with the Greenland ice sheet, GIMP stores the necessary data. Download the latest GIMP Ice Mask dataset as GeoTIFF. To access the GIMP dataset, a NASA Earthdata account is needed. Details on the GIMP dataset can be found [here](https://gitlab.awi.de/sdreutte/basemap/-/blob/master/data/gimp.md).
+Download GIMP here: https://nsidc.org/data/nsidc-0714/versions/1
+The geoTIFF grid `GimpIceMask_90m_v1.1.tif` should be placed in `~data/gimp/`.
+
+##### SCAR Antarctic Digital Database (ADD)
+To overlay the map with the Antarctic ice sheet and shelf ice, the ADD stores the necessary data. Download the coastline data and the rock outcrop data (medium resolution version is sufficient. Details on the ADD dataset can be found [here](https://gitlab.awi.de/sdreutte/basemap/-/blob/master/data/add.md).
+Download ADD coastlines here: https://data.bas.ac.uk/download/f18d638d4-3c71-442f-9415-01d17842d98d
+And ADD rock outcrops here: https://data.bas.ac.uk/download/fe90d6aec-b53e-40c1-ad52-c05e03a58c1d
+The GPKG files `add_coastline_medium_res_polygon_v7.3.gpkg` and `add_rock_outcrop_medium_res_polygon_v7.3.gpkg` should be placed in `~data/add/`.
+
+# 2. Adjust scripts
+
+Most of the necessary adjustments are happening in `~/process/config.py`. There are two minor things that need adjustment in other scripts.
+
+## 2.1 `~/process/config.py`
+
+##### Input Files
+After downloading the input data you need to make sure the input files in the configuration fit the filenames of the downloaded datasets. New versions of the datasets might lead to new filenames. Go through the `INPUT FILES` section in the configuration and check them all.
+
+##### General Create Options
+The default settings for resolution (degree/meter), the z factor for shading and the extent of the polar grids are set here. If you want different settings, adjust the variables in this section. However, read the notes, expecially on resolution.
+
+##### Operational Parameters
+If your Python 3.x installation has a different executable than OSGeo4W (`python3`), adjust this here.
+
+## 2.2 `~/process/1_process_add.py`
+
+Some of the vector processes require the GPKG layer names in their SQL commands. In ADD there is an inconsistency with the layer bames, which is why the layer name varibale for rock outcrops needs to be overwritten at the beginning of the process. In case this gets fixed or other incositencies arise in newer versions, make sure these are taken care of here.
+
+## 2.3 `~/process/run.bat`
+
+Again, if your Python 3.x executable is different to `python3`, adjust this at the beginning of the run script.
+
+# 3. Run
+
+Now run `~/process/run.bat` (or run all python scripts manually in order) and check the results, once finished.

README.md

@@ -8,65 +8,26 @@ The above mentioned websites use a global basemap with a colored and shaded digi
 
 At this point, the **AWI Basemap**s is available in three color sets. All three are based on the AWI blues (RGB bright: 7/172/231, RGB dark: 0/62/110). Yet, the "default" color palette reduces these blues to half their saturation, the "bright" color palette uses the blues directly, and the "dark" color palette uses the blues with a quarter of their saturation.
 
-
-<!---
 ## Instructions
 
-`INSTRUCTIONS.md` holds the instructions for manually creating the **AWI Basemap** with QGIS tools, as well as instructions on the usage of the scripts.
--->
+`INSTRUCTIONS.md` holds the instructions on the usage of the scripts to create the AWI Basemaps.
 
 ## Data
 
 In `data`, all the necessary data sources are listed. As data is not stored in this git, textfiles give information on where to download the data and in which format it needs to be available for the process.
 
+## Preview
+
+In `preview`, you can get a glimps on the final output.
+
 ## Process
 
 In `process`, all the required scripts are stored.
 
-
 ## Styles
 
-In `styles`, all the required color palettes and QGIS layer style files (QML) are stored.
+In `styles`, all the required color palettes and auxiliary QGIS layer style files (QML) are stored.
 
 ## Sources
 
 References and citations for data sources are listed in `SOURCES.md`.
-
-<!--
-## Results
-
-This is how the results look like:
-
-##### World
-
-Default:
-![AWI Basemap World Default](https://gitlab.awi.de/sdreutte/basemap/-/raw/master/images/AWI_Basemap_World.png)
-
-Bright:
-![AWI Basemap World Bright](https://gitlab.awi.de/sdreutte/basemap/-/raw/master/images/AWI_Basemap_World_bright.png)
-
-Dark:
-![AWI Basemap World Dark](https://gitlab.awi.de/sdreutte/basemap/-/raw/master/images/AWI_Basemap_World_dark.png)
-
-##### Arctic
-
-Default:
-![AWI Basemap Arctic Default](https://gitlab.awi.de/sdreutte/basemap/-/raw/master/images/AWI_Basemap_Arctic.png)
-
-Bright:
-![AWI Basemap Arctic Bright](https://gitlab.awi.de/sdreutte/basemap/-/raw/master/images/AWI_Basemap_Arctic_bright.png)
-
-Dark:
-![AWI Basemap Arctic Dark](https://gitlab.awi.de/sdreutte/basemap/-/raw/master/images/AWI_Basemap_Arctic_dark.png)
-
-##### Antarctic
-
-Default:
-![AWI Basemap Antarctic Default](https://gitlab.awi.de/sdreutte/basemap/-/raw/master/images/AWI_Basemap_Antarctic.png)
-
-Bright:
-![AWI Basemap Antarctic Bright](https://gitlab.awi.de/sdreutte/basemap/-/raw/master/images/AWI_BasemapAntarctic_bright.png)
-
-Dark:
-![AWI Basemap Antarctic Dark](https://gitlab.awi.de/sdreutte/basemap/-/raw/master/images/AWI_Basemap_Antarctic_dark.png)
--->
\ No newline at end of file

SOURCES.md

-**General Bathymetric Chart of the Oceans (GEBCO_2020)**
+**General Bathymetric Chart of the Oceans (GEBCO_2020)**  
 GEBCO Compilation Group (2020) GEBCO 2020 Grid (doi:10.5285/a29c5465-b138-234d-e053-6c86abc040b9)
 
-**Global Land Ice Measurements from Space (GLIMS version 20200630)**
+**Global Land Ice Measurements from Space (GLIMS version 20200630)**  
 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
 
-**Greenland Ice Mapping Project (GIMP version 1.1)**
+**Greenland Ice Mapping Project (GIMP version 1.1)**  
 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].
 
-**SCAR Antarctic Digital Database (ADD v7.3)**
-Gerrish, L., Fretwell, P., & Cooper, P. (2020). Medium resolution vector polygons of the Antarctic coastline (7.3) [Data set]. UK Polar Data Centre, Natural Environment Research Council, UK Research & Innovation. https://doi.org/10.5285/ed0a7b70-5adc-4c1e-8d8a-0bb5ee659d18
-
+**SCAR Antarctic Digital Database (ADD v7.3)**  
+Gerrish, L., Fretwell, P., & Cooper, P. (2020). Medium resolution vector polygons of the Antarctic coastline (7.3) [Data set]. UK Polar Data Centre, Natural Environment Research Council, UK Research & Innovation. https://doi.org/10.5285/ed0a7b70-5adc-4c1e-8d8a-0bb5ee659d18  
 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

preview/PREVIEW.md

+# Preview
+
+Here are some previews of the final results:
+
+##### World
+
+Default:
+![AWI Basemap World Default](https://gitlab.awi.de/sdreutte/basemap/-/raw/master/preview/AWI_Basemap_World.png)
+
+Bright:
+![AWI Basemap World Bright](https://gitlab.awi.de/sdreutte/basemap/-/raw/master/preview/AWI_Basemap_World_bright.png)
+
+Dark:
+![AWI Basemap World Dark](https://gitlab.awi.de/sdreutte/basemap/-/raw/master/preview/AWI_Basemap_World_dark.png)
+
+##### Arctic
+
+Default:
+![AWI Basemap Arctic Default](https://gitlab.awi.de/sdreutte/basemap/-/raw/master/preview/AWI_Basemap_Arctic.png)
+
+Bright:
+![AWI Basemap Arctic Bright](https://gitlab.awi.de/sdreutte/basemap/-/raw/master/preview/AWI_Basemap_Arctic_bright.png)
+
+Dark:
+![AWI Basemap Arctic Dark](https://gitlab.awi.de/sdreutte/basemap/-/raw/master/preview/AWI_Basemap_Arctic_dark.png)
+
+##### Antarctic
+
+Default:
+![AWI Basemap Antarctic Default](https://gitlab.awi.de/sdreutte/basemap/-/raw/master/preview/AWI_Basemap_Antarctic.png)
+
+Bright:
+![AWI Basemap Antarctic Bright](https://gitlab.awi.de/sdreutte/basemap/-/raw/master/preview/AWI_BasemapAntarctic_bright.png)
+
+Dark:
+![AWI Basemap Antarctic Dark](https://gitlab.awi.de/sdreutte/basemap/-/raw/master/preview/AWI_Basemap_Antarctic_dark.png)

process/1_process_add.py

@@ -40,22 +40,32 @@ from config import *
 
 
 #=============================================================================
-#    PROCESS
+#    GLOBAL VARIABLES
 #=============================================================================
 
+# projection of ADD data
 EPSG_ADD = 'EPSG:3031'
 
-# name of geometry column
-#geometry = 'geometry'
+
+#=============================================================================
+#    PROCESS
+#=============================================================================
 
 # layer names
 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]
-add_rock_outcrop_in_name = 'add_rock_outcrop_medium_res_polygon_v7' # inconsistancy with add 7.3 layer names
+add_rock_outcrop_in_name = os.path.splitext(os.path.basename(ADD_ROCK_OUTCROP_IN))[0]
+# inconsistancy with add 7.3 layer names, variable needs to be overwritten
+add_rock_outcrop_in_name = 'add_rock_outcrop_medium_res_polygon_v7'
 add_rock_outcrop_name = os.path.splitext(os.path.basename(ADD_ROCK_OUTCROP))[0]
 
+# heading
+print()
+print('============================================================')
+print('              ANTARCTIC DIGITAL DATANASE / ADD              ')
+
+
 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\'"'
 run(msg,cmd)

process/1_process_add_backup1.py

-# -*- coding: utf-8 -*-
-
-"""
-####################################################################
-#                                                                  #
-#                           AWI Basemap                            #
-#             PROCESS ANTARCTIC DIGITAL DATANASE DATA              #
-#                                                                  #
-####################################################################
-
-This script is processing the Antarctic Digital Database (ADD) input
-for the AWI Basemap. Input layers are the ADD coastlines and the ADD
-rock outcrops. Both layers are reprojected to a global CRS
-(EPSG_WORLD), coastlines are separated by shelf ice and ice sheet,
-and the rock outcrop geometries are fixed.
-
-"""
-
-
-#=============================================================================
-#    SCRIPT INFO
-#=============================================================================
-
-__author__ = 'Simon Dreutter'
-__version__ = '0.1'
-__date__ = '2020-12-17'
-__email__ = 'simon.dreutter@awi.de'
-__status__ = 'Developement'
-
-
-#=============================================================================
-#    IMPORT
-#=============================================================================
-
-# import AWI Basemap configuration
-from config import *
-
-
-#=============================================================================
-#    PROCESS
-#=============================================================================
-
-EPSG_ADD = 'EPSG:3031'
-
-msg = 'Reprojecting and extracting land mask from ADD coastlines...'
-cmd = f'ogr2ogr -s_srs {EPSG_ADD} -t_srs {EPSG_WORLD} -nln "add_coastline_land" -where "\"surface\" = \'land\'" -wrapdateline -f GPKG {ADD_COASTLINE_LAND} {ADD_COASTLINE_IN}'
-run(msg,cmd)
-
-msg = 'Reprojecting and extracting ice mask from ADD coastlines...'
-cmd = f'ogr2ogr -s_srs {EPSG_ADD} -t_srs {EPSG_WORLD} -nln "add_coastline_ice" -where "\"surface\" = \'ice shelf\' OR \"surface\" = \'rumple\' OR \"surface\" = \'ocean\' OR \"surface\" = \'ice tongue\'" -wrapdateline -f GPKG {ADD_SHELF_ICE} {ADD_COASTLINE_IN}'
-run(msg,cmd)
-
-msg = 'Reprojecting ADD rock outcrops...'
-cmd = f'ogr2ogr -s_srs {EPSG_ADD} -t_srs {EPSG_WORLD} -nln "add_rock_outcrop" -wrapdateline -f GPKG {ADD_ROCK_OUTCROP} {ADD_ROCK_OUTCROP_IN}'
-run(msg,cmd)
-
-msg = 'Fixing ADD rock outcrops geometries...'
-cmd = f'ogr2ogr -dialect sqlite -sql "SELECT ST_Buffer(geom, 0.0) AS geom,* FROM \"add_rock_outcrop\"" -nln "add_rock_outcrop_fixed" -f GPKG {ADD_ROCK_OUTCROP_FIXED} {ADD_ROCK_OUTCROP}'
-run(msg,cmd)
-
-print()
-print('ADD done!')
-print()
-print('Now do not forget to manually cut the rock outcrops from the ADD land layer!')
-print()
-input('Press Enter to exit!')
-print()
-
-
-#=============================================================================
-#    TODO
-#=============================================================================
-
-"""
-- find a scripted solution for 'Difference'
-- once found, add cleanup
-"""
\ No newline at end of file

process/1_process_add_backup2.py

-# -*- coding: utf-8 -*-
-
-"""
-####################################################################
-#                                                                  #
-#                           AWI Basemap                            #
-#             PROCESS ANTARCTIC DIGITAL DATANASE DATA              #
-#                                                                  #
-####################################################################
-
-This script is processing the Antarctic Digital Database (ADD) input
-for the AWI Basemap. Input layers are the ADD coastlines and the ADD
-rock outcrops. Both layers are reprojected to a global CRS
-(EPSG_WORLD), coastlines are separated by shelf ice and ice sheet,
-and the rock outcrop geometries are fixed.
-
-"""
-
-
-#=============================================================================
-#    SCRIPT INFO
-#=============================================================================
-
-__author__ = 'Simon Dreutter'
-__version__ = '0.1'
-__date__ = '2020-12-17'
-__email__ = 'simon.dreutter@awi.de'
-__status__ = 'Developement'
-
-
-#=============================================================================
-#    IMPORT
-#=============================================================================
-
-# import standard libraries
-import os
-
-# import AWI Basemap configuration
-from config import *
-
-
-#=============================================================================
-#    PROCESS
-#=============================================================================
-
-EPSG_ADD = 'EPSG:3031'
-
-# name of geometry column
-geometry = 'geometry'
-
-# layer names
-add_coastline_in_name = os.path.splitext(os.path.basename(ADD_COASTLINE_IN))[0]
-add_coastline_land_name = os.path.splitext(os.path.basename(ADD_COASTLINE_LAND))[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]
-add_rock_outcrop_in_name = 'add_rock_outcrop_medium_res_polygon_v7' # inconsistancy with add layer names
-add_rock_outcrop_name = os.path.splitext(os.path.basename(ADD_ROCK_OUTCROP))[0]
-add_rock_outcrop_fixed_name = os.path.splitext(os.path.basename(ADD_ROCK_OUTCROP_FIXED))[0]
-
-msg = 'Reprojecting and extracting land mask from ADD coastlines...'
-cmd = f'ogr2ogr -s_srs {EPSG_ADD} -t_srs {EPSG_WORLD} -wrapdateline -nln "{add_coastline_land_name}" -nlt  MULTIPOLYGON -f "ESRI Shapefile" -overwrite {ADD_COASTLINE_LAND} {ADD_COASTLINE_IN} -dialect SQLite -sql "SELECT ST_Multi(geom) 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_Multi(geom) from \'{add_coastline_in_name}\' WHERE surface IN (\'ice shelf\',\'rumple\',\'ocean\',\'ice tongue\')"'
-run(msg,cmd)
-
-msg = 'Reprojecting and fixing ADD rock outcrops...'
-cmd = f'ogr2ogr -s_srs {EPSG_ADD} -t_srs {EPSG_WORLD} -wrapdateline -nln "{add_rock_outcrop_name}" -nlt  MULTIPOLYGON -f "ESRI Shapefile" -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}\'"'
-run(msg,cmd)
-
-msg = 'Cutting rock outcrops from ice sheet...'
-cmd = f'ogr2ogr -dialect SQLite -sql "select ST_Difference(a.{geometry}, ST_Simplify(b.{geometry},500)) as {geometry} from {add_coastline_land_name} as a, \'{ADD_ROCK_OUTCROP}\'.{add_rock_outcrop_name} as b" -overwrite -f GPKG {ADD_ICE_SHEET} {ADD_COASTLINE_LAND}'
-run(msg,cmd)
-
-
-#msg = 'Fixing ADD rock outcrops geometries...'
-#cmd = f'ogr2ogr -dialect sqlite -sql "SELECT ST_Buffer({geometry}, 0.0) AS {geometry},* FROM \"{add_rock_outcrop_name}\"" -nln "{add_rock_outcrop_fixed_name}" -overwrite -f "ESRI Shapefile" {ADD_ROCK_OUTCROP_FIXED} {ADD_ROCK_OUTCROP}'
-#run(msg,cmd)
-
-#msg = 'Cutting rock outcrops from ice sheet...'
-#cmd = f'ogr2ogr -dialect SQLite -sql "select ST_Difference(a.{geometry}, ST_Simplify(b.{geometry},500)) as {geometry} from {add_coastline_land_name} as a, \'{ADD_ROCK_OUTCROP_FIXED}\'.{add_rock_outcrop_fixed_name} as b" -overwrite -f "ESRI Shapefile" {ADD_ICE_SHEET} {ADD_COASTLINE_LAND}'
-#run(msg,cmd)
-
-
-
-#=============================================================================
-#    TODO
-#=============================================================================
-
-"""
-- find a scripted solution for 'Difference'
-- once found, add cleanup
-"""
\ No newline at end of file

process/2_process_gimp.py

@@ -44,9 +44,13 @@ from config import *
 # layer names
 gimp_name = os.path.splitext(os.path.basename(GIMP))[0]
 
+# heading
+print()
+print('============================================================')
+print('            GREENLAND ICE MAPPING PROJECT / GIMP            ')
+
 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}'
-#{GDAL_CREATE_OPTIONS}
 run(msg,cmd)
 
 msg = 'Polygonizing GIMP raster...'

process/3_process_glims.py

@@ -10,8 +10,8 @@
 
 This script is processing the Global Land Ice Measurements from
 Space (GLIMS) input for the AWI Basemap. The GLIMS geometries are
-fixed and the global dataset is clipped for the Arctic and the
-Antarctic.
+fixed (with zero buffer) and the global dataset is clipped for the
+Arctic and the Antarctic.
 
 """
 
@@ -48,6 +48,11 @@ glims_world_name = os.path.splitext(os.path.basename(GLIMS_WORLD))[0]
 glims_arctic_name = os.path.splitext(os.path.basename(GLIMS_ARCTIC))[0]
 glims_antarctic_name = os.path.splitext(os.path.basename(GLIMS_ANTARCTIC))[0]
 
+# heading
+print()
+print('============================================================')
+print('      GLOBAL LAND ICE MEASUREMENTS FROM SPACE / GLIMS       ')
+
 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"'
 run(msg,cmd)

process/4_process_gebco.py

@@ -42,49 +42,54 @@ from config import *
 #    PROCESS
 #=============================================================================
 
+# number of pixels for interpolation (fill NoData at poles)
+interpolation = 1 / RESOLUTION_M * 5000         # TEMPORARY
+
+# heading
+print()
+print('============================================================')
+print('      GENERAL BATHYMETRIC CHART OF THE OCEANS / GEBCO       ')
+
 msg = 'Resampling GEBCO...'
-cmd = f'gdalwarp {WARP_OPTIONS} -s_srs {EPSG_WORLD} -t_srs {EPSG_WORLD} -tr {RESOLUTION_DEG} {RESOLUTION_DEG} -r cubic -of GTiff {GEBCO_IN} {GEBCO_WORLD}'
-# {GDAL_CREATE_OPTIONS}
-#run(msg,cmd)
+cmd = f'gdalwarp -s_srs {EPSG_WORLD} -t_srs {EPSG_WORLD} -tr {RESOLUTION_DEG} {RESOLUTION_DEG} -r cubic -of GTiff {GDAL_CREATE_OPTIONS} {GEBCO_IN} {GEBCO_WORLD}'
+run(msg,cmd)
 
 msg = 'Clipping GEBCO for the Arctic...'
-cmd = f'gdal_translate -projwin -180.0 90.0 180.0 {ARCTIC_EXTENT_LAT} -of VRT {GEBCO_IN} {GEBCO_ARCTIC_CLIP_VRT}'
+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}'
 run(msg,cmd)
 
 msg = 'Reprojecting GEBCO for the Arctic...'
 cmd = f'gdalwarp -s_srs {EPSG_WORLD} -t_srs {EPSG_ARCTIC} -r near -of VRT {GEBCO_ARCTIC_CLIP_VRT} {GEBCO_ARCTIC_FULL_VRT}'
-#{WARP_OPTIONS}
 run(msg,cmd)
 
 msg = 'Resampling Arctic to target resolution...'
-cmd = f'gdalwarp {WARP_OPTIONS} -s_srs {EPSG_ARCTIC} -t_srs {EPSG_ARCTIC} -tr {RESOLUTION_M} {RESOLUTION_M} -r cubic -of GTiff {GEBCO_ARCTIC_FULL_VRT} {TMP_POLAR}'
+cmd = f'gdalwarp -s_srs {EPSG_ARCTIC} -t_srs {EPSG_ARCTIC} -tr {RESOLUTION_M} {RESOLUTION_M} -r cubic -of GTiff {GDAL_CREATE_OPTIONS} {GEBCO_ARCTIC_FULL_VRT} {TMP_POLAR}'
 #{GEBCO_ARCTIC}
 run(msg,cmd)
 
 # TEMPORARY
 msg = 'Filling North Pole gap...'
-cmd = f'{PYTHON} -m gdal_fillnodata -md {INTERPOLATION} -b 1 -of GTiff {TMP_POLAR} {GEBCO_ARCTIC}'
+cmd = f'{PYTHON} -m gdal_fillnodata -md {interpolation} -b 1 -of GTiff {TMP_POLAR} {GEBCO_ARCTIC}'
 run(msg,cmd)
 cleanup(TMP_POLAR)
 # /TEMPORARY
 
 msg = 'Clipping GEBCO for the Antarctic...'
-cmd = f'gdal_translate -projwin -180.0 {ANTARCTIC_EXTENT_LAT} 180.0 -90.0 -of VRT {GEBCO_IN} {GEBCO_ANTARCTIC_CLIP_VRT}'
+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}'
 run(msg,cmd)
 
 msg = 'Reprojecting GEBCO for the Antarctic...'
 cmd = f'gdalwarp -s_srs {EPSG_WORLD} -t_srs {EPSG_ANTARCTIC} -r near -of VRT {GEBCO_ANTARCTIC_CLIP_VRT} {GEBCO_ANTARCTIC_FULL_VRT}'
-#{WARP_OPTIONS}
 run(msg,cmd)
 
 msg = 'Resampling Antarctic to target resolution...'
-cmd = f'gdalwarp {WARP_OPTIONS} -s_srs {EPSG_ANTARCTIC} -t_srs {EPSG_ANTARCTIC} -tr {RESOLUTION_M} {RESOLUTION_M} -r cubic -of GTiff {GEBCO_ANTARCTIC_FULL_VRT} {TMP_POLAR}'
+cmd = f'gdalwarp -s_srs {EPSG_ANTARCTIC} -t_srs {EPSG_ANTARCTIC} -tr {RESOLUTION_M} {RESOLUTION_M} -r cubic -of GTiff {GDAL_CREATE_OPTIONS} {GEBCO_ANTARCTIC_FULL_VRT} {TMP_POLAR}'
 #{GEBCO_ANTARCTIC}
 run(msg,cmd)
 
 # TEMPORARY
 msg = 'Filling South Pole gap...'
-cmd = f'{PYTHON} -m gdal_fillnodata -md {INTERPOLATION} -b 1 -of GTiff {TMP_POLAR} {GEBCO_ANTARCTIC}'
+cmd = f'{PYTHON} -m gdal_fillnodata -md {interpolation} -b 1 -of GTiff {TMP_POLAR} {GEBCO_ANTARCTIC}'
 run(msg,cmd)
 cleanup(TMP_POLAR)
 # /TEMPORARY
@@ -101,13 +106,19 @@ cleanup(GEBCO_ANTARCTIC_FULL_VRT)
 
 """
 The current solution for creating the polar stereographic GEBCO grids (VRT)
-produces holes on the north and south pole. Another solution needs to be
-found. Normal clip and warp however produces issues:
+produces holes on the north and south pole, which are then interpolated. This
+likely comes from creating polar stereographic grids with such large extents
+(low latitudes). Everything larger than down to 70°N / up to 70°S creates a
+hole.
+The interpolation (gdal_fillnodata), however, extrapolates
+the grid around the minmax latitude circle.
+Another solution needs to be found. Normal clip and warp, however, produces
+issues:
 
 - Nearest Neighbor: raster artifacts
 - Cubic (and others): dateline blur
 
 What to do?
 
-- get rid of interpolation (because of the ugly extrapolation)
-"""
\ No newline at end of file
+- somehow get rid of interpolation
+"""

process/5_create_shading_layers.py

@@ -4,11 +4,11 @@
 ####################################################################
 #                                                                  #
 #                           AWI Basemap                            #
-#                      CREATE SHADING LAYERS                       #
+#                      COMPUTE SHADING LAYERS                      #
 #                                                                  #
 ####################################################################
 
-This scipt creates the hillshade layers for the World, the Arctic
+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
@@ -41,14 +41,19 @@ from config import *
 #    PROCESS
 #=============================================================================
 
+# heading
+print()
+print('============================================================')
+print('                   COMPUTE SHADING LAYERS                   ')
+
 msg = 'Creating World Hillshade...'
-cmd = f'gdaldem hillshade {GEBCO_WORLD} {GEBCO_WORLD_HILLSHADE} -of GTiff -b 1 -z {Z_FACTOR} -s 111120.0 -az 315.0 -alt 45.0 -compute_edges -combined {GDAL_CREATE_OPTIONS}'
+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}'
 run(msg,cmd)
 
 msg = 'Creating Arctic Hillshade...'
-cmd = f'gdaldem hillshade {GEBCO_ARCTIC} {GEBCO_ARCTIC_HILLSHADE} -of GTiff -b 1 -z {Z_FACTOR} -s 1.0 -az 315.0 -alt 45.0 -compute_edges -combined {GDAL_CREATE_OPTIONS}'
+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}'
 run(msg,cmd)
 
 msg = 'Creating Antarctic Hillshade...'
-cmd = f'gdaldem hillshade {GEBCO_ANTARCTIC} {GEBCO_ANTARCTIC_HILLSHADE} -of GTiff -b 1 -z {Z_FACTOR} -s 1.0 -az 315.0 -alt 45.0 -compute_edges -combined {GDAL_CREATE_OPTIONS}'
+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}'
 run(msg,cmd)