GenericMappingTools / GenericMappingTools/pygmt
Problem using option maxradius in pygmt.surface
- Dominant language
- Python
- Stars
- 874
- Forks
- 255
- Avg merge
- 1d 21h
- Merged PRs (30d)
- 40
Description
### 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
```python
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
```python-traceback
There was either no message or "The kernel for test.ipynb appears to have died. It will restart automatically."
```
### System information
```bash
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
```
Contributor guide
Assessment
This issue has not been assessed yet.