GenericMappingTools / pygmt

A Python interface for the Generic Mapping Tools.
https://www.pygmt.org
BSD 3-Clause "New" or "Revised" License
745 stars 216 forks source link

Function to load NASA Blue and Black marble mosaic images #1442

Open weiji14 opened 3 years ago

weiji14 commented 3 years ago

Description of the desired feature

Provide access to cloud-free colour images of the Earth! Specifically, the @earth_day_ and @earth_night_ GMT remote data files (see https://docs.generic-mapping-tools.org/6.2/datasets/remote-data.html#global-earth-day-night-images).

image

Context is that I'm using it at #1437, and it would be good to have some sample RGB images. We might as well have a proper function similar to pygmt.datasets.load_earth_relief, but for the Blue and Black Marble images instead.

Note that there is a lot of logic from the load_earth_relief function that could be reused:

https://github.com/GenericMappingTools/pygmt/blob/ebbb0cc0045454f05db5abdfb3bcb5e1f18e96a7/pygmt/datasets/earth_relief.py#L92-L140

These are some ways I'm considering:

Option 1 - rename the earth_relief.py file to earth_data.py first, and then separate the re-usable chunk of code into a function called _load_earth_data or something and have load_earth_relief/load_earth_day/load_earth_night call that. This is similar to what we did for blockmean and blockmedian at #1092. Less breakages for users already using load_earth_relief, but maybe more work for the devs.

Option 2 - have a single function called load_earth_data which accepts a parameter called 'type' that accepts arguments like 'relief'/'day'/'night'. We could then pretty much re-use the existing load_earth_relief function though it will need to be renamed to load_earth_data I guess. This would involve a deprecation notice for users using load_earth_relief, but less work for the devs in terms of not needing to write a new function?

Option 3 - any other ideas?

Are you willing to help implement and maintain this feature? Best to discuss first.

seisman commented 3 years ago

There are many differences between earth_relief and earth_day/earth_night datasets:

I think these differences mean that very few code in the earth_relief.py file can be reused.

willschlitzer commented 1 year ago

There are many differences between earth_relief and earth_day/earth_night datasets:

* earth_day/earth_night datasets are in geotiff format, not netcdf.

  * What's the return type for the `load_earth_day` function?
  * Does xarray support reading geotiff?
  * If it does support, does the returned DataArray can be passed to the `grdimage` module?

* earth_day/earth_night only provide pixel-registered version, so the `registration` parameter is not needed

* each resolution of earth_day/earth_night is stored in a single file, but earth_relief data higher than 06m resolution are stored in tiles.

I think these differences mean that very few code in the earth_relief.py file can be reused.

@seisman and @weiji14 I think I have made a fix for problems of no registration or region passed to _load_remote_dataset. Would I be able to change the engine argument for load_dataarray for something that can accept a TIFF file?

weiji14 commented 1 year ago

@seisman and @weiji14 I think I have made a fix for problems of no registration or region passed to _load_remote_dataset. Would I be able to change the engine argument for load_dataarray for something that can accept a TIFF file?

Kinda. The tricky part is that Blue/Black Marble are RGB images with 3 bands, whereas the other datasets you've been working on are 1 band datasets. So you would need to use rioxarray.open_rasterio instead.

Now one complication is that these GeoTIFF files store the RGB data in a colormap attribute, and there's no standardized way on how to read these into an xarray.DataArray via rioxarray yet. I did try for a while at https://github.com/corteva/rioxarray/pull/414, but haven't got around to finishing that Pull Request. If you've got time, it's worth testing a few things and see what works and what doesn't.