search

GenericMappingTools / gmtserver-admin

Cache data and script for managing the GMT data server
GNU Lesser General Public License v3.0
7 stars 3 forks source link

Changes to @earth_day and @earth_night grids from GMT 6.4 to 6.5? #257

Open weiji14 opened 6 months ago

weiji14 commented 6 months ago

Just checking something that I noticed at https://github.com/GenericMappingTools/pygmt/pull/2963#discussion_r1444072938. Are the Blue Marble @earth_day_01d_p and Black Marble @earth_night_01d_p datasets only downloadable with GMT 6.5.0 (and not GMT 6.4.0 or older) now?

E.g. doing gmt which @earth_day_01d_p -Vd with GMT 6.4.0 gives:

gmt [DEBUG]: GMT_Create_Session: Terminal width = 190
gmt [DEBUG]: Obtained the ppid from parent: 3763
gmt [DEBUG]: Enter: gmtinit_new_GMT_ctrl
gmt [DEBUG]: GMT->session.SHAREDIR = ~/mambaforge/envs/pygmt/share/gmt
gmt [DEBUG]: GMT->session.HOMEDIR = ~
gmt [DEBUG]: GMT->session.USERDIR = ~/.gmt [created]
gmt [DEBUG]: GMT->session.CACHEDIR = ~/.gmt/cache [created]
gmt [DEBUG]: GMT: 0. Will try to find subdir=postscriptlight stem = PSL_custom_fonts suffix=.txt
gmt [DEBUG]: GMT: 1. gmt_getsharepath trying current dir
gmt [DEBUG]: GMT: 2. gmt_getsharepath trying USERDIR ~/.gmt
gmt [DEBUG]: GMT: 3. gmt_getsharepath trying USERDIR subdir ~/.gmt/postscriptlight
gmt [DEBUG]: GMT: 4. gmt_getsharepath trying SHAREDIR subdir ~/mambaforge/envs/pygmt/share/gmt/postscriptlight
gmt [DEBUG]: GMT: 5. gmt_getsharepath trying SHAREDIR ~/mambaforge/envs/pygmt/share/gmt
gmt [DEBUG]: GMT: 6. gmt_getsharepath failed
gmt [DEBUG]: Map distance calculation will be Cartesian
gmt [DEBUG]: Exit:  gmtinit_new_GMT_ctrl
gmt [DEBUG]: Enter: New_PSL_Ctrl
gmt [DEBUG]: Exit:  New_PSL_Ctrl
gmt [DEBUG]: Enter: gmt_manage_workflow
gmt [DEBUG]: Exit : gmt_manage_workflow
gmt [DEBUG]: Enter: PSL_beginsession
gmt [DEBUG]: Exit : PSL_beginsession
gmt [DEBUG]: Enter: PSL_setdefaults
gmt [DEBUG]: Exit : PSL_setdefaults
gmt [DEBUG]: Enter: gmtlib_io_init
gmt [DEBUG]: Exit : gmtlib_io_init
gmt [DEBUG]: Enter: gmt_hash_init
gmt [DEBUG]: Exit:  gmt_hash_init
gmt [DEBUG]: Enter: gmt_hash_init
gmt [DEBUG]: Exit:  gmt_hash_init
gmt [DEBUG]: Enter: gmt_reload_settings
gmt [DEBUG]: The PROJ_GEODESIC set to Vincenty
gmt [DEBUG]: Look for file ~/gmt.conf
gmt [DEBUG]: Look for file ~/.gmt/gmt.conf
gmt [DEBUG]: Look for file ~/.gmt/server/gmt.conf
gmt [DEBUG]: Look for file ~/.gmt/cache/gmt.conf
gmt [DEBUG]: Could not find file gmt.conf
gmt [DEBUG]: Exit:  gmt_reload_settings
gmt [DEBUG]: Enter: gmtlib_plot_C_format
gmt [DEBUG]: Exit:  gmtlib_plot_C_format
gmt [DEBUG]: Enter: gmtinit_get_history
gmt [DEBUG]: Initialize FFTW with 16 threads.
gmt [DEBUG]: GMT_Create_Session initialized GMT structure
gmt [DEBUG]: Loading core GMT shared library: libgmt.so
gmt [DEBUG]: Shared Library # 0 (core). Path = libgmt.so
gmt [DEBUG]: Loading GMT plugins from: ~/mambaforge/envs/pygmt/lib/gmt/plugins
gmt [DEBUG]: Shared Library # 1 (supplements). Path = ~/mambaforge/envs/pygmt/lib/gmt/plugins/supplements.so
gmt [DEBUG]: Local file ~/.gmt/server/gmt_data_server.txt found
gmt [DEBUG]: File ~/.gmt/server/gmt_data_server.txt less than 24 hours old, refresh is premature.
gmt [DEBUG]: Load contents from ~/.gmt/server/gmt_data_server.txt
gmt [DEBUG]: Local file ~/.gmt/server/gmt_hash_server.txt found
gmt [DEBUG]: File ~/.gmt/server/gmt_hash_server.txt less than 24 hours old, refresh is premature.
gmt [WARNING]: Remote dataset given to a data processing module but no registration was specified - default to gridline registration (if available)
gmt [DEBUG]: Input remote grid modified to have registration: @earth_day_01d_p
gmt [DEBUG]: Map distance calculation will be using great circle approximation with authalic auxiliary latitudes and authalic (R_2) radius = 6371007.1809 m, in meter.
gmt [DEBUG]: Revised options: @earth_day_01d_p -Vd
gmtwhich [DEBUG]: gmtapi_init_export: Passed family = Data Table and geometry = Non-Geographical
gmtwhich [DEBUG]: Object ID 0 : Registered Data Table Stream 7fd32ddf2780 as an Output resource with geometry Non-Geographical [n_objects = 1]
gmtwhich [DEBUG]: gmtapi_init_export: Added stdout to registered destinations
gmtwhich [DEBUG]: GMT_Init_IO: Returned first Output object ID = 0
gmtwhich [DEBUG]: GMT_Begin_IO: Initialize record-by-record access for Output
gmtwhich [DEBUG]: gmtapi_next_io_source: Selected object 0
gmtwhich [INFORMATION]: Writing Data Table to Standard Output stream
gmtwhich [DEBUG]: GMT_Begin_IO: Output resource access is now enabled [record-by-record]
gmtwhich [DEBUG]: Look for file earth_day_01d_p.tif in ~/.gmt
gmtwhich [DEBUG]: Look for file earth_day_01d_p.tif in ~/.gmt/cache
gmtwhich [DEBUG]: Look for file earth_day_01d_p.tif in ~/.gmt/server
gmtwhich [ERROR]: File earth_day_01d_p.tif not found!
gmtwhich [DEBUG]: GMT_End_IO: Output resource access is now disabled
gmtwhich [DEBUG]: gmtlib_unregister_io: Unregistering object no 0 [n_objects = 0]
gmt [DEBUG]: Entering GMT_Destroy_Session

But with GMT 6.5.0, gmt which @earth_day_01d_p works fine.

Also, I noticed that the internal data structure of the Blue/Black Marble images have changed from a 1-band ColorInterp/Color Table structure to a 3-band structure:

diff --git a/gdalinfo_earth_day b/gdalinfo_earth_day
index 75e7bbc91..02d818804 100644
--- a/gdalinfo_earth_day
+++ b/gdalinfo_earth_day
@@ -1,5 +1,5 @@
 Driver: GTiff/GeoTIFF
-Files: earth_day_01d_p.tif
+Files: ~/.gmt/cache/earth_day_01d_p.tif
 Size is 360, 180
 Coordinate System is:
 GEOGCRS["unknown",
@@ -26,269 +26,40 @@ Metadata:
   AREA_OR_POINT=Area
 Image Structure Metadata:
   COMPRESSION=DEFLATE
-  INTERLEAVE=BAND
-  PREDICTOR=2
+  INTERLEAVE=PIXEL
 Corner Coordinates:
 Upper Left  (-180.0000000,  90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"N)
 Lower Left  (-180.0000000, -90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"S)
 Upper Right ( 180.0000000,  90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"N)
 Lower Right ( 180.0000000, -90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"S)
 Center      (   0.0000000,   0.0000000) (  0d 0' 0.01"E,  0d 0' 0.01"N)
-Band 1 Block=360x22 Type=Byte, ColorInterp=Palette
-  Color Table (RGB with 256 entries)
-    0: 248,224,168,255
-    1: 240,224,152,255
-    2: 104,104,96,255
-    3: 16,24,56,255
-    4: 24,28,24,255
-    5: 104,80,48,255
-    6: 220,224,228,255
-    7: 48,56,32,255
-    8: 184,192,176,255
-    9: 104,104,40,255
-   10: 16,48,8,255
-   11: 96,96,56,255
-   12: 188,184,184,255
-   13: 128,112,60,255
-   14: 48,72,24,255
-   15: 8,16,48,255
-   16: 216,200,192,255
-   17: 80,80,28,255
-   18: 56,48,44,255
-   19: 72,88,24,255
-   20: 204,168,140,255
-   21: 128,136,136,255
-   22: 72,68,40,255
-   23: 144,120,68,255
-   24: 144,144,136,255
-   25: 88,80,48,255
-   26: 168,160,152,255
-   27: 88,104,40,255
-   28: 8,16,32,255
-   29: 168,116,84,255
-   30: 120,100,48,255
-   31: 200,192,184,255
-   32: 176,168,160,255
-   33: 136,112,80,255
-   34: 204,204,200,255
-   35: 184,152,124,255
-   36: 96,76,48,255
-   37: 88,72,48,255
-   38: 184,196,192,255
-   39: 28,56,48,255
-   40: 92,100,124,255
-   41: 184,176,176,255
-   42: 160,160,152,255
-   43: 132,144,144,255
-   44: 208,192,184,255
-   45: 128,100,64,255
-   46: 72,100,24,255
-   47: 36,48,16,255
-   48: 140,140,80,255
-   49: 152,144,136,255
-   50: 56,64,16,255
-   51: 168,168,168,255
-   52: 224,200,144,255
-   53: 176,176,168,255
-   54: 80,88,28,255
-   55: 152,144,144,255
-   56: 156,148,84,255
-   57: 160,160,160,255
-   58: 16,16,48,255
-   59: 168,176,168,255
-   60: 216,200,176,255
-   61: 16,8,48,255
-   62: 72,80,24,255
-   63: 132,140,124,255
-   64: 88,80,40,255
-   65: 152,152,176,255
-   66: 16,16,56,255
-   67: 20,48,24,255
-   68: 200,200,192,255
-   69: 56,48,72,255
-   70: 40,56,24,255
-   71: 120,136,56,255
-   72: 152,160,152,255
-   73: 8,8,24,255
-   74: 192,184,176,255
-   75: 248,244,172,255
-   76: 116,124,136,255
-   77: 88,68,40,255
-   78: 64,72,32,255
-   79: 172,180,176,255
-   80: 176,152,132,255
-   81: 24,40,24,255
-   82: 60,100,24,255
-   83: 168,148,128,255
-   84: 176,132,88,255
-   85: 200,192,192,255
-   86: 76,76,48,255
-   87: 152,116,80,255
-   88: 8,16,60,255
-   89: 152,152,144,255
-   90: 240,216,152,255
-   91: 52,80,32,255
-   92: 20,24,40,255
-   93: 132,144,156,255
-   94: 88,88,48,255
-   95: 20,32,8,255
-   96: 84,68,28,255
-   97: 172,164,180,255
-   98: 152,152,152,255
-   99: 120,112,56,255
-  100: 88,88,40,255
-  101: 56,80,8,255
-  102: 80,80,40,255
-  103: 176,176,160,255
-  104: 48,56,8,255
-  105: 104,84,68,255
-  106: 48,120,76,255
-  107: 184,176,168,255
-  108: 108,96,56,255
-  109: 232,224,212,255
-  110: 16,16,40,255
-  111: 96,88,48,255
-  112: 144,152,144,255
-  113: 104,112,68,255
-  114: 220,180,148,255
-  115: 28,32,60,255
-  116: 56,72,8,255
-  117: 104,96,40,255
-  118: 56,72,24,255
-  119: 24,64,24,255
-  120: 148,164,160,255
-  121: 180,176,140,255
-  122: 216,208,196,255
-  123: 156,128,92,255
-  124: 72,84,40,255
-  125: 56,80,16,255
-  126: 72,72,24,255
-  127: 48,44,12,255
-  128: 124,108,96,255
-  129: 64,80,24,255
-  130: 132,136,148,255
-  131: 200,200,184,255
-  132: 120,100,68,255
-  133: 16,28,24,255
-  134: 144,144,144,255
-  135: 116,120,124,255
-  136: 248,224,184,255
-  137: 184,184,168,255
-  138: 72,80,32,255
-  139: 152,148,160,255
-  140: 56,76,44,255
-  141: 232,208,144,255
-  142: 176,168,168,255
-  143: 104,96,88,255
-  144: 60,56,32,255
-  145: 168,152,92,255
-  146: 184,184,176,255
-  147: 152,144,124,255
-  148: 232,180,120,255
-  149: 208,204,192,255
-  150: 152,152,128,255
-  151: 188,180,108,255
-  152: 48,56,16,255
-  153: 28,60,40,255
-  154: 56,64,24,255
-  155: 160,152,144,255
-  156: 48,72,8,255
-  157: 168,160,160,255
-  158: 228,228,184,255
-  159: 48,64,8,255
-  160: 160,152,152,255
-  161: 148,112,64,255
-  162: 32,48,8,255
-  163: 80,68,40,255
-  164: 132,140,100,255
-  165: 84,104,28,255
-  166: 36,48,28,255
-  167: 188,192,184,255
-  168: 168,168,160,255
-  169: 240,248,248,255
-  170: 60,52,24,255
-  171: 248,248,220,255
-  172: 36,36,8,255
-  173: 144,152,152,255
-  174: 32,64,8,255
-  175: 96,80,80,255
-  176: 72,140,148,255
-  177: 48,72,16,255
-  178: 200,184,176,255
-  179: 148,144,152,255
-  180: 136,140,136,255
-  181: 244,200,144,255
-  182: 180,200,212,255
-  183: 160,168,160,255
-  184: 40,64,8,255
-  185: 56,72,16,255
-  186: 48,56,24,255
-  187: 64,72,24,255
-  188: 176,184,168,255
-  189: 48,64,16,255
-  190: 104,116,56,255
-  191: 28,24,56,255
-  192: 132,124,136,255
-  193: 72,72,32,255
-  194: 80,88,40,255
-  195: 56,88,12,255
-  196: 236,204,184,255
-  197: 72,64,28,255
-  198: 64,84,32,255
-  199: 24,48,8,255
-  200: 88,80,84,255
-  201: 240,232,236,255
-  202: 148,100,60,255
-  203: 148,116,100,255
-  204: 160,144,140,255
-  205: 208,132,88,255
-  206: 36,32,24,255
-  207: 148,132,132,255
-  208: 52,84,24,255
-  209: 196,180,168,255
-  210: 204,188,140,255
-  211: 24,16,52,255
-  212: 104,120,40,255
-  213: 240,184,144,255
-  214: 72,88,32,255
-  215: 104,84,40,255
-  216: 60,52,16,255
-  217: 208,200,184,255
-  218: 32,56,8,255
-  219: 24,56,8,255
-  220: 160,168,152,255
-  221: 116,136,164,255
-  222: 204,196,172,255
-  223: 48,84,12,255
-  224: 88,128,80,255
-  225: 64,64,24,255
-  226: 248,244,192,255
-  227: 172,180,192,255
-  228: 144,144,60,255
-  229: 40,56,8,255
-  230: 76,72,88,255
-  231: 40,72,24,255
-  232: 196,180,160,255
-  233: 152,168,132,255
-  234: 64,88,24,255
-  235: 64,64,16,255
-  236: 40,48,8,255
-  237: 220,152,124,255
-  238: 92,84,28,255
-  239: 204,208,220,255
-  240: 104,128,96,255
-  241: 132,168,164,255
-  242: 40,76,8,255
-  243: 188,184,196,255
-  244: 56,76,80,255
-  245: 156,148,104,255
-  246: 244,240,156,255
-  247: 68,80,16,255
-  248: 28,64,80,255
-  249: 152,176,188,255
-  250: 132,84,56,255
-  251: 8,8,48,255
-  252: 8,8,56,255
-  253: 8,40,24,255
-  254: 8,72,96,255
-  255: 248,248,248,255
+Band 1 Block=360x7 Type=Byte, ColorInterp=Undefined
+  Min=10.000 Max=255.000
+  Minimum=10.000, Maximum=255.000, Mean=49.508, StdDev=65.299
+  NoData Value=0
+  Metadata:
+    STATISTICS_MAXIMUM=255
+    STATISTICS_MEAN=49.507654320988
+    STATISTICS_MINIMUM=10
+    STATISTICS_STDDEV=65.299010755209
+    STATISTICS_VALID_PERCENT=100
+Band 2 Block=360x7 Type=Byte, ColorInterp=Undefined
+  Min=10.000 Max=255.000
+  Minimum=10.000, Maximum=255.000, Mean=49.911, StdDev=62.129
+  NoData Value=0
+  Metadata:
+    STATISTICS_MAXIMUM=255
+    STATISTICS_MEAN=49.910941358025
+    STATISTICS_MINIMUM=10
+    STATISTICS_STDDEV=62.129078727248
+    STATISTICS_VALID_PERCENT=100
+Band 3 Block=360x7 Type=Byte, ColorInterp=Undefined
+  Min=10.000 Max=255.000
+  Minimum=10.000, Maximum=255.000, Mean=65.605, StdDev=44.820
+  NoData Value=0
+  Metadata:
+    STATISTICS_MAXIMUM=255
+    STATISTICS_MEAN=65.605154320988
+    STATISTICS_MINIMUM=10
+    STATISTICS_STDDEV=44.819618148349
+    STATISTICS_VALID_PERCENT=100

Looking at https://oceania.generic-mapping-tools.org/gmt_data_server.txt, the modification date was 2023-09-29, and seems to be applied from https://github.com/GenericMappingTools/gmtserver-admin/commit/5bcc3d6190cf88f9ccdbb9ed3ca80451e94e61b2.

Just checking if these changes are expected. Personally, I prefer the new 3-band structure since it will simplify things in PyGMT (https://github.com/GenericMappingTools/pygmt/issues/1442#issuecomment-1341079313), but wanted to know if breaking the downloads for GMT 6.4.0 are ok.

PaulWessel commented 6 months ago

So here is what happened last September. We realised our spherical filtering was incorrect with a too-short filter radius. All data sets were reprocessed with the updated gmtserver-admin scripts. For the NASA images it is true they were indexed images (at least the full resolution 30s still is I think since no filtering done there). From memory (this is years ago) I think I just used photoshop to downsample the full file to get the lower resolutions using some semi-reasonable filter width. However, with the September work we had to do actual filtering and one cannot filter indices, hence we used grdmix to pull out red, green, and blue grids, filtered these, then grdmix put these filtered r,g,b back into one image, which now was RGB format. GMT has no ability to convert that to an indexed image, but probably GDAL does. We never did but we could discuss if that is important. Biggest problem with the images is we have not tiled them like we do for the grids. Hoping @joa-quim can provide some suggestions here. I think GMT 6.4 and earlier can plot RGB images as well as indexed images (since our GDAL bridge returns index images as RGB anyway).

seisman commented 6 months ago

Are the Blue Marble @earth_day_01d_p and Black Marble @earth_night_01d_p datasets only downloadable with GMT 6.5.0 (and not GMT 6.4.0 or older) now?

E.g. doing gmt which @earth_day_01d_p -Vd with GMT 6.4.0 gives:

But with GMT 6.5.0, gmt which @earth_day_01d_p works fine.

Just to clarify that to access these datasets, you need gmt which -Ga @earth_day_01d_p instead. Without -Ga, which won't download the dataset.

joa-quim commented 6 months ago

Paul, you used a color quantization option to go from RGB to indexed from a GDAL module that I forgot the name

PaulWessel commented 6 months ago

Thinking it must be _gdaltranslate but nothing obvious stand out. One argument for quantising to 256 unique colors via an index table is that our day/night images are not tiled so to use the earth_day_30s.tif means it is a singe file of 622 Mb. In absence of tiling we could probably get that down to 1/3 with indexing. But, what is the amount of work needed to enable _gmtgdalwrite.c learn a new trick to do this for us? I guess the quantising is tricky. How would we do it? Compute difference between all pixels in RGB space and pick the 256 points the have the highest sum of distances between them? So 255 244 239 will be selected as the final RGB index entry if other nearby colors are relatively near? Perhaps better to find what the GDAL option you refer to is.. Reading in index images is way more useful (and we do) than writing them out, which is encroaching on image processing stuff.

PaulWessel commented 6 months ago

Pretty sure I used PhotoShop who as quantification/index option. However, only lets me save as PNG, not TIFF, and I am pretty sure we expect the da/night images to be TIFFs from GMT's point of view. Perhaps need to convert the PNG index image to a TIFF index image with GDAL.

PaulWessel commented 6 months ago

Experiment: The 01m day tiff (RGB) is 171M but a photoshop -> index PNG is 53M and Preview can save this to index TIF at 47 Mb. I can plot the whole thing with

gmt grdimage earth_day_01m_indexed.tiff -Dr -Rd -JQ20c -B -pdf map

since the -180/180/-90/90 got lost probably in the commercial tools. Does reduce transmission size/time by 3 though. @joa-quim, if you can tell me how to add back in the reference (-Rd) to that indexed tiff file than we can revert the day/night filtered RGB images to indexed TIFF for once via PhotoShop->Preview->add-CRS-> indexes final file 1/3 the size and works in GMT.

PaulWessel commented 6 months ago

Looks like a fun research problem to apply the k-means statistics to determine the 256 best colors but how often is this needed. I think it is tolerable to do a bi manual labor on those images to shrink and speed them up via indices, no? What do @weiji14 say?

joa-quim commented 6 months ago

I only find this now (a python script) https://manpages.opensuse.org/Tumbleweed/gdal/rgb2pct.1.en.html

I have this implemented in Mirone. Probably from an old MEX

But I don't think it worth's spending time on this.

weiji14 commented 6 months ago

Just to clarify that to access these datasets, you need gmt which -Ga @earth_day_01d_p instead. Without -Ga, which won't download the dataset.

Ah I see, thanks @seisman for clarifying! Somehow it worked without -Ga on GMT 6.5.0, but I just tested with -Ga on GMT 6.4.0 and it worked.

Personally I'd prefer the non-quantized version of the GeoTIFF right now (even though they are bigger) instead of the quantized/indexed version. Nowadays with GeoTIFFs, GDAL should be possible to read spatial subsets without downloading the whole file, as long as the data is arranged as square tiles internally. Could we convert the GeoTIFFs to Cloud-Optimized GeoTIFFs, e.g. using GDAL's COG driver?