Parsing Q to grdimage

[lsawade]

Parsing ‘Q’ to ‘grdimage’

Combing the surface function, and grdimage, I was able to interpolate a grid and plot it within a ‘rectangular’ region.

Within that region only values above a certain values are important, so I set values in the output ‘xarray’ below that threshold to np.nan.

I expected parsing Q=True would omit plotting nans much like in the original implementation of GMT, but looking at the Baseplotting it’s implementation didn’t strike me immediately.

Are you willing to help implement and maintain this feature?

Yes, I would try, but I have 0 experience with GMT and it would probably take forever.

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

@lsawade It would be helpful if you could provide an example so that we can try to reproduce your issue and give some feedbacks.

[weiji14]

Do you mean that np.nan isn't treated as NaN properly by PyGMT?

I know Q=True works when passing NetCDF/GeoTiffs (see e.g. https://nbviewer.jupyter.org/github/weiji14/deepbedmap/blob/master/paper_figures.ipynb#Thumbnail-figures, there's some small 'white' spots), and I seem to have a vague recollection that NaNs in xarray.DataArrays don't work (i.e. are still coloured), but it's been a while ago and an example would help.

[lsawade]

Hi @seisman , hi @weiji14 ,

Let x, y, z be arbitrarily scattered, geolocated data (lon, lat, z, respectively) in a region.

1. Interpolate the data using pygmt.surface and set values below threshold to NaN

This will result in grid being output as an xarray, which is fine.

import numpy as np
import pygmt as gmt

## — something that generates lon, lat, z — ##

region = [-125,-60, 10,50]  # USA ~ish
grid = gmt.surface(x=lon, y=lat, z=illpoints, spacing=0.5, region=region)

grid.where(grid < 100, np.nan)

I double checked the xarray, it’s definitely populated with NaNs

2. Plot the data skipping without NaNs

import pygmt as gmt

cmap = "magma"
region = np.array([-125,-60, 10,50])
buffer = np.array([-20, +20, -10, +10])

# init figure
fig = gmt.Figure()
fig.basemap(region=region+buffer, frame=True, projection="Cyl_stere/0/-20/8i")

# Plot. background map
fig.grdimage('@earth_relief_15m', shading='+', t=70)

# Create colormap
gmt.makecpt(cmap="magma", series=[0, np.max(grid.data)], continuous=False, reverse=True, Q=True)

# Plot data on top
fig.grdimage(grid=grid, C=True, Q=True). 

fig.show()

Result:

Expected:

Basically nothing where the black lines are (Technically also nothing in that middle white splodge, but you get the idea).

I'm probably doing something wrong either way....

[lsawade]

Hi @seisman , hi @weiji14 ,

Forget everything... it does work. I assume that I didn't have to reassign the array... I was being a dum-dum..

Example with randomly generated data:

import numpy as np
import pygmt as gmt

# Generate random data within region
region = np.array([-125,-60, 10,50])
lon = np.random.rand(1000)*(region[1]-region[0]) + region[0]
lat = np.random.rand(1000)*(region[3]-region[2]) + region[2]
z = np.cos(lon/180*np.pi*4)**2 + np.sin(4*lat/180*np.pi)**2

# Interpolate
grid = gmt.surface(x=lon, y=lat, z=z, spacing=0.5, region=[-125,-60, 10,50])

# Set values to Nan
nangrid = grid.where(grid < 1., np.nan)

# Plot things...
cmap = "magma"

# To extend map boundaries
buffer = np.array([-20, +20, -10, +10])

# init figure
fig = gmt.Figure()
fig.basemap(region=region+buffer, frame=True, projection="Cyl_stere/0/-20/8i")
fig.grdimage('@earth_relief_15m', shading='+', t=70)

# Colormap
gmt.makecpt(cmap="magma", series=[grid.min().values, grid.max().values], continuous=False, reverse=True, Q=True)

# Plot grid 
fig.grdimage(grid=nangrid, C=True, Q=True, t=30)
fig.show()

Example output map:

My map with interpolated geolocated data:

Thanks anyways for the quick response! I appreciate it!

Feel free to close since I am the issue :D

[weiji14]

Ah cool, glad to see you got it working. That region+buffer trick is quite neat, I'll definitely be using that for my next pygmt map! I'll close this now, but I think it would be nice to alias Q=transparency at some point :wink:

[lsawade]

@weiji14 I agree, or Q=mask, since t already controls overall transparency.

Yeah, at some point I had to generated maps automatically where I would know where seismic station and event are so i would use their min/max locations + buffer as region. Definitely neat, when you dont know where the data is gonna be!

Thanks again!