GenericMappingTools / gmt

The Generic Mapping Tools
https://www.generic-mapping-tools.org
Other
858 stars 359 forks source link

Surface does not accept Cartesian values for spacing #8476

Open newnorsk opened 6 months ago

newnorsk commented 6 months ago

Description of the problem

The "surface" function does not seem to work with Cartesian (projected) values for the spacing of the new grid, e.g. "e" meters, "k" kilometers. The documentation seems to indicate that it should, and indeed, the blockmean etc. functions are tolerant of projected spacing parameters.

Minimal Complete Verifiable Example

import pygmt

# This example taken almost wholesale from the PyGMT Refresher Webinar
# Load a table of ship bathymetric observations off Baja California
data = pygmt.datasets.load_sample_data('bathymetry')
# Use pygmt.info to store the x/y range of data values rounded to the nearest degree
region = pygmt.info(data, spacing=1)

# This works:
# Compute the median position and value for 10 arc-minute blocks
spacing = "10m"
data_median = pygmt.blockmedian(data, region=region, spacing=spacing)
# This was in the tutorial; I think its no longer necessary, but if it ain't broke . . .
data_median_np = data_median.to_numpy()
# Produce a grid from the x,y,z data using pygmt.surfae with 10 arc-minute grid spacing
grid = pygmt.surface(data=data_median_np, region=region, spacing=spacing)

# This doesn't work:
# Compute the position and value for 20 km blocks: approx. same resolution as 10 arc-minute
spacing = "20k"
data_median = pygmt.blockmedian(data, region=region, spacing=spacing)
data_median_np = data_median.to_numpy()
# Produce a grid from the x,y,z data using pygmt.surfae with 20 km grid spacing
grid = pygmt.surface(data=data_median_np, region=region, spacing=spacing)

Full error message

surface [WARNING]: (x_max-x_min) must equal (NX + eps) * x_inc), where NX is an integer and |eps| <= 0.0001.
surface [WARNING]: (y_max-y_min) must equal (NY + eps) * y_inc), where NY is an integer and |eps| <= 0.0001.
surface (gmtapi_init_grdheader): Please select compatible -R and -I values
surface [ERROR]: Grid must have at least 4 nodes in each direction (you have 2 by 2) - abort.

---------------------------------------------------------------------------
GMTCLibError                              Traceback (most recent call last)
Cell In[9], line 9
      6 data_median_np = data_median.to_numpy()
      8 # Produce a grid from the x,y,z data using pygmt.surfae with 20 km grid spacing
----> 9 grid = pygmt.surface(data=data_median_np, region=region, spacing=spacing)

File /opt/miniconda3/envs/pygmt/lib/python3.12/site-packages/pygmt/helpers/decorators.py:603, in use_alias.<locals>.alias_decorator.<locals>.new_module(*args, **kwargs)
    596     msg = (
    597         "Parameters 'Y' and 'yshift' are deprecated since v0.8.0. "
    598         "and will be removed in v0.12.0. "
    599         "Use Figure.shift_origin(yshift=...) instead."
    600     )
    601     warnings.warn(msg, category=SyntaxWarning, stacklevel=2)
--> 603 return module_func(*args, **kwargs)

File /opt/miniconda3/envs/pygmt/lib/python3.12/site-packages/pygmt/helpers/decorators.py:776, in kwargs_to_strings.<locals>.converter.<locals>.new_module(*args, **kwargs)
    773             bound.arguments["kwargs"][arg] = newvalue
    775 # Execute the original function and return its output
--> 776 return module_func(*bound.args, **bound.kwargs)

File /opt/miniconda3/envs/pygmt/lib/python3.12/site-packages/pygmt/src/surface.py:168, in surface(data, x, y, z, **kwargs)
    166         if (outgrid := kwargs.get("G")) is None:
    167             kwargs["G"] = outgrid = tmpfile.name  # output to tmpfile
--> 168         lib.call_module(
    169             module="surface", args=build_arg_string(kwargs, infile=infile)
    170         )
    172 return load_dataarray(outgrid) if outgrid == tmpfile.name else None

File /opt/miniconda3/envs/pygmt/lib/python3.12/site-packages/pygmt/clib/session.py:624, in Session.call_module(self, module, args)
    620 status = c_call_module(
    621     self.session_pointer, module.encode(), mode, args.encode()
    622 )
    623 if status != 0:
--> 624     raise GMTCLibError(
    625         f"Module '{module}' failed with status code {status}:\n{self._error_message}"
    626     )

GMTCLibError: Module 'surface' failed with status code 79:
surface [ERROR]: Grid must have at least 4 nodes in each direction (you have 2 by 2) - abort.

System information

PyGMT information:
  version: v0.11.0
System information:
  python: 3.12.2 | packaged by conda-forge | (main, Feb 16 2024, 21:00:12) [Clang 16.0.6 ]
  executable: /opt/miniconda3/envs/pygmt/bin/python
  machine: macOS-10.15.7-x86_64-i386-64bit
Dependency information:
  numpy: 1.26.4
  pandas: 2.2.1
  xarray: 2024.3.0
  netCDF4: 1.6.5
  packaging: 24.0
  contextily: 1.6.0
  geopandas: 0.14.3
  ipython: None
  rioxarray: 0.15.5
  ghostscript: 10.03.0
GMT library information:
  binary version: 6.5.0
  cores: 12
  grid layout: rows
  image layout: 
  library path: /opt/miniconda3/envs/pygmt/lib/libgmt.dylib
  padding: 2
  plugin dir: /opt/miniconda3/envs/pygmt/lib/gmt/plugins
  share dir: /opt/miniconda3/envs/pygmt/share/gmt
  version: 6.5.0
welcome[bot] commented 6 months ago

👋 Thanks for opening your first issue here! Please make sure you filled out the template with as much detail as possible. You might also want to take a look at our contributing guidelines and code of conduct.

seisman commented 6 months ago

I can reproduce the issue using the GMT CLI, so it's more likely an upstream bug:

gmt surface @tut_ship.xyz -R245/255/20/30 -I20k -Gabc.nc

I'm transferring the bug report to the upstream GMT repository.

welcome[bot] commented 6 months ago

👋 Thanks for opening your first issue here! Please make sure you filled out the template with as much detail as possible. We appreciate that you took the time to contribute!

Please make sure you read our Contributing Guide and abide by our Code of Conduct.