Rectilinear Grids for Data Extraction

The most straightforward of the RegionGrid types is the RectilinearGrid. This is the type that is used for most datasets on a rectilinear longitude/latitude grid. Examples of such datasets include:

• Level 3 products from the Global Precipitation Measurement Mission

• Final regridded products from reanalysis such as ERA5 and MERRA2

• Model output from simple climate models such as Isca and SpeedyWeather.jl

Basically, for each of these datasets, the data is given in such a way that the coordinates of the grid can be expressed via two vectors/ranges:

• A vector/range of longitudes

• A vector/range of latitudes

using GeoRegions
using RegionGrids
using CairoMakie

Creating Rectilinear Grids

A Rectilinear Grid can be created as follows:

ggrd = RegionGrid(geo,lon,lat)

where geo is a GeoRegion of interest that is found within the domain defined by the longitude and latitude grid vectors.

lon = collect(0:5:359);  nlon = length(lon)
lat = collect(-90:5:90); nlat = length(lat)
geo = GeoRegion([10,230,-50,10],[50,10,-40,50])

ggrd = RegionGrid(geo,lon,lat)

The RLinearMask Grid type has the following properties:
    Longitude Indices     (ilon) : [63, 64, 65, 66, 67, 68, 69, 70, 71, 72  …  38, 39, 40, 41, 42, 43, 44, 45, 46, 47]
    Latitude Indices      (ilat) : [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29]
    Longitude Points       (lon) : [-50, -45, -40, -35, -30, -25, -20, -15, -10, -5  …  185, 190, 195, 200, 205, 210, 215, 220, 225, 230]
    Latitude Points        (lat) : [-40, -35, -30, -25, -20, -15, -10, -5, 0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50]
    Rotated X Coordinates    (X)
    Rotated Y Coordinates    (Y)
    Rotation (°)             (θ) : 0.0
    RegionGrid Mask       (mask)
    RegionGrid Weights (weights)
    RegionGrid Size              : 57 lon points x 19 lat points
    RegionGrid Validity		  : 451 / 1083

The API for creating a Rectilinear Grid can be found here

What is in a Rectilinear Grid?

RegionGrids.RectilinearGrid Type

RectilinearGrid <: RegionGrid

A RectilinearGrid is a RegionGrid that is created based on rectilinear longitude/latitude grids. It has its own subtypes: RectGrid, TiltGrid and PolyGrid.

All RectilinearGrid types contain the following fields:

• lon - A Vector of Floats, defining the longitude vector describing the region.

• lat - A Vector of Floats, defining the latitude vector describing the region.

• ilon - A Vector of Ints, defining the indices used to extract the longitude vector from the input longitude vector.

• ilat - A Vector of Ints, defining the indices used to extract the latitude vector from the input latitude vector.

• mask - An Array of NaNs and 1s, defining the gridpoints in the RegionGrid where the data is valid.

• weights - A Vector of Floats, defining the latitude-weights of each valid point in the grid. Will be NaN if outside the bounds of the GeoRegion used to define this RectilinearGrid.

• X - A Vector of Floats, defining the X-coordinates (in meters) of each point in the "derotated" RegionGrid about the centroid for the shape of the GeoRegion.

• Y - A Vector of Floats, defining the Y-coordinates (in meters) of each point in the "derotated" RegionGrid about the centroid for the shape of the GeoRegion.

• θ - A Float storing the information on the angle (in degrees) about which the data was rotated in the anti-clockwise direction. Mathematically, it is rotation - geo.θ.

source

We see that in a RectilinearGrid type, we have the lon and lat fields that defined the longitude and latitude vectors that have been cropped to fit the GeoRegion bounds.

ggrd.lon, ggrd.lat

([-50, -45, -40, -35, -30, -25, -20, -15, -10, -5  …  185, 190, 195, 200, 205, 210, 215, 220, 225, 230], [-40, -35, -30, -25, -20, -15, -10, -5, 0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50])

An example of using Rectilinear Grids

Say we have some sample data, here randomly generated.

data = rand(nlon,nlat)

72×37 Matrix{Float64}:
 0.594031   0.370406   0.653784   0.658668   …  0.200389  0.7545     0.395881
 0.311228   0.775029   0.379001   0.467459      0.932197  0.109529   0.551655
 0.144867   0.499142   0.131413   0.107546      0.533819  0.357134   0.509336
 0.699063   0.259849   0.409234   0.348452      0.470889  0.574833   0.684293
 0.858819   0.856203   0.214188   0.989466      0.192973  0.598771   0.470687
 0.115923   0.523977   0.220766   0.514563   …  0.534718  0.219755   0.713796
 0.863104   0.730402   0.582665   0.307259      0.217736  0.362924   0.740164
 0.407251   0.328958   0.964691   0.395981      0.148941  0.253645   0.546946
 0.350294   0.818371   0.558631   0.265805      0.838154  0.056618   0.98289
 0.665256   0.83654    0.983255   0.631644      0.19524   0.199121   0.63007
 ⋮                                           ⋱            ⋮          
 0.914544   0.264683   0.0220539  0.200672      0.55459   0.994349   0.48262
 0.944744   0.517995   0.381262   0.178171      0.805901  0.0183933  0.638374
 0.826313   0.694451   0.844945   0.46647    …  0.239759  0.789882   0.154788
 0.0635768  0.654077   0.756509   0.899491      0.661919  0.698493   0.401084
 0.0483262  0.105485   0.106881   0.879114      0.583252  0.741341   0.681865
 0.0554238  0.988399   0.584018   0.0258239     0.571837  0.684024   0.818872
 0.231348   0.626622   0.865014   0.311129      0.229121  0.911919   0.565592
 0.595772   0.0834849  0.0515404  0.522331   …  0.365634  0.0191625  0.669174
 0.938593   0.857992   0.539042   0.400569      0.325006  0.733879   0.281915

We extract the valid data within the GeoRegion of interest that we defined above:

ndata = extract(data,ggrd)

57×19 Matrix{Float64}:
   0.750341  NaN          NaN          …  NaN         NaN  NaN  NaN  NaN  NaN
 NaN           0.432133   NaN             NaN         NaN  NaN  NaN  NaN  NaN
 NaN           0.175487     0.133564      NaN         NaN  NaN  NaN  NaN  NaN
 NaN           0.0891875    0.66338       NaN         NaN  NaN  NaN  NaN  NaN
 NaN           0.748393     0.0624212     NaN         NaN  NaN  NaN  NaN  NaN
 NaN           0.687155     0.488382   …  NaN         NaN  NaN  NaN  NaN  NaN
 NaN         NaN            0.90279       NaN         NaN  NaN  NaN  NaN  NaN
 NaN         NaN            0.765002      NaN         NaN  NaN  NaN  NaN  NaN
 NaN         NaN            0.0876515     NaN         NaN  NaN  NaN  NaN  NaN
 NaN         NaN            0.100495        0.584289  NaN  NaN  NaN  NaN  NaN
   ⋮                                   ⋱                     ⋮            
 NaN         NaN          NaN             NaN         NaN  NaN  NaN  NaN  NaN
 NaN         NaN          NaN             NaN         NaN  NaN  NaN  NaN  NaN
 NaN         NaN          NaN          …  NaN         NaN  NaN  NaN  NaN  NaN
 NaN         NaN          NaN             NaN         NaN  NaN  NaN  NaN  NaN
 NaN         NaN          NaN             NaN         NaN  NaN  NaN  NaN  NaN
 NaN         NaN          NaN             NaN         NaN  NaN  NaN  NaN  NaN
 NaN         NaN          NaN             NaN         NaN  NaN  NaN  NaN  NaN
 NaN         NaN          NaN          …  NaN         NaN  NaN  NaN  NaN  NaN
 NaN         NaN          NaN             NaN         NaN  NaN  NaN  NaN  NaN

And now let us visualize the results.

slon,slat = coordinates(geo) # extract the coordinates
fig = Figure()

ax1 = Axis(
    fig[1,1],width=450,height=150,
    limits=(-180,360,-90,90)
)
heatmap!(ax1,lon,lat,data,colorrange=(-1,1))
lines!(ax1,slon,slat,color=:black,linewidth=2)
lines!(ax1,slon.+360,slat,color=:black,linewidth=2,linestyle=:dash)

hidexdecorations!(ax1,ticks=false,grid=false)

ax2 = Axis(
    fig[2,1],width=450,height=150,
    limits=(-180,360,-90,90)
)
heatmap!(ax2,ggrd.lon,ggrd.lat,ndata,colorrange=(-1,1))
lines!(ax2,slon,slat,color=:black,linewidth=2)

Label(fig[3,:],"Longitude / º")
Label(fig[:,0],"Latitude / º",rotation=pi/2)

resize_to_layout!(fig)
fig