Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Open
weiji14 opened this issue Jan 7, 2024 · 9 comments
Open

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

weiji14 opened this issue Jan 7, 2024 · 9 comments

Comments

@weiji14
Copy link
Member

weiji14 commented Jan 7, 2024

Just checking something that I noticed at GenericMappingTools/pygmt#2963 (comment). 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 5bcc3d6.

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

@PaulWessel
Copy link
Member

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
Copy link
Member

seisman commented Jan 8, 2024

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
Copy link
Member

joa-quim commented Jan 8, 2024

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

@PaulWessel
Copy link
Member

Thinking it must be gdal_translate 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 gmt_gdalwrite.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
Copy link
Member

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
Copy link
Member

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
Copy link
Member

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
Copy link
Member

joa-quim commented Jan 8, 2024

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
Copy link
Member Author

weiji14 commented Jan 8, 2024

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?

weiji14 added a commit to GenericMappingTools/pygmt that referenced this issue Jan 8, 2024
Downloading @earth_day_01d_p using GMT 6.4.0 actually works. The new earth_day_01d_p.tif file processed on 2023-09-29 doesn't have the quantized colormap, so no need for the extra handling. See GenericMappingTools/gmtserver-admin#257 for more info.
weiji14 added a commit to GenericMappingTools/pygmt that referenced this issue Jan 11, 2024
…2963)

Rioxarray can now read the 3-bands from `@earth_day_01d_p` directly
without needing to parse the colorinterp information from the GeoTIFF file.

* Update baseline images for test_tilemap_* with GMT 6.5.0 and gs 10.02.1

* Trigger re-caching of earth_day_01d_p file

Modify pygmt/helpers/caching.py file slightly to create new cache.

* Remove elif block for manually creating RGB image

Downloading @earth_day_01d_p using GMT 6.4.0 actually works.
The new earth_day_01d_p.tif file processed on 2023-09-29 doesn't have
the quantized colormap, so no need for the extra handling. See
GenericMappingTools/gmtserver-admin#257 for more info.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

4 participants