Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Surface does not accept Cartesian values for spacing #8476

Open
newnorsk opened this issue May 6, 2024 · 3 comments
Open

Surface does not accept Cartesian values for spacing #8476

newnorsk opened this issue May 6, 2024 · 3 comments
Labels
bug Something isn't working

Comments

@newnorsk
Copy link

newnorsk commented May 6, 2024

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
@newnorsk newnorsk added the bug Something isn't working label May 6, 2024
Copy link

welcome bot commented May 6, 2024

👋 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
Copy link
Member

seisman commented May 6, 2024

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.

@seisman seisman transferred this issue from GenericMappingTools/pygmt May 6, 2024
Copy link

welcome bot commented May 6, 2024

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working
Projects
None yet
Development

No branches or pull requests

2 participants