Skip to content

The Shape of a GeoRegion

In this section, we go through process of extracting information on the shape of the GeoRegion.

Let us now setup the example

julia
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

We load a predefined AR6 GeoRegion for Northwest North America:

julia
geo = GeoRegion("AR6_NWN")
The GeoRegion AR6_NWN has the following properties:
    Region ID     (ID) : AR6_NWN
    Parent ID    (pID) : GLB
    Name        (name) : Northwest North America
    Bounds   (N,S,E,W) : 81.0, 50.0, -105.0, -168.0
    Rotation       (θ) : 0.0
    File Path   (path) : /home/runner/work/GeoRegions.jl/GeoRegions.jl/src/.files/AR6/AR6_NWN.json
    Centroid (geometry.centroid) : [-132.5882898173895, 63.942314590781585]
    Shape       (geometry.shape) : Vector{Point{2, Float64}}(9)

Retrieving the Bounds of the GeoRegion

The bounds of the GeoRegion are given by the .N, .S, .E and .W fields in the GeoRegion struct, that denote the north and south latitude bounds, and the east and west longitude bounds.

julia
geo.N,
geo.S,
geo.E,
geo.W
(81.0, 50.0, -105.0, -168.0)

Retrieving the Rotation

By default, there is no rotation projection passed onto the GeoRegion. However, if you have passed on a rotation projection, you can retrieve it via the field , which is given in degrees.

julia
geo.θ
0.0

Retrieving the Centroid

We can also retrieve the longitude/latitude coordinates of the centroid of the GeoRegion using the .geometry.centroid field.

julia
lonc,latc = geo.geometry.centroid
2-element Point{2, Float64} with indices SOneTo(2):
 -132.5882898173895
   63.942314590781585

Retrieving the coordinates of a GeoRegion

Using the function coordinates(), we are able to retrieve the coordinates of the vertices that define the shape of the GeoRegion. In the below example, we plot the longitude and latitude points.

julia
lon,lat = coordinates(geo)

aspect = (geo.E-geo.W+4)/(geo.N-geo.S+4)
fig = Figure()
ax = Axis(
    fig[1,1],width=750,height=750/aspect,
    limits=(geo.W-2,geo.E+2,geo.S-2,geo.N+2),
    xlabel = "Longitude / º", ylabel = "Latitude / º"
)
lines!(ax,clon,clat,color=:black,linewidth=3)
lines!(ax,lon,lat,linewidth=5)
scatter!(ax,lonc,latc,markersize=20)
resize_to_layout!(fig)
fig

This plotting works on maps where the projection is equirectangular. However, there are many different map projections, and so we can specify the number of segments n on each side of the shape using the keyword argument n in order to ensure during plotting that is plotted correctly on the projection.

julia
nlon,nlat = coordinates(geo,n=50)
length(lon), length(nlon)
(9, 401)

Number of Points

The shape is defined by 8 sides. Therefore there are by default 8 + 1 = 9 coordinate points to close the polygon. If we specify 50 segments a side, we therefore have 400 + 1 = 401 coordinate points to close the polygon.

See the API here