Function to load NASA Blue and Black marble mosaic images

[weiji14]

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).

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:

pygmt/pygmt/datasets/earth_relief.py

Lines 92 to 140 in ebbb0cc

	# earth relief data stored as single grids for low resolutions
	non_tiled_resolutions = ["01d", "30m", "20m", "15m", "10m", "06m"]
	# earth relief data stored as tiles for high resolutions
	tiled_resolutions = ["05m", "04m", "03m", "02m", "01m", "30s", "15s", "03s", "01s"]
	# resolutions of original land-only SRTM tiles from NASA
	land_only_srtm_resolutions = ["03s", "01s"]
	if registration in ("pixel", "gridline", None):
	# If None, let GMT decide on Pixel/Gridline type
	reg = f"_{registration[0]}" if registration else ""
	else:
	raise GMTInvalidInput(
	f"Invalid grid registration: '{registration}', should be either "
	"'pixel', 'gridline' or None. Default is None, where a "
	"pixel-registered grid is returned unless only the "
	"gridline-registered grid is available."
	)
	if resolution not in non_tiled_resolutions + tiled_resolutions:
	raise GMTInvalidInput(f"Invalid Earth relief resolution '{resolution}'.")
	# Check combination of resolution and registration.
	if (resolution == "15s" and registration == "gridline") or (
	resolution in ("03s", "01s") and registration == "pixel"
	):
	raise GMTInvalidInput(
	f"{registration}-registered Earth relief data for "
	f"resolution '{resolution}' is not supported."
	)
	# Choose earth relief data prefix
	earth_relief_prefix = "earth_relief_"
	if use_srtm and resolution in land_only_srtm_resolutions:
	earth_relief_prefix = "srtm_relief_"
	# different ways to load tiled and non-tiled earth relief data
	# Known issue: tiled grids don't support slice operation
	# See https://github.com/GenericMappingTools/pygmt/issues/524
	if region is None:
	if resolution not in non_tiled_resolutions:
	raise GMTInvalidInput(
	f"'region' is required for Earth relief resolution '{resolution}'."
	)
	fname = which(f"@earth_relief_{resolution}{reg}", download="a")
	with xr.open_dataarray(fname, engine="netcdf4") as dataarray:
	grid = dataarray.load()
	_ = grid.gmt # load GMTDataArray accessor information
	else:
	grid = grdcut(f"@{earth_relief_prefix}{resolution}{reg}", region=region)

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]

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.

[willschlitzer]

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]

@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 corteva/rioxarray#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.

[weiji14]

FYI, with GMT_IMAGE implemented in #3338 and #3128, we no longer need rioxarray to open 3-band GeoTIFF images 🎉. I.e., the Blue Marble GeoTIFFs can be read using GMT directly with some special logic that will be implemented in the load_blue_marble function that I'll try to finalize in #2235.