I am new to Julia and work on satellite remote sensing of the atmosphere and ocean. I am looking for some advice on mapping point and polygon (vector) data in geographic coordinates using Julia.
I produce relatively coarse resolution data on scales of 5 to 50 km. The swaths are irregular since the size of the sensor field-of-view changes across the flight track of the orbit. In other words, the data are not raster, but geographic data are in simple arrays (“raw” coordinates, if you will).
I want to make color maps of several scalar fields that I’ve generated on a irregular grid and do so interactively using a relatively fast plotting backend rather than write output to postscript.
My challenge has been finding a convenient way to plot my spatial fields using packages in the Julia ecosystem. Since my task is common in GIS, I looked for Julia Geo packages which map Polygons, or MultiPolygons. But those I found appear to expect shapefiles or GeoJSON inputs and I want to stay away from that.
I simply want to put lots of polygon data that are in simple arrays on the Earth efficiently. I’m guessing this will require putting my polygon arrays into one of the geo data structures that has been developed for Julia, but I can’t figure out how that should be done.
Does anyone have suggestions of how to using some of the available Julia packages to make maps in this way? Or perhaps there is an alternative? I’m happy to combine a few steps using low level interfaces to accomplish the task.
Hi @haffner , welcome to the community. I sympathize with your requirements and I think that the GeoStats.jl framework may suit your needs. Given that you are new to Julia, please consider watching this JuliaCon talk to get an idea of the current functionality:
Regarding your specific visualization needs, I understood that you want to plot data over a general unstructured mesh? The MeshViz.jl submodule was developed for that. Please read the test suite for examples of visualizations.
Hi, I confess that I didn’t fully understand your requests, namely all that polygonal data coming from satellites, but you may be interested into looking at RemoteS.jl
However, I have to confess that it doesn’t meet one of your conditions as it creates the images in postscript under the hood. (but for vector data, polygons, postscript/pdf gives the highest quality images).
@joa-quim,
I looked at the example for reading and mapping MODIS L2 data in the documentation on RemoteS.grid_at_sensor and that’s very nice. Sometimes I work with the final QC’d MODIS data. However I want to map the FOV of each satellite pixel like this instead,
With this type of figure I can see each measurement at original resolution which varies across the swath and with location in certain projections (distortion). It helps for detailed analysis because it’s a truer representation (spatially) of the retrieved field. Noise and outliers are easier to detect in this type of image than with a map made by interpolating between center points of the measurements.
Re: figure quality, I agree there’s no contest with postscript. I see RemoteS uses GMT and I sometimes use that to make ps figures for publications.
That’s good to know @Raf. I think I’ve been brainwashed by cartopy that each plotting command needs a crs specified in the call. Glad that’s not the case w/ Julia.
Easy. Instead of interpolating ask it to return the x,y,z and do a scatterplot.
using RemoteS, GMT
# Return the ``lon lat sst`` in a GMTdataset type
D = grid_at_sensor("AQUA_MODIS.20220310T130000.L2.SST.NRT.nc", "sst", xyz=true);
# Create a colormap (fields 5 & 6 of the global bounding box have already the sst_min/sst_max)
C = makecpt(range=(D.ds_bbox[5], D.ds_bbox[6]),);
# show it
imshow(D, C=C, zcolor=D[:,3], marker=:square, ms="1p", proj=:guess, colorbar=true, coast=true)
Thanks, this looks convenient and is very helpful. If I have 1000 polygons of various shapes, each one as vectors x,y with value v, would I use GMT.plot to put them on a map? In GMT (not GMT.jl) I use the plot command giving it an input file with the polygon coordinates and values to assign color. Can this be done in GMT.jl without first writing the polygon data to file or stdout? Doing it all by passing the data in memory?
Sure. You’ll need to create a vector of GMTdataset (see the mat2ds function) object and then pass it as first argument to the plot and use option Z (also aliased to level or levels) much like you’d use it in plain GMT.
What you are describing is very similar to choropleth maps. See this example.
@dlaklan also built a couple of very nice US examples using dataframes but he didn’t yet finish the tutorial
