Problem using option maxradius in pygmt.surface

[cbouliga]

Description of the problem

Hi,

I have been having difficulties using the option 'maxradius' in pygmt.surface. I had this problem using data with x and y data in projected coordinates (utm in my case) collected in a small area. The problems were either that the Kernel crashes for certain choice of sampling or that the masked area is erroneous, sometimes shifted to the north or to the south with respect to the data and sometimes totally wrong. The non-masked grid seems however correct. After several tests, the problem disappeared when i removed the averaged value of my coordinates so that the coordinates are centered on the point of coordinates (0,0) and are small (so this is a way to solve this problem). I have the impression this happen when the sampling is small with respect to the values of the coordinates (although the data are in a small area). So it could just be a numerical procession problem.

In order to explain and reproduce this errors with data that are already online, i used the sample bathymetry data from pygmt and made as if the data are in m and i added 500000m to the x and y coordinates so that the coordinates look like projected coordinates.

Here is below the script that i run.

Thank you in advance for your help, Claire.

Minimal Complete Verifiable Example

import numpy as np
import matplotlib.pyplot as plt
import pygmt

#Figure parameters
plt.rcParams['figure.dpi'] = 600

# Read bathy data
bathy = pygmt.datasets.load_sample_data(name='bathymetry')

# Select dataset
imin=74500
imax=78300
indices=np.arange(imin,imax+1,1)
bathy=bathy.loc[indices]

# Modify coordinates
bathy.longitude=bathy.longitude+500000
bathy.latitude=bathy.latitude+500000

# Sampling 
dl=0.1 # kernel crashes
dl=0.15

# Perform gridding of topography data
xmin = np.floor((np.min(bathy.longitude)-dl)/dl)*dl
xmax = np.ceil((np.max(bathy.longitude)+dl)/dl)*dl
ymin = np.floor((np.min(bathy.latitude)-dl)/dl)*dl
ymax = np.ceil((np.max(bathy.latitude)+dl)/dl)*dl
grid = pygmt.surface(data=bathy, spacing=dl, tension=0, region=[xmin,xmax,ymin,ymax],maxradius="0c")
#grid = pygmt.surface(data=bathy, spacing=dl, tension=0, region=[xmin,xmax,ymin,ymax]) #no problem when using this

# Plot
plt.subplot(1,2,1)
plt.scatter(np.array(bathy.longitude),np.array(bathy.latitude), c=np.array(bathy.bathymetry), cmap='jet',s=0.1)
plt.axis('scaled')
plt.colorbar(shrink=0.4)
plt.subplot(1,2,2)
XV, YV = np.meshgrid(grid.x.values, grid.y.values)
plt.contourf(XV, YV, grid.values,cmap='jet',levels=256) 
plt.colorbar(shrink=0.4)
plt.plot(XV.flatten(),YV.flatten(),'.k',markersize=0.1)
plt.plot(np.array(bathy.longitude),np.array(bathy.latitude), '.g',markersize=0.2)
plt.axis('scaled')

Full error message

There was either no message or "The kernel for test.ipynb appears to have died. It will restart automatically."

System information

PyGMT information:
  version: v0.13.0
System information:
  python: 3.12.3 (main, Sep 11 2024, 14:17:37) [GCC 13.2.0]
  executable: /home/bouligac/venv/processing/bin/python3
  machine: Linux-5.15.133.1-microsoft-standard-WSL2-x86_64-with-glibc2.39
Dependency information:
  numpy: 2.1.1
  pandas: 2.2.3
  xarray: 2024.9.0
  netCDF4: 1.7.1.post2
  packaging: 24.1
  contextily: 1.6.2
  geopandas: None
  IPython: 8.28.0
  rioxarray: 0.17.0
  gdal: None
  ghostscript: 10.02.1
GMT library information:
  version: 6.5.0
  padding: 2
  share dir: /usr/share/gmt
  plugin dir: /usr/lib/x86_64-linux-gnu/gmt/plugins
  library path: /usr/lib/x86_64-linux-gnu/libgmt.so
  cores: 8
  grid layout: rows
  image layout: 
  binary version: 6.5.0

[welcome[bot]]

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

I can't reproduce your issue on Linux. Your script produces the following image which I assume is correct:

[yvonnefroehlich]

I can't reproduce your issue on Linux. Your script produces the following image which I assume is correct:

Hm. As I understand the issue the script works for a sampling of dl=0.15:

# Sampling 
dl=0.1 # kernel crashes
dl=0.15

But if one uncomments dl=0.15 and dl=0.1 is used, the script crashes. However, it works when maxradius="0c" is omitted in pygmt.surface.

[seisman]

Thanks, I misunderstood the original post. Now I can reproduce it and will find some time to debug it.

[cbouliga]

Hi,

Thank you very much for your prompt reply. Sorry, my message may have not been clear. I obtain the same image as you. The problem is that the masking of the grid is shifted with respect to the data (when i use dl=0.15). If i use dl=0.1, it's worse since the kernel crashes. On the other hand if you remove the two following lines then the masking of the grid is well centered on the data:

bathy.longitude=bathy.longitude+500000 bathy.latitude=bathy.latitude+500000

So i have the impression that the issue arises only when the ratio of the spacing with respect to the coordinates become very small. The issue may have not come up before possibly because pyGMT is more often used for large scale studies.

Sincerely, Claire.

Claire Bouligand Maître de conférences

---

Institut des Sciences de la Terre Université de Grenoble Alpes CS 40700 38058 Grenoble cedex 9 France ☏ +33 4 76 63 51 77

De: "Dongdong Tian" @.> À: "GenericMappingTools/pygmt" @.> Cc: "CLAIRE BOULIGAND" @.>, "Author" @.> Envoyé: Mercredi 9 Octobre 2024 19:26:10 Objet: Re: [GenericMappingTools/pygmt] Problem using option maxradius in pygmt.surface (Issue #3500)

I can't reproduce your issue on Linux. Your script produces the following image which I assume is correct: [ https://github.com/user-attachments/assets/dbec3099-30c5-4c4c-b444-300734e15a0f | map.png (view on web) ]

— Reply to this email directly, [ https://github.com/GenericMappingTools/pygmt/issues/3500#issuecomment-2403791520 | view it on GitHub ] , or [ https://github.com/notifications/unsubscribe-auth/BL7BSYB7I5QQXOOQ3KKVVFDZ2XQUFAVCNFSM6AAAAABPVAK6W6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDIMBTG44TCNJSGA | unsubscribe ] . You are receiving this because you authored the thread. Message ID: @.***>