GMT_IMAGE: Implement the GMT_IMAGE.to_dataarray method for 3-band images #3128

seisman · 2024-03-20T03:01:58Z

This PR implements the to_dataarray method for GMT_IMAGE.

Similar to #2398 but for images.

~~Need to wait for #3446.~~

Minimal xarray.DataArray output with data and coordinates, no metadata yet.

pygmt/clib/session.py

seisman · 2024-04-14T04:46:54Z

xref GenericMappingTools/gmt#3299

pygmt/datatypes/image.py

Extra metadata from the _GMT_GRID_HEADER struct.

weiji14 · 2024-06-18T10:13:04Z

pygmt/datatypes/image.py

+            name=header.name,
+            attrs=header.data_attrs,


The image name is currently hardcoded to z, is that ok for an RGB image?

pygmt/pygmt/datatypes/header.py

Lines 193 to 198 in ac44706

@property

def name(self) -> str:

"""

Name of the grid.

"""

return "z"

The attrs fields might need some work. I'm getting 'actual_range': array([ 1.79769313e+308, -1.79769313e+308])} when loading the @earth_day_01d image.

I agree they make no sense, but they're consistent with the behavior in GMT.

gmt grdinfo @earth_day_01d /Users/seisman/.gmt/server/earth/earth_day/earth_day_01d_p.tif: Title: Grid imported via GDAL /Users/seisman/.gmt/server/earth/earth_day/earth_day_01d_p.tif: Command: /Users/seisman/.gmt/server/earth/earth_day/earth_day_01d_p.tif: Remark: /Users/seisman/.gmt/server/earth/earth_day/earth_day_01d_p.tif: Pixel node registration used [Geographic grid] /Users/seisman/.gmt/server/earth/earth_day/earth_day_01d_p.tif: Grid file format: gd = Import/export through GDAL /Users/seisman/.gmt/server/earth/earth_day/earth_day_01d_p.tif: x_min: -180 x_max: 180 x_inc: 1 name: x n_columns: 360 /Users/seisman/.gmt/server/earth/earth_day/earth_day_01d_p.tif: y_min: -90 y_max: 90 y_inc: 1 name: y n_rows: 180 /Users/seisman/.gmt/server/earth/earth_day/earth_day_01d_p.tif: v_min: 1.79769313486e+308 v_max: -1.79769313486e+308 name: z /Users/seisman/.gmt/server/earth/earth_day/earth_day_01d_p.tif: scale_factor: 1 add_offset: 0 /Users/seisman/.gmt/server/earth/earth_day/earth_day_01d_p.tif: Default CPT: +proj=longlat +R=6378137 +no_defs

The GMT's image support was likely added by Joaquim so that you may ping him for more information.

Reorder the dimensions to follow Channel, Height, Width (CHW) convention. Also added doctest checking output DataArray object and the image's x and y coordinates.

Get the registration and gtype info from the grid header and apply it to the GMT accessor attributes.

Trying to match some of the doctests in _GMT_GRID.

Remove hardcoded attribute that was only meant for NetCDF files, so that GeoTIFF files won't have it.

pygmt/datatypes/header.py

pygmt/clib/session.py

pygmt/datatypes/image.py

seisman · 2024-09-28T14:25:22Z

I think this PR is almost ready for review. Here are the known limitations/issues:

Currently, only 3-band images are supported. [We can work on 1-band and 4-band images in a separate PR to make it easier to review]
As mentioned in GMT_IMAGE: Implement the GMT_IMAGE.to_dataarray method for 3-band images #3128 (comment), the name is hardcoded to z and the actual_range doesn't work. [For the hardcoded name z, I'm not sure if we can do better; for actual_range, I think we can remove it because (1) actual_range makes no sense for 3-band RGB images; (2) actual_range is the standard attribution defined by CF-1.7 convention, but it's only for netCDF files, not for images; (3) the xarray.DataArray returned by rioxarray.open_rasterio doesn't have the actual_range attribution].
For the y coordinates, both GMT API and rioxarray return a descending y array, but in _GMT_GRID.to_dataarray() method, we flip the y coordinate to make it ascending. Do we want to do the same thing for images?

seisman · 2024-09-28T15:39:57Z

pygmt/datatypes/image.py

+            title:
+            history:
+            description:


Maybe don't add these attributes if they're empty?

These come from the header no? For GeoTIFF files, they don't make sense since GMT is not parsing the TIFF tags, but the title/history/description might show up for 3-band NetCDF files (if those can be passed into this GMT_IMAGE struct).

Right, title/history/description are recommended attributes of CF-1.7. So, for netCDF files, I guess we should always show them even if they're empty.

If those are only needed for CF-1.7, then maybe we can indent the 3 attrs["..."] lines here to be inside the if self.type in {GridFormat.NX} block?

pygmt/pygmt/datatypes/header.py

Lines 207 to 217 in cf48764

if self.type in {

GridFormat.NB,

GridFormat.NS,

GridFormat.NI,

GridFormat.NF,

GridFormat.ND,

}: # Only set the 'Conventions' attribute for netCDF.

attrs["Conventions"] = "CF-1.7"

attrs["title"] = self.title.decode()

attrs["history"] = self.command.decode()

attrs["description"] = self.remark.decode()

(2) actual_range is the standard attribution defined by CF-1.7 convention, but it's only for netCDF files, not for images;

Same here, we can put attrs["actual_range"] under the if-block above:

pygmt/pygmt/datatypes/header.py

Line 223 in cf48764

attrs["actual_range"] = np.array([self.z_min, self.z_max])

Done at 62f0ce0.

weiji14

Currently, only 3-band images are supported. [We can work on 1-band and 4-band images in a separate PR to make it easier to review]

Yep, just do 3-band images here for now.

For the y coordinates, both GMT API and rioxarray return a descending y array, but in _GMT_GRID.to_dataarray() method, we flip the y coordinate to make it ascending. Do we want to do the same thing for images?

I'd prefer descending y to match that the output of rioxarray.open_rasterio. In #2398 (comment), we kept the ascending y for GMT_GRID for pragmatic reasons, namely to match the output of xr.load_datarray so that the unit tests won't break. But I think the convention for most xarray.DataArray objects is to have descending y and ascending x.

weiji14 · 2024-09-29T03:47:15Z

pygmt/datatypes/image.py

+        >>> da.gmt.registration, da.gmt.gtype
+        (1, 0)


The gtype for earth_day_01d should be Geographic (1), not Cartesian (0).

Suggested change

>>> da.gmt.registration, da.gmt.gtype

(1, 0)

>>> da.gmt.registration, da.gmt.gtype

(1, 1)

according to:

$ gmt grdinfo @earth_day_01d ~/.gmt/server/earth/earth_day/earth_day_01d_p.tif: Title: Grid imported via GDAL ~/.gmt/server/earth/earth_day/earth_day_01d_p.tif: Command: ~/.gmt/server/earth/earth_day/earth_day_01d_p.tif: Remark: ~/.gmt/server/earth/earth_day/earth_day_01d_p.tif: Pixel node registration used [Geographic grid] ~/.gmt/server/earth/earth_day/earth_day_01d_p.tif: Grid file format: gd = Import/export through GDAL ~/.gmt/server/earth/earth_day/earth_day_01d_p.tif: x_min: -180 x_max: 180 x_inc: 1 name: x n_columns: 360 ~/.gmt/server/earth/earth_day/earth_day_01d_p.tif: y_min: -90 y_max: 90 y_inc: 1 name: y n_rows: 180 ~/.gmt/server/earth/earth_day/earth_day_01d_p.tif: v_min: 1.79769313486e+308 v_max: -1.79769313486e+308 name: z ~/.gmt/server/earth/earth_day/earth_day_01d_p.tif: scale_factor: 1 add_offset: 0 ~/.gmt/server/earth/earth_day/earth_day_01d_p.tif: Default CPT:

We either need to modify the logic here:

pygmt/pygmt/datatypes/header.py

Line 254 in cf48764

gtype = 1 if dims[0] == "lat" and dims[1] == "lon" else 0

Or I can manually override the gtype in the load_blue_marble function which I'm trying to finish up at #2235.

https://github.com/GenericMappingTools/gmt/blob/7809736ba32d87a4a96b15444419eb176c6a35f3/src/gmt_grdio.c#L1360

This is the function (again, not a public API function) that GMT uses to determine the grid type. The logic here is not that complicated to duplicate in Python.

How easy would it be to make that C function public actually?

How easy would it be to make that C function public actually?

It doesn't seem too difficult. But we still need to implement our own functions if we want to be compatible with GMT 6.4-6.5.

Yeah, it'll be some time before GMT 6.6 comes out. I just want the gtype determination to be consistent between GMT and PyGMT, and the best case is to bind to the GMT C function. But we can also reimplement it in PyGMT for now (in a separate PR), maybe in header.py?

Another possibility is to wrap the "hidden", private GMT_GRID_HEADER_HIDDEN structure https://github.com/GenericMappingTools/gmt/blob/7809736ba32d87a4a96b15444419eb176c6a35f3/src/gmt_hidden.h#L151. The structure has many members but we can only wrap the first few members. The members that are most useful to us are: grdtype (Cartesian or geographic), BC (boundary condition), varname (the actual variable name so that we don't have to hardcode the grid name to z), and cpt (the default CPT for this grid).

In f715aee, I've updated the logic for determining the grid/image gtype based on ProjRefPROJ4.

The logic of codes come from: https://github.com/GenericMappingTools/gmt/blob/7809736ba32d87a4a96b15444419eb176c6a35f3/src/gmt_grdio.c#L3778

weiji14 · 2024-09-29T03:50:21Z

pygmt/datatypes/image.py

+            title:
+            history:
+            description:


If those are only needed for CF-1.7, then maybe we can indent the 3 attrs["..."] lines here to be inside the if self.type in {GridFormat.NX} block?

pygmt/pygmt/datatypes/header.py

Lines 207 to 217 in cf48764

if self.type in {

GridFormat.NB,

GridFormat.NS,

GridFormat.NI,

GridFormat.NF,

GridFormat.ND,

}: # Only set the 'Conventions' attribute for netCDF.

attrs["Conventions"] = "CF-1.7"

attrs["title"] = self.title.decode()

attrs["history"] = self.command.decode()

attrs["description"] = self.remark.decode()

(2) actual_range is the standard attribution defined by CF-1.7 convention, but it's only for netCDF files, not for images;

Same here, we can put attrs["actual_range"] under the if-block above:

pygmt/pygmt/datatypes/header.py

Line 223 in cf48764

attrs["actual_range"] = np.array([self.z_min, self.z_max])

…T codes

weiji14

Wonderful, thanks again @seisman! Excited to have direct 3-band read support in PyGMT without rioxarray! There's still a few tiny details that could be improved (around the attrs/metadata), but this looks really good already as an initial implementation. 🚀

Wrap GMT's standard data type GMT_IMAGE for images

bcf43f0

seisman assigned weiji14 Mar 20, 2024

seisman added the enhancement Improving an existing feature label Mar 20, 2024

seisman mentioned this pull request Mar 20, 2024

Wrap GMT's standard data type GMT_GRID for grids #2398

Merged

5 tasks

Initial implementation of to_dataarray method for _GMT_IMAGE class

a052a1a

Minimal xarray.DataArray output with data and coordinates, no metadata yet.

seisman commented Mar 20, 2024

View reviewed changes

pygmt/clib/session.py Outdated Show resolved Hide resolved

weiji14 mentioned this pull request Apr 1, 2024

GMT_GRID_HEADER: Parse grid header and add grid properties #3134

Merged

Merge branch 'main' into datatypes/gmtimage

56a6d65

seisman mentioned this pull request Apr 19, 2024

Add placeholder codes for wrapping GMT_IMAGE/GMT_CUBE #3181

Closed

weiji14 added 2 commits June 18, 2024 21:16

Merge branch 'main' into datatypes/gmtimage

f71e79c

Small typo fixes and add output type-hint for to_dataarray

4cce4a2

weiji14 reviewed Jun 18, 2024

View reviewed changes

pygmt/datatypes/image.py Outdated Show resolved Hide resolved

weiji14 added 2 commits June 18, 2024 21:47

Fix mypy error using np.array([0, 1, 2]) instead of np.arange

e02b650

Parse name and data_attrs from grid/image header

f3d4b1f

Extra metadata from the _GMT_GRID_HEADER struct.

weiji14 reviewed Jun 18, 2024

View reviewed changes

weiji14 added 4 commits June 20, 2024 21:36

Transpose array to (band, y, x) order and add doctest for to_dataarray

4390136

Reorder the dimensions to follow Channel, Height, Width (CHW) convention. Also added doctest checking output DataArray object and the image's x and y coordinates.

Set registration and gtype from header

5f25669

Get the registration and gtype info from the grid header and apply it to the GMT accessor attributes.

Print basic shape and padding info in _GMT_IMAGE doctest

a3c6c14

Trying to match some of the doctests in _GMT_GRID.

Only set Conventions = CF-1.7 attribute for NetCDF grid type

5888e10

Remove hardcoded attribute that was only meant for NetCDF files, so that GeoTIFF files won't have it.

weiji14 reviewed Jun 20, 2024

View reviewed changes

pygmt/datatypes/header.py Outdated Show resolved Hide resolved

weiji14 added 2 commits June 20, 2024 22:00

Merge branch 'main' into datatypes/gmtimage

798e658

Remove rioxarray import

3dbf2f2

weiji14 mentioned this pull request Jun 20, 2024

**BREAKING** pygmt.grdcut: Refactor to store output in virtualfiles for grids #3115

Open

seisman commented Jun 20, 2024

View reviewed changes

pygmt/clib/session.py Outdated Show resolved Hide resolved

seisman mentioned this pull request Jul 8, 2024

POC: Wrap GMT_Read_Data and read datasets/grids/images into GMT data container #3318

Closed

5 tasks

seisman commented Jul 19, 2024

View reviewed changes

pygmt/datatypes/image.py Outdated Show resolved Hide resolved

seisman mentioned this pull request Jul 19, 2024

Wrap GMT's standard data type GMT_IMAGE for images #3338

Merged

7 tasks

Merge branch 'main' into datatypes/gmtimage

6b860bf

seisman changed the title ~~WIP: Wrap GMT's standard data type GMT_IMAGE for images~~ WIP: Implemention GMT_IMAGE.to_dataarray method Jul 27, 2024

Merge branch 'ctypesarray' into datatypes/gmtimage

8433e78

seisman marked this pull request as draft September 20, 2024 15:19

Update the to_dataarray method to match the codes in GMT_GRID

3e3a6f3

seisman changed the base branch from main to ctypesarray September 20, 2024 16:02

seisman added 2 commits September 21, 2024 00:31

image data should has uint8 dtype

12ef40a

Further improve the tests

f64fbb8

Base automatically changed from ctypesarray to main September 24, 2024 02:10

seisman added 2 commits September 24, 2024 10:14

Merge branch 'main' into datatypes/gmtimage

4f2ae48

Add a note that currently only 3-band images are supported

d49afed

This was referenced Sep 24, 2024

Add enums GridFormat for GMT grid format ID #3449

Merged

Set the 'Conventions' attribute to 'CF-1.7' for netCDF grids only #3463

Merged

seisman added 2 commits September 28, 2024 21:18

Merge branch 'main' into datatypes/gmtimage

f70bec0

Remove the old GMTGridID enums from pygmt/datatypes/header.py

2fd13fb

seisman changed the title ~~GMT_IMAGE: Implement the GMT_IMAGE.to_dataarray method~~ GMT_IMAGE: Implement the GMT_IMAGE.to_dataarray method for 3-band images Sep 28, 2024

seisman added the needs review This PR has higher priority and needs review. label Sep 28, 2024

seisman marked this pull request as ready for review September 28, 2024 14:31

A minor fix

9972ba1

seisman commented Sep 28, 2024

View reviewed changes

weiji14 reviewed Sep 29, 2024

View reviewed changes

seisman added 2 commits September 29, 2024 12:28

Set CF-1.7-specific attributes for netCDF formats only

62f0ce0

Improve the logic for determine grid/image gtype based on upstream GM…

f715aee

…T codes

weiji14 approved these changes Sep 29, 2024

View reviewed changes

seisman removed the needs review This PR has higher priority and needs review. label Sep 29, 2024

seisman merged commit 6606845 into main Sep 29, 2024
20 of 21 checks passed

seisman deleted the datatypes/gmtimage branch September 29, 2024 08:11

seisman mentioned this pull request Sep 29, 2024

Add function to load Blue Marble dataset #2235

Merged

6 tasks

This was referenced Sep 29, 2024

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

Closed

Integrate with rioxarray and implement GMT_IMAGE to read RGB GeoTIFFs #1555

Open

Get projection/crs from GMT_GRID_HEADER #3470

Open

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

GMT_IMAGE: Implement the GMT_IMAGE.to_dataarray method for 3-band images #3128

GMT_IMAGE: Implement the GMT_IMAGE.to_dataarray method for 3-band images #3128

seisman commented Mar 20, 2024 •

edited

Loading

seisman commented Apr 14, 2024

weiji14 Jun 18, 2024

seisman Jun 18, 2024

seisman commented Sep 28, 2024 •

edited

Loading

seisman Sep 28, 2024

weiji14 Sep 28, 2024

seisman Sep 29, 2024

weiji14 Sep 29, 2024

seisman Sep 29, 2024

weiji14 left a comment

weiji14 Sep 29, 2024 •

edited

Loading

seisman Sep 29, 2024

weiji14 Sep 29, 2024

seisman Sep 29, 2024

weiji14 Sep 29, 2024

seisman Sep 29, 2024 •

edited

Loading

seisman Sep 29, 2024

weiji14 Sep 29, 2024

weiji14 left a comment

	@property
	def name(self) -> str:
	"""
	Name of the grid.
	"""
	return "z"

	if self.type in {
	GridFormat.NB,
	GridFormat.NS,
	GridFormat.NI,
	GridFormat.NF,
	GridFormat.ND,
	}: # Only set the 'Conventions' attribute for netCDF.
	attrs["Conventions"] = "CF-1.7"
	attrs["title"] = self.title.decode()
	attrs["history"] = self.command.decode()
	attrs["description"] = self.remark.decode()

GMT_IMAGE: Implement the GMT_IMAGE.to_dataarray method for 3-band images #3128

GMT_IMAGE: Implement the GMT_IMAGE.to_dataarray method for 3-band images #3128

Conversation

seisman commented Mar 20, 2024 • edited Loading

seisman commented Apr 14, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

seisman commented Sep 28, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

weiji14 left a comment

Choose a reason for hiding this comment

weiji14 Sep 29, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

seisman Sep 29, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

weiji14 left a comment

Choose a reason for hiding this comment

seisman commented Mar 20, 2024 •

edited

Loading

seisman commented Sep 28, 2024 •

edited

Loading

weiji14 Sep 29, 2024 •

edited

Loading

seisman Sep 29, 2024 •

edited

Loading