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)