Grid Interpolation of Data from Polar-Orbiting Satellites

@naufalwx can you share the data so that I can try to reproduce the interpolation with GeoStats.jl? I am still confused about what you want to achieve. It looks like we can do it already with the current version of the project.

Please feel free to reach out in our gitter channel: https://gitter.im/JuliaEarth/GeoStats.jl

Following the documentation of Scipy’s griddata, they first generate a tesselation for the given input knots and then interpolate for the target points. I think the building blocks to do this in Julia are there, you can make use of

to create the tesselation. The package also provides a locate function to tell you which triangle your target point is in. The only remaining step would then be the linear interpolation inside the triangle.

In a single step by GMT’s triangulate that is available via GMT.jl

I don’t think there is even a need to perform tesselation in this case. The grid he is plotting is apparently a masked regular grid? If that is the case, this could be literally solved in 4 lines of GeoStats.jl

Anyways, let’s wait for a reply, it may take a couple of days again…

Hi everyone,

Thank you for getting back to me so soon. Here is an example file: https://github.com/naufalwx/TC-Microwave-Database/blob/master/Matthew_test.nc

My initial plan was to try to do a bilinear interpolation. But, I don’t know how, given the data structure. So I looked at packages in other languages and figured out I could treat the data points as scattered data points, and use another interpolation method instead of bilinear, which brought me to scipy’s griddata.

@naufalwx can you explain how the data can be loaded in Julia? Tried FileIO but it failed.

@juliohm It’s a NetCDF file. So you should just be able to use ncread. The three variables on there are labeled “sfc_precip”, “x”, and “y” :slight_smile:

Sorry but this file has no valid data in neither of its 3 vars: x, y, z

Could you maybe share a small script to load the data? Help us help you :slight_smile:

Definitely :slight_smile: I am able to read-in the file in the same kernel that I wrote the file, but if I open a new kernel and try to read-in the same file, it failed. Give me just a moment to figure it out and I’ll send out another notification when I’ve resolved this issue. Sorry for the delay, but thank you for your continued interest in helping me! :slight_smile:

I’ve resolved the issue. The updated file should be on the GitHub link. To read-in the file:


using NetCDF;

x = ncread("Matthew_test.nc", "x");
y = ncread("Matthew_test.nc", "y");
precip = ncread("Matthew_test.nc", "sfc_precip");

# If you want to look at the variable attributes
ncinfo("Matthew_test.nc")

Just make sure that the file is in the right directory :slight_smile: All files should have the dimension of (221, 366). There are NaNs in “sfc_precip”.

Thanks, now I get the format. There is a feature request open for this NetCDF convention: https://github.com/juliohm/GeoStats.jl/issues/22 @Balinus was working on it, but I don’t know he moved forward with the implementation.

The issue will be addressed in future releases.

1 Like
using GMT

# Assume that the Matthew_test.nc is at current dir

# Load the x,y, z arrays from file
x = gmt("read -Tg Matthew_test.nc?x");
y = gmt("read -Tg Matthew_test.nc?y");
z = gmt("read -Tg Matthew_test.nc?sfc_precip");

# Check the range of the coordinates. Need to know x_min, x_max, y_min, y_max
# Those values are the 5th and 6th of
x.range
y.range

# Compute a grid. Don't know original point separation so from dimensions and x values
# esrtimated that they are about 4 km.
G = gmt("nearneighbor -R-1428/1215/-2729/2326 -I4 -S20", [x.z[:] y.z[:] z.z[:]]);

# Quick image
grdimage(G, shade="+ne0.8+a100", proj="X10c/0", frame="a", fmt="jpg", show=true)

3 Likes

GMT visualizations look amazing as usual.

Still haven’t found the time! :confused: 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 :slight_smile: 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.