Spatial analysis for a list of coordinates

Q&A
1

Dear GMT experts,

I have a list of coordinates and need average values for altitude, precipitation, etc. that are stored in grid files. I would like to calculate the average values for a circular window, with the coordinates representing the center point.

My current workflow is a combination of grdfilter and grdtrack. However, grdfilter becomes slow with large grids and large window sizes. This is inconvenient because the grid data and coordinate lists represent the solutions of time-dependent models and the workflow has to be executed very often (> 1000 times for 2000 km x 2000 km grids and window sizes > 50 km).
Perhaps I overlooked it among the many functions—is there a flag in grdtrack that allows a search radius and spatial statistics, or another tool that could speed up the workflow?

Thank you very much for your answers.
Jörg

2

Hi, this is how I would probably do it with GMT.jl
It takes a polygon (the circle in your case) and uses gmtselect (or similar) to find the x,y,z grid points that are inside that polygon and than do some stats with them.The following example computes the average height of Swisserland as given by the 6m global grid.

using GMT

julia> G = gmtread("@earth_relief_10m");

julia> suiss= coast(DCW="CH", dump=true, minpts=300);

julia> zonal_statistics(G, suiss, mean)

BoundingBox: [1207.8861083984375, 1207.8861083984375]

1×1 GMTdataset{Float64, 2}
 Row │       X
─────┼─────────
   1 │ 1207.89
3

Hi Joaquim,

I still use GMT’s Classic Mode. That’s probably the legacy of over 20 years of GMT and code snippets for everything imaginable.

My first thought was also to use gmtselect and calculate the statistics from there. However, writing so many files put me off. With several hundred coordinates per time step and five raster data sets for the statistics, 100 time steps per scenario and 64 scenarios, that’s quite a lot of files that need to be written. Do you think the julia approach could be faster than grdfilter? The grdfilter commands are very easy to parallelize, but as far as I understand, your calculation does not require writing temporary files because the data remains in RAM?

Our IT department blocks julia codes on Windows computers, but I could try running julia - gmt on the Linux cluster if it pays off in terms of performance.

Thanks for your solution suggestion
Jörg

4

The scenario you describe seems perfectly suited for the type f solution I mentioned. Though ofc, details may prove to be a bit more hard to solve. It should be a very fast procedure as everything is in fact done with in memory only.

1- Yes, it can be run in Linux
2- I don’t believe taht they can stop you for using JUlia on Windows. Even if you are not allowed to run an installer, you can always use the portable version and add the ...\julia\bin to the path.

Note, whatever you go, Linux or Win, always insrtall only the Julia 1.10.? version. Later versions are much slower to start.

1 Like