Blockmean w/ x,y coordinates in km

Q&A PyGMT Q&A
1

I’m attempting to use blockmean in a workflow to produce a heatmap of earthquakes within a cross-section. The input data is a projected x-coordinate in km and a depth in km.

Here is some example data (both in km):
projx = [5.03380,4.73689,5.06110,4.85559,4.87003,2.48495,3.02243,4.50333,3.91557,3.67545,4.46045]
depth = [3.900,2.420,4.300,8.480,6.340,3.000,0.810,-0.120,2.100,-0.670,4.970]

Ideally, I’d like to use a spacing in x and y of 0.5 km (spacing=‘0.5k’) and count the number of earthquakes within each 0.5km x 0.5km block. The pygmt command looks like this:

hmgrid = pygmt.blockmean(x=projx,
y=depth,
z=np.zeros(len(depth)),
spacing=‘0.5k’,
region=’‘0/5.55/-4/15’,
S=‘n’)

The command runs successfully, but the resulting xyz values are confusing. I expect a grid of points every 0.5 km within the stated region, but I get something much finer.

I think this line in the blockmedian documentation is critical, " e (meter), f (foot), k (km), M (mile), n (nautical mile) or u (US survey foot), in which case the increment will be converted to the equivalent degrees longitude at the middle latitude of the region". I’d like to avoid that conversion all together…

So I attempted to define the incoming columns with coltypes=‘ic’, but I still don’t get the expected result. I also tried to define the region explicitly with a ‘+uk’ at the end of the region string without success.

Details:
pygmt version 0.5.0
gmt version 6.3.0

Thanks in advance for any help,

Wes

2

The increment units etc are only used for geographic (lon,lat) data. Once you have projected data (km, meters, inches) they are essentially a Cartesian case. So you should use 0.5 as the Cartesian increment. Since your data is not in longitude and latitude the section you copied do not apply to your case.

1 Like