GitHub

Extracting Gridded Data using RegionGrid

Suppose we have Global Data. However, we are only interested in a specific region (say, the North Central American region as defined in AR6), how do we extract data for this region?

The simple answer is, we use the extractGrid() function, which takes in a RegionGrid and a data array, and returns a new data array for the GeoRegion.

Setup

using GeoRegions
using DelimitedFiles
using CairoMakie

download("https://raw.githubusercontent.com/natgeo-wong/GeoPlottingData/main/coastline_resl.txt","coast.cst")
coast = readdlm("coast.cst",comments=true)
clon  = coast[:,1]
clat  = coast[:,2]
nothing

Let us define the GeoRegion of interest

geo  = GeoRegion("AR6_NCA")
The Polygonal Region AR6_NCA has the following properties:
    Region ID    (ID) : AR6_NCA
    Parent ID   (pID) : GLB
    Name       (name) : Northern Central America
    Bounds  (N,S,E,W) : [33.8, 16.0, -90.0, -122.5]
    Shape     (shape) : Point2{Float64}[[-90.0, 25.0], [-104.5, 16.0], [-122.5, 33.8], [-105.0, 33.8], [-90.0, 25.0]]
        (is180,is360) : (true, false)

We also define some sample test data in the global domain

lon = collect(0:360); pop!(lon); nlon = length(lon)
lat = collect(-90:90);           nlat = length(lat)
odata = randn((nlon,nlat))
360×181 Matrix{Float64}:
  0.743878    0.215347   0.797505    …   2.38487    -1.41452     -0.487161
 -0.836976   -0.522279  -0.742974       -0.0774914  -0.285505    -0.672039
  1.17802     0.407393  -0.783909        0.494963    2.11024     -2.19828
  0.0610333  -0.319797  -0.230227        0.33881     1.35005      0.656233
  1.0184     -0.102823   0.683836        1.06714    -0.00687154   0.311811
  1.07802    -0.783349   0.259099    …   0.51524     0.114212    -0.546989
 -1.23991     1.18968    0.443209       -1.8809     -0.0532506    2.22403
  0.249049    0.346461  -0.995266       -0.598945   -1.08407     -0.563606
  2.41308    -0.453135  -0.366251       -0.236916    2.03501      0.132898
  0.31055    -0.188774   0.511434       -1.89519    -0.557563     0.85635
  ⋮                                  ⋱                            ⋮
 -0.363373    1.37689   -0.0365136       0.707786    0.0865657    1.41909
 -0.847941    1.29591   -0.114976        0.383235    0.0452663   -0.736994
  0.440844    0.021256  -1.19499         2.27946    -2.28355      0.254091
  0.973469   -0.707987   1.65242         1.64485     0.483906    -0.381885
  1.5878      0.590023   0.125085    …  -3.03593    -0.305305     2.008
  1.68234    -1.34705   -1.06532        -2.41418    -0.36214      0.113718
 -0.700984    1.40048   -0.499192       -0.0660869  -0.285646    -0.730803
  0.970599    0.3441    -0.00305265     -0.149035   -0.840992    -1.10148
  0.767426   -0.495731  -0.237186       -1.30129     1.29208      0.458902

Our next step is to define the RegionGrid based on the longitude and latitude vectors and their intersection with the GeoRegion

ggrd = RegionGrid(geo,lon,lat)
The Polygonal Grid has the following properties:
    Grid Bounds       (grid) : [124, 107, 271, 239]
    Longitude Indices (ilon) : [239, 240, 241, 242, 243, 244, 245, 246, 247, 248  …  262, 263, 264, 265, 266, 267, 268, 269, 270, 271]
    Latitude Indices  (ilat) : [107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124]
    Longitude Points   (lon) : [-122.0, -121.0, -120.0, -119.0, -118.0, -117.0, -116.0, -115.0, -114.0, -113.0  …  -99.0, -98.0, -97.0, -96.0, -95.0, -94.0, -93.0, -92.0, -91.0, -90.0]
    Latitude Points    (lat) : [16.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0, 26.0, 27.0, 28.0, 29.0, 30.0, 31.0, 32.0, 33.0]
    Region Size (nlon * nlat) : 33 lon points x 18 lat points
    Region Mask (sum(mask) / (nlon * nlat)) : 281 / 594

Then we extract the data

ndata = extractGrid(odata,ggrd)
33×18 Matrix{Float64}:
 NaN  NaN  NaN  NaN  NaN  NaN         …  NaN         NaN          NaN
 NaN  NaN  NaN  NaN  NaN  NaN            NaN         NaN           -1.42362
 NaN  NaN  NaN  NaN  NaN  NaN            NaN           0.456088    -0.417757
 NaN  NaN  NaN  NaN  NaN  NaN              0.583188   -0.306318     0.48176
 NaN  NaN  NaN  NaN  NaN  NaN             -0.483981   -0.115201     0.268305
 NaN  NaN  NaN  NaN  NaN  NaN         …   -2.18194    -0.683028    -0.135998
 NaN  NaN  NaN  NaN  NaN  NaN              0.517423    0.487337     1.31885
 NaN  NaN  NaN  NaN  NaN  NaN              0.825275   -1.53006     -0.484056
 NaN  NaN  NaN  NaN  NaN  NaN             -0.307652    0.0673038   -1.92951
 NaN  NaN  NaN  NaN  NaN  NaN              0.699475   -0.542391    -1.49064
   ⋮                        ⋮         ⋱    ⋮                      
 NaN  NaN  NaN  NaN  NaN   -0.426891     NaN         NaN          NaN
 NaN  NaN  NaN  NaN  NaN   -0.948739  …  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

Let us plot the old and new data

fig = Figure()
_,_,slon,slat = coordGeoRegion(geo); slon = slon
aspect = (maximum(slon)-minimum(slon))/(maximum(slat)-minimum(slat))

ax = Axis(
    fig[1,1],width=350,height=350/aspect,
    limits=(minimum(slon)-2+360,maximum(slon)+2+360,minimum(slat)-2,maximum(slat)+2)
)
contourf!(
    ax,lon,lat,odata,
    levels=range(-1,1,length=11),extendlow=:auto,extendhigh=:auto
)
lines!(ax,clon,clat,color=:black)
lines!(ax,slon.+360,slat.+360,color=:red,linewidth=5)

ax = Axis(
    fig[1,2],width=350,height=350/aspect,
    limits=(minimum(slon)-2,maximum(slon)+2,minimum(slat)-2,maximum(slat)+2)
)
contourf!(
    ax,ggrd.lon,ggrd.lat,ndata,
    levels=range(-1,1,length=11),extendlow=:auto,extendhigh=:auto
)
lines!(ax,clon,clat,color=:black)
lines!(ax,slon,slat,color=:red,linewidth=5)

hideydecorations!(ax, ticks = false,grid=false)

resize_to_layout!(fig)
fig
Example block output

API

GeoRegions.extractGrid — Function
extractGrid(
    odata :: AbstractArray{<:Real,N},
    ggrd  :: RegionGrid
) where N <: Int -> Array{<:Real,N}

Extracts data from odata, an Array of dimension N (where N ∈ 2,3,4) that contains data of a Parent GeoRegion, into another Array of dimension N, containing only within a sub GeoRegion we are interested in.

Warning

Please ensure that the 1st dimension is longitude and 2nd dimension is latitude before proceeding. The order of the 3rd and 4th dimensions (when used), however, is not significant.

Arguments

  • odata : An array of dimension N containing gridded data in the region we are interesting in extracting from
  • ggrd : A RegionGrid containing detailed information on what to extract
source
GeoRegions.extractGrid! — Function
extractGrid!(
    odata :: AbstractArray{<:Real,N},
    ndata :: AbstractArray{<:Real,N},
    ggrd  :: RegionGrid
) where N <: Int -> nothing

Extracts data from odata, an Array of dimension N (where N ∈ 2,3,4) that contains data of a Parent GeoRegion, into ndata, another Array of dimension N, containing only within a sub GeoRegion we are interested in.

This allows for iterable in-place modification to save memory space and reduce allocations if the dimensions are fixed.

Warning

Please ensure that the 1st dimension is longitude and 2nd dimension is latitude before proceeding. The order of the 3rd and 4th dimensions (when used), however, is not significant.

Arguments

  • odata : An array of dimension N containing gridded data in the region we are interesting in extracting from
  • ndata : An array of dimension N meant as a placeholder for the extracted data to minimize allocations
  • ggrd : A RegionGrid containing detailed information on what to extract
source

This page was generated using Literate.jl.