Plotting a georeferenced field on a 3d globe

[aasdelat]

I have georeferenced data in longitude/latitude, and I want to plot them on a three-dimensional globe with the WGS84 datum, using the GLMakie and GeoMakie libraries.

You can download the data from here.

The grid is irregular in longitude-latitude. Each georeferenced data point corresponds to a cell whose corners or vertices have known longitude and latitude coordinates. I want each cell to be colored with a constant color (no color gradients) that corresponds to the data of the grid point in the cell.

asinghvi17 gave me a code that makes it but with non-constant color and with no cells, and I wonder even if it is possible what I want to do.

I have slightly modified it, but I get an error. See the comments in the code:

using NCDatasets # io
using GLMakie, GeoMakie # visualization
using Delaunator # fast meshing of 2d point clouds
using ProgressMeter

ds = NCDataset("/home/antonio/Documentos/codigo/grid_utils/N320_weights.nc")

## From here =======================================================
## Code that runs, but no constant cell color:
#t1 = Delaunator.triangulate(tuple.(ds["longitude"], ds["latitude"]))
#fig, ax, plt = mesh(
#    t1.points, GLMakie.TriangleFace.(t1.triangles); 
#    color = replace(ds["spatial_weights"][:, 1], missing => NaN),
#    shading = NoShading, # flat lighting on the mesh 
#    axis = (; type = GeoMakie.GlobeAxis, show_axis = false),
#    figure = (; backgroundcolor = :black)
#)
## Until here, the code that runs
## Until here ========================================================


# From here =======================================================
# From here, the code that does not run. It is an attemt to make cell constant
#colors:

# For each grid point we have the cell surrounding it that is delimited by
# the longitudes bounds and latitudes bounds of that grid point.
# So, there is a two dimensional array. One dimension corresponds to the grid
# point, so that it is as long as the number of grid points. The other dimension
# is two, so that the elemen [1,1] is the upper latitude of the first cell and
# [1,2] is the lower latitude of the same cell. For the grid point number 238,
# the latitude boundaries are at [238,1] (upper) and [238,2] (lower), and so on.
# And analogously for the longitudes.
# So, we form the vertices and the faces as follows:
vertex = []
faces  = []
@showprogress for i in 1:length(ds["longitude"])
    push!(   vertex, [ ds["lon_bnds"][i,1], ds["lat_bnds"][i,1] ]   )
    push!(   vertex, [ ds["lon_bnds"][i,1], ds["lat_bnds"][i,2] ]   )
    push!(   vertex, [ ds["lon_bnds"][i,2], ds["lat_bnds"][i,2] ]   )
    push!(   vertex, [ ds["lon_bnds"][i,2], ds["lat_bnds"][i,1] ]   )
    
    vertex_len = length(vertex)
    push!( faces, [ vertex_len-3, vertex_len-2, vertex_len-1 ])
    push!( faces, [ vertex_len-3, vertex_len-1, vertex_len   ])
end # for
vertex = vcat(vertex'...)
faces  = vcat(faces'...)

fig, ax, plt = mesh(
    vertex, faces; 
    color = replace(ds["spatial_weights"][:, 1], missing => NaN),
    shading = NoShading,  flat lighting on the mesh 
    axis = (; type = GeoMakie.GlobeAxis, show_axis = false),
    figure = (; backgroundcolor = :black)
)
# Throws the error:
#ERROR: Buffer vertices does not have the same length as the other buffers.
#	vertex_color has length 542080
#	vertices has length 2168320
# Until here, the attempt that does nor run
# Until here ========================================================


# Plot coastlines, color black with 10% opacity, 20km above sea level
coastlines_plt = lines!(ax, GeoMakie.coastlines(); color = (:black, 0.1), zlevel = 20_000)

# Display the figure so you can see it and interact with it!
# GlobeAxis user controls are not the greatest at the moment but that
# will improve this summer.
fig

# Animate the whole thing
## first, add a title to the top
#title = Label(fig[0, 1]; text = "N320 spatial weights", tellwidth = false, tellheight = true, color = :white)
#sw = ds["spatial_weights"]
## then, simply loop over the frames and update the color of the mesh plot, and the title!
## Animation in Makie is really simple using Observables.
#record(fig, "over_time.mp4", 1:length(ds["valid_time"]), framerate = 60, px_per_unit = 2, update = false) do i
#    plt.color[] = replace(msl[:, i], missing => NaN)
#    title.text[] = string(ds["valid_time"][i])
#end

[asinghvi17]

I’ll take a look! You can ping people so they know about your post on Discourse like so: @asinghvi17

but in general I’m more active on Slack and Github, just FYI

So it looks like you have quad faces here that you want to visualize. IIRC GeometryBasics should just support those and we have a FaceView attribute now to provide per face meshing. I’ll try to come up with an example in a bit.

@Raf just FYI, this is another file that just doesn’t load the bounds in Rasters…it’s not fully CF compliant though. @aasdelat where are you getting these from? If they are manually generated you should consider adding a CF “coordinates” attribute to your value array otherwise it’s a bit hard to figure out what is going on…

[asinghvi17]

Off topic, it is useful when performing computations, to materialize all your data in Julia beforehand. So instead of referring to ds["lon_bnds"][i,1] in a loop, you should instead collect it into Julia:

lon_bnds = ds["lon_bnds"][:, :]

then accessing it should be infinitely faster. That’s the reason your loop is so slow right now - it has to do a C library call and then decompression for every indexing event on the dataset!

If we get Rasters working with this data it will read it in automatically if the data is small enough

[asinghvi17]

also, I’m not sure your data even makes sense when plotted this way. All the bounding boxes seem to intersect quite a bit, here is the first one:

so in general for data in this format there is a way to plot it, but I have no idea how this got generated in the first place. The plot doesn’t make much sense because the bounding boxes overlap so much.

For what it’s worth, here is the code:

using GeometryBasics
using GLMakie, GeoMakie
using NCDatasets

ds = NCDataset("/Users/anshul/Downloads/N320_weights.nc")
# Access the bounds variables and "flatten" them, such that 
lons = ds["lon_bnds"][:, :]
lats = ds["lat_bnds"][:, :]
# Initialize vectors but don't explicitly store anything in them yet
# these will be filled with random bits right now
vs = Vector{Point2d}(undef, 4 * length(ds["longitude"]))
fs = Vector{QuadFace{Int}}(undef, length(ds["longitude"]))
# Now we can fill the vectors with the values we are interested in.
@showprogress Threads.@threads for i in 1:length(ds["longitude"])
    idx = (i - 1) * 4 + 1 # 4 points per cell, so we store 4 points in `vs` for every point in `ds["longitude"]`
    # Store the 4 points that make up the cell in `vs`
    vs[idx] = Point2d(lons[i,1], lats[i,1])
    vs[idx + 1] = Point2d(lons[i,1], lats[i,2])
    vs[idx + 2] = Point2d(lons[i,2], lats[i,2])
    vs[idx + 3] = Point2d(lons[i,2], lats[i,1])
    # Store the 4 faces that make up the cell in `fs`
    fs[i] = QuadFace(idx, idx + 1, idx + 2, idx + 3)
end # for

msh = GeometryBasics.Mesh(
    vs, fs; 
    color = FaceView(
        ds["spatial_weights"][:],                  # the data to color by
        QuadFace.(1:length(ds["spatial_weights"])) # the faces to apply it to (as indices)
    )
)
msh = GeometryBasics.normal_mesh(msh)

mesh(msh; axis = (; type = GeoMakie.GlobeAxis), shading = NoShading)


# Reference plot, just to make sure what we are doing is correct.
using Delaunator
tri = Delaunator.triangulate(Point2d.(ds["longitude"][:], ds["latitude"][:]))
mesh(tri.points, TriangleFace.(tri.triangles); color = ds["spatial_weights"][:], axis = (; type = GeoMakie.GlobeAxis), shading = NoShading)

[aasdelat]

Thank you, @asinghvi17, I was not sure if you nicname was the same in Julia Discourse as in GitHub.

It is me who has generated the file, so I have ammended it according to your indications (I hope). I have added the coordinates = "longitude latitude" attribute to the spatial_weight NetCDF variable. I have also changed the name of the NetCDF variable form spatial_weights (plural) to spatial_weight (singular).

The new file can be downloaded form here (I cannot edit my first post in order to update the link).

[aasdelat]

Thank you for your advice, @asinghvi17.

Perhaps we have to take into account that the cells corresponding to the northernmost and southernmost latitudes, are not quadrilaterals, but triangles, so they have a repeated vertex.

So, for example, the first grid point, which is in the northernmost latitude circle, has latitude bounds: 90.0, 89.6455; and longitude bounds: -10.0, 10.0.
The vertices (that your code generates very well) are:

 [-10.0, 90.0]
 [-10.0, 89.6455398227125]
 [10.0, 89.6455398227125]
 [10.0, 90.0]

where the first element is the longitude and the second one is the latitude.

But you can note that the points [-10.0, 90.0] and [10.0, 90.0], are the same, because both of them are exactly in the North Pole. And any point with 90º latitude is at the North pole regardless of its longitude, and the same applies to the South Pole for points with -90º latitude.

Is this a problem for the making of the graphic?

I have run your code (by first commenting out the last lines that draw the old triangle mesh) and I can see the sphere very well, and I do not see any defect on it, as the ones you show in one of the images in your post.

Is there a way to draw the wires of the mesh?
And also, can I apply a custom colormap?
The help on the Mesh object is not very clear for me. Is there a friendly documentation about plotting meshes like the one we are playing with?
And, of course, thank you very much for your help, @asinghvi17.