Pygmt.surface - issue with spacing

Q&A PyGMT Q&A
1

Hi all! I’m new to using pygmt and have ran into an issue trying to plot using pygmt.surface() using spacing in a format other than degrees. I am trying to plot the following gridded modeled data that has a 4km x 4km resolution:

image
image669×509 62.7 KB

To do this, I first use pygmt.blockmean(), using the data’s lon, lat, and NOx concentrations as my x,y, and z, and spacing= [‘4000e’,‘4000e’]. Then I use the blockmean output as ‘data’ in the pygmt.surface() function, along with my specified region, and the spacing defined as [‘4000e’,‘4000e’]. When doing this I receive the following error:

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

I’ve tried changing the spacing to be a larger resolution (e.g, 16000e,16000e and 5.0k,5.0k) and smaller. I’ve also updated to the most recent version of pygmt. However, when I rerun the script using spacing in the unit degrees, the script works just fine. I believe the pygmt.surface() function is not registering that I am passing a spacing in meters opposed to degrees, and it results in this error. Is there any type of fix for this? Thank you in advance!

2

Hello @vlang,

Welcome to the GMT forum :slightly_smiling_face:!

Can you please provide your code and ideally your data? This will make it significantly easier for people to help you out. You can format your script as code by placing three backticks in the line before and after the block with the script:

```
your script formated as code
```
3

Thank you, @yvonnefroehlich! Unfortunately the data I am working with is quite large but I can provide a working example of it.

data = xr.open_dataset('data.nc')
## print statements of data 
data.XLONG
xarray.DataArray'XLONG'south_north: 1008west_east: 1332
array([[-121.76193 , -121.72748 , -121.69299 , ...,  -73.55261 ,  -73.51788 ,
         -73.483154],
       [-121.771576, -121.73709 , -121.70261 , ...,  -73.54343 ,  -73.5087  ,
         -73.47397 ],
       [-121.78122 , -121.74672 , -121.71222 , ...,  -73.53424 ,  -73.49951 ,
         -73.46475 ],
       ...,
       [-137.24739 , -137.1959  , -137.14438 , ...,  -58.728455,  -58.675934,
         -58.623474],
       [-137.27182 , -137.2203  , -137.16876 , ...,  -58.704895,  -58.652374,
         -58.599884],
       [-137.2963  , -137.24475 , -137.19318 , ...,  -58.681335,  -58.628784,
         -58.576263]], dtype=float32)
data.XLAT
xarray.DataArray'XLAT'south_north: 1008west_east: 1332
array([[18.191376, 18.200508, 18.209648, ..., 18.530296, 18.521599, 18.512878],
       [18.224121, 18.233261, 18.242393, ..., 18.563217, 18.554527, 18.545807],
       [18.256859, 18.266014, 18.275154, ..., 18.596169, 18.587456, 18.578735],
       ...,
       [51.81338 , 51.828457, 51.84353 , ..., 52.373074, 52.358704, 52.34431 ],
       [51.845192, 51.860294, 51.875366, ..., 52.405136, 52.39076 , 52.376354],
       [51.877014, 51.892117, 51.907192, ..., 52.43719 , 52.42281 , 52.4084  ]],
      dtype=float32)
data.NOx_all
xarray.DataArray'NOx_all'Time: 1south_north: 1008west_east: 1332
array([[[1., 1.2., 1., ..., 1., 1.6., 1.3.],
        [1.2., 1.3., 1.4., ..., 1.4., 1.6., 1.],
        [0., .75., .4., ..., .4., .6., .1.],
        ...,
        [1.2., 1.3., 1.4., ..., 1.4., 1.6., 1.],
        [0., .75., .4., ..., .4., .6., .1.],
        [1., 1.2., 1.2., ..., 1.5., 1.5., 1.5.]]], dtype=float32)

# Plotting of data file
space = ['4000e','4000e']
nox_blockmean=pygmt.blockmean(x=np.array(data.XLONG).flatten(),
                       y=np.array(data.XLAT).flatten(),
                       z=np.array(data.NOx[0]).flatten(),
                        region= [data.XLONG.min(), data.XLONG.max(), data.XLAT.min(), data.XLAT.max()] ,
                       spacing= space)

nox_surf= pygmt.surface(data=nox_blockmean,
                       region=[data.XLONG.min(), data.XLONG.max(), data.XLAT.min(), data.XLAT.max()],
                        #spacing= 0.05)  #does not error when using degrees 
                       spacing= space)  # errors when using space in meters or km

I hope this was easy enough to follow. I appreciate any suggestions that can be offered. Thank you!

-Victoria

4

You say your data is already gridded at 400x400m. If so, why grid it again? Maybe just reformat your table into a grid?

5

Thank you for the reply. The grid I am using is curvilinear. I am unable to plot the data with grdimage because my map projection differs from the coordinate system in which my data was gridded.

6

OK, what you do is fine. I wonder if puGMT’s surface method handles the increment correctly. Seems like [4000e, 4000e] is seen as 4000,4000 (degrees) which would explain why you get a 1x1 grid.