Still haven’t found the time! 
It’s on my to-do list!
Thanks, and it can even be made nicer if one translates this example with SSTs to Julia. It would be almost a ono-to-one operation.
The thing I wish GMT.jl could do is hide this complex syntax from the user by wrapping the low-level gmt commands into a Julian API. The raw commands are too low-level for daily usage. I came to the conclusion that nothing can beat the high-quality map visualizations that gmt produces, but only a few people are capable of generating them.
4 Likes
I mentioned that before. Given the level of details that GMT offers to parameterize the plots/calculus it’s virtually impossible to completely replace its condensed syntax by a verbose one. But for the more common options we can (and have done) do it. I see that you have not consulted the GMT.jl man and examples lately.
Our (long, ~30 years) historical record tells us that many tens of thousands of people have learn to do it.
1 Like
Just sharing what I think as someone that is not following the GMT trends nor history 
if it is a tool for experts on GMT only, that is fine. But it is great to hear that some parts of it are being ported to a more Julian API.
1 Like
Hi @naufalwx,
I would like to share an update regarding structured grids in GeoStats.jl.
Here is how you can load and plot the data:
using GeoStats
using NetCDF
using Plots
# load 2D matrices containing data
X = ncread("Matthew.nc", "x");
Y = ncread("Matthew.nc", "y");
P = ncread("Matthew.nc", "sfc_precip")
# spatial data
sdata = StructuredGridData(Dict(:precipitation => P), X, Y)
# may take some time depending on the Plots.jl backend
# good idea to set gr(format=:png) for quick plots
plot(sdata, ms=0.1, size=(1000,1000))
Notice that this plot recipe is sub-optimal. It is just plotting the structured grid data as a collection of points with color. This is a fallback recipe for all spatial data types in GeoStats.jl. We cannot do better until Plots.jl backends implement something like pcolormesh from matplotlib.
You can plot subsets of spatial data (or any spatial object in GeoStats.jl) for example:
plot(view(sdata, 1:1000), ms=0.1, size=(500,500))
In the example above I am plotting the first 1000 points in the satellite trajectory. You can use view to take any subset of points in the range 1:npoints(sdata).
Now we can start defining our estimation problem. Say for example we want to estimate values in a regular grid covering the data. We first define the grid with:
# collect coordinates extrema
xmin, xmax = extrema(X)
ymin, ymax = extrema(Y)
# spatial domain
sdomain = RegularGrid((xmin,ymin), (xmax,ymax), dims=(50,50))
plot(sdomain, ms=0.1, size=(500,500)
plot!(sdata, ms=0.1)
After the estimation domain is defined, you can pick any solver from the framework that is compatible with the domain type to solve the estimation problem:
problem = EstimationProblem(sdata, sdomain, :precipitation)
Please let me know if something is not clear. The codebase is under heavy development.