Refactor _load_remote_dataset to write CRS information for raster images

[seisman]

Write CRS information in _load_remote_dataset, so that we don't have to repeat it in load_blue_marble/load_black_marblue.

I have very little knowledge about CRS, datums or geodesy stuff, so here are my questions:

1. Does it makes sense to also write CRS for grids, rather than just images?
2. Currently, we define attributes horizontal_datum and vertical_datum. These attributes were added since the start of the project (Add metadata to Earth relief grids #248). As far as I can see, these attributes are not standard attributes in CF-convention (https://cfconventions.org/Data/cf-conventions/cf-conventions-1.8/cf-conventions.html), and they are part of the CRS. So, perhaps we should remove these two attributes?

[seisman]

1. Does it makes sense to also write CRS for grids, rather than just images?
2. Currently, we define attributes horizontal_datum and vertical_datum. These attributes were added since the start of the project (Add metadata to Earth relief grids #248). As far as I can see, these attributes are not standard attributes in CF-convention (https://cfconventions.org/Data/cf-conventions/cf-conventions-1.8/cf-conventions.html), and they are part of the CRS. So, perhaps we should remove these two attributes?

ChatGPT said yes to both questions. @weiji14 I guess you're more familiar with CRS, so what do you think?

[weiji14]

1. Does it makes sense to also write CRS for grids, rather than just images?

Grids can have a CRS attached yes. But how it should be set in xarray.DataArray grids is still a topic of debate (pydata/xarray#2288). TLDR is that rioxarray's .rio.crs accessor is somewhat the de-facto default today, but there's discussions on whether to have a more 'global' .crs accessor that can be used by other libraries (including PyGMT) without depending on rioxarray that has a heavy dependency on rasterio/GDAL.

So I would try to avoid adding CRS to grids in PyGMT for now (until the xarray community agrees on a standard) to avoid having to rewrite it later. The reason we have added CRS to 3-band images (for tile_map) is so we can do reprojection from EPSG:3857 (Spherical Mercator) to OGC:CRS84 (lon/lat).

2. Currently, we define attributes horizontal_datum and vertical_datum. These attributes were added since the start of the project (Add metadata to Earth relief grids #248).

Probably fine to keep those two attributes, even if we decide to set the CRS (more metadata is better than less metadata). A CRS doesn't necessarily specify a horizontal and/or vertical datum if I'm not mistaken (looking at some PROJ examples). More specifically, the CRS (assuming PROJ.4 syntax) should have the +datum and +ellps attributes to map to CF grid_mapping (see https://github.com/cf-convention/cf-conventions/wiki/Mapping-from-CF-Grid-Mapping-Attributes-to-CRS-WKT-Elements#table-2---future-cf-17-cf-grid-mapping-attributes)

As far as I can see, these attributes are not standard attributes in CF-convention (https://cfconventions.org/Data/cf-conventions/cf-conventions-1.8/cf-conventions.html), and they are part of the CRS. So, perhaps we should remove these two attributes?

CF-compliance is tricky... There is a horizontal_datum_name attribute under grid_mapping (see https://cfconventions.org/Data/cf-conventions/cf-conventions-1.10/cf-conventions.html#grid-mappings-and-projections and https://cf-xarray.readthedocs.io/en/stable/grid_mappings.html). Closest thing to vertical_datum would be reference_ellipsoid_name I think. The CF-convention docs does have this note:

reference_ellipsoid_name, prime_meridian_name, horizontal_datum_name and geographic_crs_name must be all defined if any one is defined, and if projected_crs_name is defined then geographic_crs_name must be also.

So we would need to set all 4 attributes if following CF-convention 🙂

[seisman]

OK, we can wait to see when the community can agree about the ".crs" accessor.

Currently, this PR adds the crs parameter to the internal GMTRemoteDataset class, and it defaults to None. For image datasets (i.e., earth_day and earth_night), it's set to OGC:CRS84 and the rio.crs metadata is written.

OK with these changes?

[weiji14]

Yep, ok with the changes keeping things DRY.