From 96b90e0957c89aadaa25707ee0055349775fd85e Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Mon, 9 Jan 2023 14:56:58 +0100 Subject: [PATCH 1/3] Optimize planetry recipes for keeping highest resolution See #171 for background. I have solved the issue by determining the exact increment and tile size required for the highest resolution, which gets a name based on the floor of the actual increment. Need floor since 2.56 cannot be rint to 3 since we also will make a 03m grid. Testing this now but seems to work so want to get these into the repo for further testing. --- recipes/mars_relief.recipe | 5 +++-- recipes/mercury_relief.recipe | 3 ++- recipes/moon_relief.recipe | 5 +++-- recipes/pluto_relief.recipe | 3 ++- recipes/venus_relief.recipe | 12 ++++++++++-- scripts/get_prefix_func.sh | 17 ++++++++++------- scripts/srv_downsampler_grid.sh | 14 +++++++++++--- scripts/srv_downsampler_image.sh | 12 ++++++++++-- scripts/srv_tiler.sh | 20 +++++++++++++++----- 9 files changed, 66 insertions(+), 25 deletions(-) diff --git a/recipes/mars_relief.recipe b/recipes/mars_relief.recipe index 4da39c8..1669198 100644 --- a/recipes/mars_relief.recipe +++ b/recipes/mars_relief.recipe @@ -1,5 +1,5 @@ # Recipe file for down-filtering MOLA -# 2023-01-07 PW +# 2023-01-09 PW # # We use a precision of 0.5 m with a 6000 m offset to fit the range of -8528 to +21226 in 16-bit ints # @@ -28,8 +28,9 @@ # DST_SRTM=no # # List of desired output resolution and chunk size. Flag the source resolution with code == master +# Given dimension 106694 x 53347 we have spacing of 12.1468873601 arc seconds as master and 14x7 tiles # res unit tile chunk code -12 s 10 4096 master +12.1468873601 s 25.7142857143 4096 master 15 s 10 4096 30 s 15 4096 01 m 30 4096 diff --git a/recipes/mercury_relief.recipe b/recipes/mercury_relief.recipe index f8b05c1..36a9d07 100644 --- a/recipes/mercury_relief.recipe +++ b/recipes/mercury_relief.recipe @@ -28,8 +28,9 @@ # DST_SRTM=no # # List of desired output resolution and chunk size. Flag the source resolution with code == master +# Given dimension 23040 x 11520 we have spacing of 56.25 arc seconds as master and 12x6 tiles # res unit tile chunk code -56 s 30 4096 master +56.25 s 30 4096 master 01 m 30 4096 02 m 60 4096 03 m 90 2048 diff --git a/recipes/moon_relief.recipe b/recipes/moon_relief.recipe index 2dd4381..9f06e9e 100644 --- a/recipes/moon_relief.recipe +++ b/recipes/moon_relief.recipe @@ -1,5 +1,5 @@ # Recipe file for down-filtering LOLA -# 2023-01-03 PW +# 2023-01-09 PW # # We use a precision of 0.5 m with a zero offset to fit the range of -9128.5 to 10781.5 in 16-bit ints # @@ -28,8 +28,9 @@ # DST_SRTM=no # # List of desired output resolution and chunk size. Flag the source resolution with code == master +# Given dimension 92160 x 46080 we have spacing of 114.0625 arc seconds as master and 36x18 tiles # res unit tile chunk code -14 s 10 4096 master +14.0625 s 10 4096 master 15 s 10 4096 30 s 15 4096 01 m 30 4096 diff --git a/recipes/pluto_relief.recipe b/recipes/pluto_relief.recipe index 3779497..5574c1e 100644 --- a/recipes/pluto_relief.recipe +++ b/recipes/pluto_relief.recipe @@ -28,8 +28,9 @@ # DST_SRTM=no # # List of desired output resolution and chunk size. Flag the source resolution with code == master +# Given dimension 24888 x 12444 we have spacing of 52.07... arc seconds as master and 12x6 tiles # res unit tile chunk code -52 s 15 4096 master +52.0732883317 s 30 4096 master 01 m 30 4096 02 m 60 4096 03 m 90 2048 diff --git a/recipes/venus_relief.recipe b/recipes/venus_relief.recipe index c7fec7c..768ff08 100644 --- a/recipes/venus_relief.recipe +++ b/recipes/venus_relief.recipe @@ -1,5 +1,5 @@ # Recipe file for down-filtering Magellan -# 2023-01-07 PW +# 2023-01-09 PW # # We use a precision of 0.5 m with a 4000 m offset to fit the range of -2951 to 11687 in 16-bit ints # @@ -28,7 +28,15 @@ # DST_SRTM=no # # List of desired output resolution and chunk size. Flag the source resolution with code == master +# Given dimension 8192 x 4096 we have spacing of 2.63... arc minutes as master and 6x4 tiles # res unit tile chunk code -23 m 0 4096 master +2.63671875 m 60 4096 master +03 m 90 2048 +04 m 180 2048 +05 m 180 2048 +06 m 0 4096 +10 m 0 4096 +15 m 0 4096 +20 m 0 4096 30 m 0 4096 01 d 0 4096 diff --git a/scripts/get_prefix_func.sh b/scripts/get_prefix_func.sh index d0c9159..a7349ac 100644 --- a/scripts/get_prefix_func.sh +++ b/scripts/get_prefix_func.sh @@ -1,17 +1,20 @@ #!/bin/bash -# Function that creates tile label +# Function that creates tile label in integer degrees W/S function get_prefix () { # Takes west (in -180/180 range) and south and makes the {N|S}yy{W|E}xxx prefix - if [ $1 -ge 0 ]; then - X=$(printf "E%03d" $1) + # Get nearest integer degrees + W=$(gmt math -Q $1 RINT =) + S=$(gmt math -Q $2 RINT =) + if [ $W -ge 0 ]; then + X=$(printf "E%03d" $W) else - t=$(gmt math -Q $1 NEG =) + t=$(gmt math -Q $W NEG =) X=$(printf "W%03d" $t) fi - if [ $2 -ge 0 ]; then - Y=$(printf "N%02d" $2) + if [ $S -ge 0 ]; then + Y=$(printf "N%02d" $S) else - t=$(gmt math -Q $2 NEG =) + t=$(gmt math -Q $S NEG =) Y=$(printf "S%02d" $t) fi echo ${Y}${X} diff --git a/scripts/srv_downsampler_grid.sh b/scripts/srv_downsampler_grid.sh index 9302d0d..27c559e 100755 --- a/scripts/srv_downsampler_grid.sh +++ b/scripts/srv_downsampler_grid.sh @@ -10,6 +10,13 @@ # title, radius of the planetary body, desired node registration and resolutions, # desired output grid format and name prefix, and filter type, etc. Thus, this # script should handle data from different planets. +# Note: If the highest resolution grid is not an integer unit then some exploration +# needs to be done to determine what increment and tile size give an integer number +# of tiles over 360 and 180 ranges. E.g., below is the master line for mars_relief +# (which had 200 m pixels on Mars spheroid) and earth_relief (which as 15s exactly): +# 12.1468873601 s 25.7142857143 4096 master +# 15 s 10 4096 master +# Easiest to work with number of rows and find suitable common factors. if [ $# -eq 0 ]; then echo "usage: srv_downsampler_grid.sh recipefile" @@ -161,11 +168,12 @@ while read RES UNIT DST_TILE_SIZE CHUNK MASTER; do echo "Bad resolution $RES - aborting" exit -1 fi - if [ ! ${RES} = "01" ]; then # Use plural unit + IRES=$(gmt math -Q ${RES} FLOOR =) + if [ ${IRES} -gt 1 ]; then # Use plural unit UNIT_NAME="${UNIT_NAME}s" fi for REG in ${DST_NODES}; do # Probably doing both pixel and gridline registered output, except for master */ - DST_FILE=${DST_PLANET}/${DST_PREFIX}/${DST_PREFIX}_${RES}${UNIT}_${REG}.grd + DST_FILE=${DST_PLANET}/${DST_PREFIX}/${DST_PREFIX}_${IRES}${UNIT}_${REG}.grd grdtitle="${TITLE} at ${RES} arc ${UNIT_NAME}" # Note: The ${SRC_ORIG/+/\\+} below is to escape any plus-symbols in the file name with a backslash so grdedit -D will work if [ -f ${DST_FILE} ]; then # Do nothing @@ -177,7 +185,7 @@ while read RES UNIT DST_TILE_SIZE CHUNK MASTER; do remark="Reformatted from master file ${SRC_ORIG/+/\\+} [${REMARK}]" gmt grdedit ${DST_FILE} -D+t"${grdtitle}"+r"${remark}"+z"${SRC_NAME} (${SRC_UNIT})" fi - elif [ "${UNIT}" = "s" ] && [ ${RES} -le 30 ]; then # Special handling of xxs -> 15s or 30s filter due to 64-bit bug in grdfilter? + elif [ "${UNIT}" = "s" ] && [ ${IRES} -le 30 ]; then # Special handling of xxs -> 15s or 30s filter due to 64-bit bug in grdfilter? # See https://github.com/GenericMappingTools/remote-datasets/issues/32 - we do south and north hemisphere separately # Get suitable Gaussian full-width filter rounded to nearest 0.01 km after adding 50 meters for noise echo "Down-filter ${SRC_FILE} to ${DST_FILE}=${DST_MODIFY}" diff --git a/scripts/srv_downsampler_image.sh b/scripts/srv_downsampler_image.sh index 3e1b979..c955508 100755 --- a/scripts/srv_downsampler_image.sh +++ b/scripts/srv_downsampler_image.sh @@ -10,6 +10,13 @@ # format, radius of the planetary body, desired node registration and resolutions, # name prefix, and filter type, etc. Thus, this script should handle images from # different planets. +# Note: If the highest resolution image is not an integer unit then some exploration +# needs to be done to determine what increment and tile size give an integer number +# of tiles over 360 and 180 ranges. E.g., below is the master line for mars_relief +# (which had 200 m pixels on Mars spheroid) and earth_relief (which as 15s exactly): +# 12.1468873601 s 25.7142857143 4096 master +# 15 s 10 4096 master +# Easiest to work with number of rows and find suitable common factors. if [ $# -eq 0 ]; then echo "usage: srv_downsampler_image.sh recipefile" @@ -102,10 +109,11 @@ while read RES UNIT TILE MASTER; do echo "Bad resolution $RES - aborting" exit -1 fi - if [ ! ${RES} = "01" ]; then # Use plural unit + IRES=$(gmt math -Q ${RES} FLOOR =) + if [ ${IRES} -gt 1 ]; then # Use plural unit UNIT_NAME="${UNIT_NAME}s" fi - DST_FILE=${DST_PLANET}/${DST_PREFIX}/${DST_PREFIX}_${RES}${UNIT}.tif + DST_FILE=${DST_PLANET}/${DST_PREFIX}/${DST_PREFIX}_${IRES}${UNIT}.tif if [ -f ${DST_FILE} ]; then # Do nothing echo "${DST_FILE} exist - skipping" elif [ "X${MASTER}" = "Xmaster" ]; then # Just make a copy of the master to a new output file diff --git a/scripts/srv_tiler.sh b/scripts/srv_tiler.sh index a2f160f..558aa8a 100755 --- a/scripts/srv_tiler.sh +++ b/scripts/srv_tiler.sh @@ -16,6 +16,14 @@ # NOTE: We will ONLY look for the global files on this local machine. We first # look in the staging/ directory, and if not there then we look in the # users server directory. +# +# Note: If the highest resolution grid is not an integer unit then some exploration +# needs to be done to determine what increment and tile size give an integer number +# of tiles over 360 and 180 ranges. E.g., below is the master line for mars_relief +# (which had 200 m pixels on Mars spheroid) and earth_relief (which as 15s exactly): +# 12.1468873601 s 25.7142857143 4096 master +# 15 s 10 4096 master +# Easiest to work with number of rows and find suitable common factors. if [ $# -eq 0 ]; then echo "usage: srv_tiler.sh recipe [-f]" @@ -110,12 +118,13 @@ while read RES UNIT DST_TILE_SIZE CHUNK MASTER ; do fi for REG in ${DST_NODES}; do # Probably doing both pixel and gridline registered output, except for master */ # Name and path of grid we wish to tile - DST_TILE_TAG=${DST_PREFIX}_${RES}${UNIT}_${REG} + IRES=$(gmt math -Q ${RES} FLOOR =) + DST_TILE_TAG=${DST_PREFIX}_${IRES}${UNIT}_${REG} DST_FILE=${DST_TILE_TAG}.grd if [ -f ${DATADIR}/${DST_FILE} ]; then # found locally DATAGRID=${DATADIR}/${DST_FILE} else # Get it via local server files - DATAGRID=~/.gmt/server/${DST_PREFIX}_${RES}${UNIT}_${REG}.grd + DATAGRID=~/.gmt/server/${DST_PREFIX}_${IRES}${UNIT}_${REG}.grd fi if [ ! -f ${DATAGRID} ]; then # No echo "No such file to tile: ${DATAGRID}" @@ -126,7 +135,8 @@ while read RES UNIT DST_TILE_SIZE CHUNK MASTER ; do FILTER_WIDTH=$(gmt math -Q ${SRC_RADIUS} 2 MUL PI MUL 360 DIV $INC MUL 0.05 ADD 10 MUL RINT 10 DIV =) # Compute number of tiles required for this grid given nominal tile size. # We enforce square tiles by only solving for ny and doubling it for nx - if [ $DST_TILE_SIZE -gt 0 ]; then # OK, we need to split the file into separate tiles + IDST_TILE_SIZE=$(gmt math -Q ${DST_TILE_SIZE} RINT =) + if [ ${IDST_TILE_SIZE} -gt 0 ]; then # OK, we need to split the file into separate tiles ny=$(gmt math -Q 180 ${DST_TILE_SIZE} DIV =) nx=$(gmt math -Q ${ny} 2 MUL =) n_tiles=$(gmt math -Q $nx $ny MUL =) @@ -135,7 +145,7 @@ while read RES UNIT DST_TILE_SIZE CHUNK MASTER ; do # Build the list of w/e/s/n for the tiles gmt grdinfo ${DATAGRID} -I${DST_TILE_SIZE} -D -C > ${TMP}/wesn.txt # Determine local temporary tile directory for this product - TILEDIR=./${DST_PLANET}/${DST_PREFIX}/${DST_PREFIX}_${RES}${UNIT}_${REG} + TILEDIR=./${DST_PLANET}/${DST_PREFIX}/${DST_PREFIX}_${IRES}${UNIT}_${REG} rm -rf ${TILEDIR} mkdir -p ${TILEDIR} while read w e s n; do @@ -156,7 +166,7 @@ while read RES UNIT DST_TILE_SIZE CHUNK MASTER ; do MSG="${TITLE} at ${RES}x${RES} arc ${UNAME} reduced by Gaussian ${DST_MODE} filtering (${FILTER_WIDTH} km fullwidth)" fi printf "/server/%s/%s/\t%s_%s_%s/\t%s\t%s\t%s\t%s\t%4s\t%s\t%s\t-\t-\t%s\t%s [%s]\n" \ - ${DST_PLANET} ${DST_PREFIX} ${DST_PREFIX} ${TAG} ${REG} ${TAG} ${REG} ${DST_SCALE} ${DST_OFFSET} ${SIZE} ${DST_TILE_SIZE} ${creation_date} ${DST_CPT} "${MSG}" "${CITE}" >> ${DST_PREFIX}_server.txt + ${DST_PLANET} ${DST_PREFIX} ${DST_PREFIX} ${TAG} ${IREG} ${TAG} ${REG} ${DST_SCALE} ${DST_OFFSET} ${SIZE} ${DST_TILE_SIZE} ${creation_date} ${DST_CPT} "${MSG}" "${CITE}" >> ${DST_PREFIX}_server.txt # Move the tiled grid away from this tree mkdir -p ${TOPDIR}/staging/tiled From 24f4b2bbedc19c74e3eb23d6a2ef55bd4fab664a Mon Sep 17 00:00:00 2001 From: Federico Esteban Date: Mon, 9 Jan 2023 11:09:58 -0300 Subject: [PATCH 2/3] Update date --- recipes/pluto_relief.recipe | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/recipes/pluto_relief.recipe b/recipes/pluto_relief.recipe index 5574c1e..2b8daa3 100644 --- a/recipes/pluto_relief.recipe +++ b/recipes/pluto_relief.recipe @@ -1,5 +1,5 @@ # Recipe file for down-filtering New Horizon -# 2023-01-07 PW +# 2023-01-09 PW # # We use a precision of 0.25 m with a 1000 m offset to fit the range of -4101 to +6491 in 16-bit ints # From 57abfd7b41716bb941531892084cc2861e9eb901 Mon Sep 17 00:00:00 2001 From: Federico Esteban Date: Mon, 9 Jan 2023 11:10:49 -0300 Subject: [PATCH 3/3] Update mercury_relief.recipe --- recipes/mercury_relief.recipe | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/recipes/mercury_relief.recipe b/recipes/mercury_relief.recipe index 36a9d07..e764da0 100644 --- a/recipes/mercury_relief.recipe +++ b/recipes/mercury_relief.recipe @@ -1,5 +1,5 @@ # Recipe file for down-filtering Messenger -# 2023-01-07 PW +# 2023-01-09 PW # # We use a precision of 0.5 m with a zero offset to fit the range of -9128.5 to 10781.5 in 16-bit ints #