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 ofFloat
s, defining the longitude vector describing the region.lat
- A Vector ofFloat
s, defining the latitude vector describing the region.ilon
- A Vector ofInt
s, defining the indices used to extract the longitude vector from the input longitude vector.ilat
- A Vector ofInt
s, 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 ofFloat
s, 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 ofFloat
s, 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 ofFloat
s, defining the Y-coordinates (in meters) of each point in the "derotated" RegionGrid about the centroid for the shape of the GeoRegion.θ
- AFloat
storing the information on the angle (in degrees) about which the data was rotated in the anti-clockwise direction. Mathematically, it isrotation - geo.θ
.
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.976239 0.663169 0.927763 … 0.910719 0.490349 0.0775779
0.670964 0.891505 0.076615 0.108977 0.967316 0.948079
0.00942298 0.363765 0.728947 0.073005 0.299825 0.592984
0.830007 0.568142 0.162029 0.525964 0.251173 0.511037
0.951316 0.993413 0.100056 0.686316 0.013777 0.22441
0.041036 0.132841 0.133508 … 0.741237 0.582849 0.686391
0.135486 0.249522 0.464281 0.863572 0.795542 0.456205
0.47052 0.661186 0.511609 0.347628 0.826045 0.510852
0.264533 0.915567 0.960266 0.323391 0.33149 0.402225
0.032835 0.456556 0.457584 0.158352 0.800315 0.898252
⋮ ⋱ ⋮
0.746238 0.00421001 0.568895 0.102521 0.928948 0.278281
0.500542 0.422493 0.917065 0.857643 0.965416 0.708134
0.320116 0.571266 0.714261 … 0.542388 0.524515 0.534764
0.427021 0.945495 0.526565 0.839656 0.76869 0.799369
0.513393 0.553173 0.356839 0.625251 0.419471 0.920456
0.189636 0.226471 0.374043 0.228017 0.115303 0.00335935
0.441508 0.142631 0.160511 0.910355 0.419984 0.0626689
0.477382 0.558011 0.728735 … 0.341774 0.0366777 0.614395
0.458544 0.170001 0.931644 0.551886 0.504976 0.420585
We extract the valid data within the GeoRegion of interest that we defined above:
ndata = extract(data,ggrd)
57×19 Matrix{Float64}:
0.947924 NaN NaN … NaN NaN NaN NaN NaN NaN
NaN 0.975769 NaN NaN NaN NaN NaN NaN NaN
NaN 0.173636 0.0982291 NaN NaN NaN NaN NaN NaN
NaN 0.8645 0.542555 NaN NaN NaN NaN NaN NaN
NaN 0.945763 0.0691678 NaN NaN NaN NaN NaN NaN
NaN 0.785554 0.344456 … NaN NaN NaN NaN NaN NaN
NaN NaN 0.574551 NaN NaN NaN NaN NaN NaN
NaN NaN 0.78429 NaN NaN NaN NaN NaN NaN
NaN NaN 0.82661 NaN NaN NaN NaN NaN NaN
NaN NaN 0.346985 0.651541 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