5  Raster-vector interactions

Prerequisites

This chapter requires importing the following packages:

using GeoDataFrames, DataFrames
using Rasters, ArchGDAL # Raster I/O and operations
using Proj # activate reprojection capabilities
import GeoInterface as GI, GeometryOps as GO, LibGEOS as LG # Vector operations
import GeoFormatTypes as GFT # for CRS types
using GeoMakie, CairoMakie # plotting
WARNING: Method definition convert_arguments(MakieCore.CellGrid, Rasters.AbstractRaster{var"#s115", 2, D, A} where A where D where var"#s115") in module RastersMakieExt at /home/runner/.julia/packages/Rasters/PZMX6/ext/RastersMakieExt/plotrecipes.jl:304 overwritten at /home/runner/.julia/packages/Rasters/PZMX6/ext/RastersMakieExt/plotrecipes.jl:308.
ERROR: Method overwriting is not permitted during Module precompilation. Use `__precompile__(false)` to opt-out of precompilation.

Here, we will also set the theme in Makie.jl to ensure that all surface and heatmap plots are represented in file as bitmap images, rather than the raw data.
This is necessary only because Quarto breaks on SVGs above a certain size.
You can feel free to set this for your own work, but it is not required.

Makie.set_theme!(
    Heatmap = (; rasterize = 2),
    Surface = (; rasterize = 2),
)

It also relies on the following data files:

src_srtm = Raster("data/srtm.tif")
src_nlcd = Raster("data/nlcd.tif")
src_grain = Raster("output/grain.tif")
src_elev = Raster("output/elev.tif")
src_dem = Raster("data/dem.tif")
zion = GeoDataFrames.read("data/zion.gpkg")
zion_points = GeoDataFrames.read("data/zion_points.gpkg")
cycle_hire_osm = GeoDataFrames.read("data/cycle_hire_osm.gpkg")
us_states = GeoDataFrames.read("data/us_states.gpkg")
nz = GeoDataFrames.read("data/nz.gpkg")
src_nz_elev = Raster("data/nz_elev.tif")
โ•ญโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ•ฎ
โ”‚ 1115ร—1450 Raster{Float32,2} โ”‚
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”ดโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ dims โ”
  โ†“ X Projected{Float64} LinRange{Float64}(995456.5461976258, 2.109456546197626e6, 1115) ForwardOrdered Regular Intervals{Start},
  โ†’ Y Projected{Float64} LinRange{Float64}(6.190960927303582e6, 4.741960927303582e6, 1450) ReverseOrdered Regular Intervals{Start}
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ metadata โ”ค
  Metadata{Rasters.GDALsource} of Dict{String, Any} with 4 entries:
  "units"    => ""
  "offset"   => 0.0
  "filepath" => "data/nz_elev.tif"
  "scale"    => 1.0
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ raster โ”ค
  extent: Extent(X = (995456.5461976258, 2.110456546197626e6), Y = (4.741960927303582e6, 6.191960927303582e6))
  missingval: NaN32
  crs: PROJCS["unknown",GEOGCS["unknown",DATUM["Unknown based on GRS80 ellipsoid",SPHEROID["GRS 1980",6378137,298.257222101004,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",173],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",1600000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]
โ””โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”˜
 โ†“ โ†’          6.19096e6    6.18996e6  โ€ฆ    4.74296e6    4.74196e6
 9.95457e5  NaN          NaN             NaN          NaN
 โ‹ฎ                                    โ‹ฑ                 โ‹ฎ
 2.10946e6  NaN          NaN          โ€ฆ  NaN          NaN

5.1 Introduction

This chapter focuses on interactions between raster and vector geographic data models, both introduced in Chapter 1. It includes four main techniques:

  • Raster cropping and masking using vector objects (Section 5.2)
  • Extracting raster values using different types of vector data (Section 5.3)
  • Raster to vector conversion (Section 5.4)
  • Vector to raster conversion (Section 5.5)

These concepts are demonstrated using data from previous chapters, to understand their potential real-world applications.

5.2 Raster masking and cropping

Many geographic data projects involve integrating data from many different sources, such as remote sensing images (rasters) and administrative boundaries (vectors). Often the extent of input raster datasets is larger than the area of interest. In this case, raster masking, cropping, or both, are useful for unifying the spatial extent of input data (Figure 5.1 (b) and (c), and the following two examples, illustrate the difference between masking and cropping). Both operations reduce memory use and computational demand for subsequent analysis, and may be a necessary preprocessing step before creating attractive maps involving raster data.

We will use two layers to illustrate raster cropping:

  • The srtm.tif raster representing elevation, in meters above sea level, in south-western Utah: a Rasters.jl file connection named src_srtm (see Figure 5.1 (a))
  • The zion.gpkg vector layer representing the Zion National Park boundaries (a DataFrame named zion)

Both target and cropping objects must have the same projection. Since it is easier and more precise to reproject vector layers, compared to rasters, we use the following expression to reproject (?sec-reprojecting-vector-geometries) the vector layer zion into the coordinate reference system (CRS) of the raster src_srtm. The CRS defines how the coordinates of the geometry relate to locations on the surface of the Earth.

zion = GO.reproject(zion; target_crs = GI.crs(src_srtm))
1ร—12 DataFrame
Row geom UNIT_CODE GIS_Notes UNIT_NAME DATE_EDIT STATE REGION GNIS_ID UNIT_TYPE CREATED_BY METADATA PARKNAME
Polygonโ€ฆ String String String DateTime String String String String String String String
1 Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{CRS}}}, Nothing, WellKnownText{CRS}}(LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{CRS}}[LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{CRS}}([(-113.084, 37.1713), (-113.086, 37.1716), (-113.091, 37.1722), (-113.091, 37.1744), (-113.091, 37.178), (-113.091, 37.1817), (-113.091, 37.1853), (-113.091, 37.1889), (-113.091, 37.1926), (-113.091, 37.1962) โ€ฆ (-113.06, 37.1713), (-113.062, 37.1713), (-113.064, 37.1713), (-113.069, 37.1713), (-113.071, 37.1713), (-113.073, 37.1713), (-113.078, 37.1713), (-113.08, 37.1713), (-113.082, 37.1713), (-113.084, 37.1713)], nothing, WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], nothing, WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]")) ZION Lands - http://landsnet.nps.gov/tractsnet/documents/ZION/Metadata/zion_metadata.xml Zion National Park 2017-06-22T00:00:00 UT IM 1455157 National Park Lands https://irma.nps.gov/App/Reference/Profile/2181118#Zion National Monument Zion

To mask the image, i.e., convert all pixels which do not intersect with the zion polygon to missing, we use the Rasters.mask function.
mask supports any geometry, vector of geometries, feature collection, or table with a geometry column!

The tabset below shows all the different ways to mask a raster. Weโ€™ll go forward with the approach of using a DataFrame.

out_image_mask = Rasters.mask(src_srtm; with = zion)
โ•ญโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ•ฎ
โ”‚ 465ร—457 Raster{UInt16,2} โ”‚
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”ดโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ dims โ”
  โ†“ X Projected{Float64} LinRange{Float64}(-113.23958321278403, -112.85291654614313, 465) ForwardOrdered Regular Intervals{Start},
  โ†’ Y Projected{Float64} LinRange{Float64}(37.51208342983253, 37.13208342985786, 457) ReverseOrdered Regular Intervals{Start}
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ metadata โ”ค
  Metadata{Rasters.GDALsource} of Dict{String, Any} with 4 entries:
  "units"    => ""
  "offset"   => 0.0
  "filepath" => "data/srtm.tif"
  "scale"    => 1.0
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ raster โ”ค
  extent: Extent(X = (-113.23958321278403, -112.85208321280986), Y = (37.13208342985786, 37.512916763165805))
  missingval: 0xffff
  crs: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
โ””โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”˜
    โ†“ โ†’        37.5121      37.5113      37.5104  โ€ฆ      37.1329      37.1321
 -113.24   0xffff       0xffff       0xffff          0xffff       0xffff
    โ‹ฎ                                             โ‹ฑ               
 -112.853  0xffff       0xffff       0xffff       โ€ฆ  0xffff       0xffff
masker = zion.geom[1]
GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}([(-113.08438684283465, 37.17130388117393), (-113.0859893948904, 37.17162779353809), (-113.09132351665151, 37.17215207288513), (-113.09132530000656, 37.17440980000173), (-113.09132830000658, 37.17803750000173), (-113.09133150000655, 37.18166840000168), (-113.09133470000658, 37.18529940000173), (-113.09134070000657, 37.188925000001696), (-113.09134660000655, 37.19255060000173), (-113.09135260000657, 37.19617620000171)  โ€ฆ  (-113.05962130000604, 37.17131270000172), (-113.06185720000607, 37.171310800001734), (-113.06415140000612, 37.17131020000169), (-113.06868130000619, 37.17130650000175), (-113.07101920000618, 37.17130460000174), (-113.07321130000625, 37.17130450000173), (-113.07773970000633, 37.17130420000172), (-113.08035690000636, 37.17130410000175), (-113.0822681000064, 37.171304000001705), (-113.08438684283465, 37.17130388117393)], nothing, WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], nothing, WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))
Rasters.mask(src_srtm; with = masker)
โ•ญโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ•ฎ
โ”‚ 465ร—457 Raster{UInt16,2} โ”‚
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”ดโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ dims โ”
  โ†“ X Projected{Float64} LinRange{Float64}(-113.23958321278403, -112.85291654614313, 465) ForwardOrdered Regular Intervals{Start},
  โ†’ Y Projected{Float64} LinRange{Float64}(37.51208342983253, 37.13208342985786, 457) ReverseOrdered Regular Intervals{Start}
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ metadata โ”ค
  Metadata{Rasters.GDALsource} of Dict{String, Any} with 4 entries:
  "units"    => ""
  "offset"   => 0.0
  "filepath" => "data/srtm.tif"
  "scale"    => 1.0
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ raster โ”ค
  extent: Extent(X = (-113.23958321278403, -112.85208321280986), Y = (37.13208342985786, 37.512916763165805))
  missingval: 0xffff
  crs: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
โ””โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”˜
    โ†“ โ†’        37.5121      37.5113      37.5104  โ€ฆ      37.1329      37.1321
 -113.24   0xffff       0xffff       0xffff          0xffff       0xffff
    โ‹ฎ                                             โ‹ฑ               
 -112.853  0xffff       0xffff       0xffff       โ€ฆ  0xffff       0xffff
masker = zion.geom
1-element Vector{GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}}:
 GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}([(-113.08438684283465, 37.17130388117393), (-113.0859893948904, 37.17162779353809), (-113.09132351665151, 37.17215207288513), (-113.09132530000656, 37.17440980000173), (-113.09132830000658, 37.17803750000173), (-113.09133150000655, 37.18166840000168), (-113.09133470000658, 37.18529940000173), (-113.09134070000657, 37.188925000001696), (-113.09134660000655, 37.19255060000173), (-113.09135260000657, 37.19617620000171)  โ€ฆ  (-113.05962130000604, 37.17131270000172), (-113.06185720000607, 37.171310800001734), (-113.06415140000612, 37.17131020000169), (-113.06868130000619, 37.17130650000175), (-113.07101920000618, 37.17130460000174), (-113.07321130000625, 37.17130450000173), (-113.07773970000633, 37.17130420000172), (-113.08035690000636, 37.17130410000175), (-113.0822681000064, 37.171304000001705), (-113.08438684283465, 37.17130388117393)], nothing, WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], nothing, WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))
Rasters.mask(src_srtm; with = masker)
โ•ญโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ•ฎ
โ”‚ 465ร—457 Raster{UInt16,2} โ”‚
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”ดโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ dims โ”
  โ†“ X Projected{Float64} LinRange{Float64}(-113.23958321278403, -112.85291654614313, 465) ForwardOrdered Regular Intervals{Start},
  โ†’ Y Projected{Float64} LinRange{Float64}(37.51208342983253, 37.13208342985786, 457) ReverseOrdered Regular Intervals{Start}
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ metadata โ”ค
  Metadata{Rasters.GDALsource} of Dict{String, Any} with 4 entries:
  "units"    => ""
  "offset"   => 0.0
  "filepath" => "data/srtm.tif"
  "scale"    => 1.0
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ raster โ”ค
  extent: Extent(X = (-113.23958321278403, -112.85208321280986), Y = (37.13208342985786, 37.512916763165805))
  missingval: 0xffff
  crs: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
โ””โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”˜
    โ†“ โ†’        37.5121      37.5113      37.5104  โ€ฆ      37.1329      37.1321
 -113.24   0xffff       0xffff       0xffff          0xffff       0xffff
    โ‹ฎ                                             โ‹ฑ               
 -112.853  0xffff       0xffff       0xffff       โ€ฆ  0xffff       0xffff
Note

Note that since Julia has a native missing/NODATA value type, we donโ€™t need to specify a NODATA value for the mask function.

However, it can sometimes be useful and more efficient to specify a sentinel value which Rasters treats as missing.

You can do this by specifying the missingval keyword argument, like so:

out_image_mask = Rasters.mask(src_srtm; with = zion, missingval = 9999)

We can write this masked raster to file with Rasters.write:

Rasters.write("output/srtm_masked.tif", out_image_mask; force = true
)
"output/srtm_masked.tif"

In Rasters.jl, cropping and masking are distinct operations. Cropping, which reduces the raster extent to the extent of the vector layer, is accomplished with the crop function.

Here, we simply pass the zion feature table to the to keyword argument, which indicates what to crop the raster โ€œtoโ€. We also set the touches keyword argument to true, to specify that pixels that partially overlap with the vector layer are included in the output.

out_image_crop = Rasters.crop(src_srtm; to = zion, touches = true)
โ•ญโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ•ฎ
โ”‚ 439ร—436 Raster{UInt16,2} โ”‚
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”ดโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ dims โ”
  โ†“ X Projected{Float64} LinRange{Float64}(-113.22874987945141, -112.86374987947575, 439) ForwardOrdered Regular Intervals{Start},
  โ†’ Y Projected{Float64} LinRange{Float64}(37.50375009649975, 37.14125009652391, 436) ReverseOrdered Regular Intervals{Start}
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ metadata โ”ค
  Metadata{Rasters.GDALsource} of Dict{String, Any} with 4 entries:
  "units"    => ""
  "offset"   => 0.0
  "filepath" => "data/srtm.tif"
  "scale"    => 1.0
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ raster โ”ค
  extent: Extent(X = (-113.22874987945141, -112.86291654614247), Y = (37.14125009652391, 37.504583429833026))
  missingval: 0xffff
  crs: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
โ””โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”˜
    โ†“ โ†’        37.5038      37.5029      37.5021  โ€ฆ      37.1421      37.1413
 -113.229  0x0668       0x065f       0x0659          0x05a0       0x05df
    โ‹ฎ                                             โ‹ฑ       โ‹ฎ       
 -112.864  0x0a00       0x0a1a       0x0a24          0x06da       0x06d4

You can also assemble an extent manually, using Extents.Extent, or extract one using GI.extent.

We can crop our masked raster as well:

out_image_mask_crop = Rasters.crop(out_image_mask; to = zion, touches = true)
โ•ญโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ•ฎ
โ”‚ 439ร—436 Raster{UInt16,2} โ”‚
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”ดโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ dims โ”
  โ†“ X Projected{Float64} LinRange{Float64}(-113.22874987945141, -112.86374987947575, 439) ForwardOrdered Regular Intervals{Start},
  โ†’ Y Projected{Float64} LinRange{Float64}(37.50375009649975, 37.14125009652391, 436) ReverseOrdered Regular Intervals{Start}
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ metadata โ”ค
  Metadata{Rasters.GDALsource} of Dict{String, Any} with 4 entries:
  "units"    => ""
  "offset"   => 0.0
  "filepath" => "data/srtm.tif"
  "scale"    => 1.0
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ raster โ”ค
  extent: Extent(X = (-113.22874987945141, -112.86291654614247), Y = (37.14125009652391, 37.504583429833026))
  missingval: 0xffff
  crs: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
โ””โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”˜
    โ†“ โ†’        37.5038      37.5029      37.5021  โ€ฆ      37.1421      37.1413
 -113.229  0xffff       0xffff       0xffff          0xffff       0xffff
    โ‹ฎ                                             โ‹ฑ       โ‹ฎ       
 -112.864  0xffff       0xffff       0xffff          0xffff       0xffff

and we write it to file using Rasters.write:

Rasters.write("output/srtm_masked_cropped.tif", out_image_mask_crop; force = true)
"output/srtm_masked_cropped.tif"

Figure 5.1 shows the original raster, and the three masking and/or cropping results.

fig = Figure(size = (600, 600))

ax1 = Axis(fig[1, 1]; title = "Original")
plot!(ax1, src_srtm)
poly!(ax1, zion.geom; color = :transparent, strokecolor = :black, strokewidth = 0.75)

ax2 = Axis(fig[1, 2]; title = "Masked")
plot!(ax2, out_image_mask)
poly!(ax2, zion.geom; color = :transparent, strokecolor = :black, strokewidth = 0.75)

ax3 = Axis(fig[2, 1]; title = "Cropped")
plot!(ax3, out_image_crop)
poly!(ax3, zion.geom; color = :transparent, strokecolor = :black, strokewidth = 0.75)

ax4 = Axis(fig[2, 2]; title = "Masked+Cropped")
plot!(ax4, out_image_mask_crop)
poly!(ax4, zion.geom; color = :transparent, strokecolor = :black, strokewidth = 0.75)

display(fig)
(a) Raster masking and cropping
CairoMakie.Screen{PDF}
(b)
Figure 5.1

5.3 Raster extraction

Raster extraction is the process of identifying and returning the values associated with a โ€˜targetโ€™ raster at specific locations, based on a (typically vector) geographic โ€˜selectorโ€™ object. The reverse of raster extractionโ€”assigning raster cell values based on vector objectsโ€”is rasterization, described in Section 5.4.

Rasters.jl provides modular raster extraction and statistics functions, and we use this package in the following examples.

5.3.1 Extraction to points

The simplest type of raster extraction is getting the values of raster cells at specific points. To demonstrate extraction to points, we will use zion_points, which contains a sample of 30 locations within the Zion National Park (Figure 5.2).

fig, ax, plt = plot(src_srtm)
scatter!(ax, zion_points.geom, color=:black, strokecolor=:white, strokewidth = 1);
display(fig)
(a) 30 point locations within the Zion National Park, with elevation in the background
CairoMakie.Screen{PDF}
(b)
Figure 5.2

The following expression extracts elevation values from srtm.tif according to zion_points, using Rasters.extract.

The API here is not great, can we do better? It currently returns a vector of named tuples, which is not very convenient.

One thought is to use a Tables.jl materializer to convert the result if possible. I understand the desire to return the geometry values. But there must be a better way than this.

result1 = DataFrame(Rasters.extract(src_srtm, zion_points; geometry = false))
30ร—1 DataFrame
5 rows omitted
Row
UInt16?
1 1802
2 2433
3 1886
4 1370
5 1452
6 1635
7 1380
8 2032
9 1830
10 1860
11 1440
12 2145
13 1942
โ‹ฎ โ‹ฎ
19 1758
20 1424
21 2159
22 1809
23 1826
24 1550
25 1799
26 2102
27 2118
28 1372
29 1905
30 1574

The first argument is the raster from which to extract values, and the second is the vector object (or collection of objects) according to which to extract the values.

Note

Rasters.jl does not yet support interpolation in extraction, so the values extracted are the values of the nearest cell center.
This corresponds to interpolate='nearest' in the Python rasterstats package.

Either way, the resulting object is a vector of raster values, corresponding to zion_points. For example, here are the elevations of the first five points.

result1[1:5, ""]
5-element Vector{Union{Missing, UInt16}}:
 0x070a
 0x0981
 0x075e
 0x055a
 0x05ac

To get a DataFrame with the original points geometries (and other attributes, if any), as well as the extracted raster values, we can assign the extraction result into a new column.

zion_points[!, "elev1"] = result1[!, ""]
zion_points
30ร—2 DataFrame
5 rows omitted
Row geom elev1
IGeometrโ€ฆ UInt16?
1 Geometry: wkbPoint 1802
2 Geometry: wkbPoint 2433
3 Geometry: wkbPoint 1886
4 Geometry: wkbPoint 1370
5 Geometry: wkbPoint 1452
6 Geometry: wkbPoint 1635
7 Geometry: wkbPoint 1380
8 Geometry: wkbPoint 2032
9 Geometry: wkbPoint 1830
10 Geometry: wkbPoint 1860
11 Geometry: wkbPoint 1440
12 Geometry: wkbPoint 2145
13 Geometry: wkbPoint 1942
โ‹ฎ โ‹ฎ โ‹ฎ
19 Geometry: wkbPoint 1758
20 Geometry: wkbPoint 1424
21 Geometry: wkbPoint 2159
22 Geometry: wkbPoint 1809
23 Geometry: wkbPoint 1826
24 Geometry: wkbPoint 1550
25 Geometry: wkbPoint 1799
26 Geometry: wkbPoint 2102
27 Geometry: wkbPoint 2118
28 Geometry: wkbPoint 1372
29 Geometry: wkbPoint 1905
30 Geometry: wkbPoint 1574

You can read from a single band by selecting the band in the Raster. TODO finish this text

5.3.2 Extraction to lines

Raster extraction is also applicable with line selectors. The typical line extraction algorithm is to extract one value for each raster cell touched by a line. However, this particular approach is not recommended to obtain values along the transects, as it is hard to get the correct distance between each pair of extracted raster values.

For line extraction, a better approach is to split the line into many points (at equal distances along the line) and then extract the values for these points using the โ€œextraction to pointsโ€ technique (Section 5.3.1). To demonstrate this, the code below creates (see ?sec-vector-data for recap) zion_transect, a straight line going from northwest to southeast of the Zion National Park.

coords = [[-113.2, 37.45], [-112.9, 37.2]]
zion_transect = GI.LineString(coords)
GeoInterface.Wrappers.LineString{false, false, Vector{Vector{Float64}}, Nothing, Nothing}([[-113.2, 37.45], [-112.9, 37.2]], nothing, nothing)

The utility of extracting heights from a linear selector is illustrated by imagining that you are planning a hike. The method demonstrated below provides an โ€˜elevation profileโ€™ of the route (the line does not need to be straight), useful for estimating how long it will take by determining the cumulative elevation gain of your journey.

First, we need to create a layer consisting of points along our line (zion_transect), at specified intervals (e.g., 250). To do that, we need to transform the line into a projected CRS (so that we work with true distances, in \(m\)), such as Universal Transverse Mercator.

zion_transect_utm = GO.reproject(zion_transect; target_crs = GFT.EPSG(32612), source_crs = GFT.EPSG(4326))
GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, EPSG{1}}([(305399.67208180577, 4.147066650206682e6), (331380.8917453843, 4.1187500947884847e6)], nothing, EPSG{1}((32612,)))

The printout of the new geometry shows this is still a straight line between two points, only with coordinates in a projected CRS.

Iโ€™ve chosen to differ from the Python treatment here - instead of selecting some number of points along the line explicitly, I will segmentize the line and extract the points. This is less precise, but we donโ€™t have the API to do arclength interpolation in GeometryOps yet. Hopefully this will be added soon.

cf. https://github.com/JuliaGeo/GeometryOps.jl/issues/210

Here, we interpolate points along the line using GO.segmentize. This operation is sometimes called line densification.

We first compute the length of the line, and then use this to segmentize the line into approximately 250 points.

_centroid, linelen = GO.centroid_and_length(zion_transect_utm)
((318390.28191359504, 4.1329083724975833e6), 38429.82026966749)

Now that we have the length of the line, we can choose a distance such that we get around 250 points along the line.

zion_transect_line = GO.segmentize(zion_transect_utm; max_distance = linelen / 250)
GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, EPSG{1}}([(305399.67208180577, 4.147066650206682e6), (305503.59696046007, 4.146953383985009e6), (305607.5218391144, 4.1468401177633363e6), (305711.44671776873, 4.1467268515416635e6), (305815.37159642304, 4.1466135853199908e6), (305919.29647507734, 4.146500319098318e6), (306023.22135373164, 4.1463870528766452e6), (306127.14623238595, 4.1462737866549725e6), (306231.0711110403, 4.1461605204332997e6), (306334.9959896946, 4.1460472542116265e6)  โ€ฆ  (330445.5678374955, 4.11976949078354e6), (330549.4927161498, 4.119656224561867e6), (330653.41759480414, 4.119542958340194e6), (330757.34247345844, 4.1194296921185213e6), (330861.26735211274, 4.1193164258968486e6), (330965.19223076705, 4.119203159675176e6), (331069.11710942135, 4.119089893453503e6), (331173.0419880757, 4.1189766272318303e6), (331276.96686673, 4.1188633610101575e6), (331380.8917453843, 4.1187500947884847e6)], nothing, EPSG{1}((32612,)))

This gives us a collection of 251 points along the line. We can extract the points that define the line segments by using GI.getpoint on the line, and then reproject the points to the CRS of the raster.

line_points = GI.getpoint(zion_transect_line)
zion_transect_pnt = GO.reproject(line_points; target_crs = GI.crs(src_srtm), source_crs = GI.crs(zion_transect_line))
251-element Vector{Tuple{Float64, Float64}}:
 (-113.2, 37.45)
 (-113.1987960226292, 37.449001661297636)
 (-113.19759207733796, 37.448003309143594)
 (-113.1963881641247, 37.44700494353919)
 (-113.1951842829878, 37.44600656448574)
 (-113.19398043392566, 37.44500817198453)
 (-113.19277661693668, 37.444009766036885)
 (-113.19157283201925, 37.4430113466441)
 (-113.19036907917175, 37.4420129138075)
 (-113.18916535839257, 37.44101446752837)
 โ‹ฎ
 (-112.90956920045984, 37.20801281544639)
 (-112.90837293949568, 37.20701125948384)
 (-112.90717671022236, 37.20600969038492)
 (-112.90598051263827, 37.205008108150906)
 (-112.90478434674178, 37.2040065127831)
 (-112.90358821253133, 37.20300490428282)
 (-112.9023921100053, 37.20200328265134)
 (-112.90119603916204, 37.201001647889974)
 (-112.90000000000002, 37.20000000000002)

Finally, we extract the elevation values for each point in our transect and combine the information with zion_transect_pnt (after โ€œpromotingโ€ it to a DataFrame, to accommodate extra attributes), using the point extraction method shown earlier (Section 5.3.1). We also attach the distance of each point along the line, to be used to plot an elevation profile.

zion_transect_pnt = DataFrame(geometry = zion_transect_pnt)

result = Rasters.extract(src_srtm, zion_transect_pnt; geometry = false)
# `Rasters.extract` returns a vector of tuples, so we take the first element of each tuple
result = first.(result)
# Add elevation data to a new data frame column, named `elevation`
zion_transect_pnt[!, "elevation"] = result
251-element Vector{UInt16}:
 0x07d1
 0x07f1
 0x07e4
 0x07c5
 0x076e
 0x074d
 0x06fe
 0x06e7
 0x06d3
 0x06ba
      โ‹ฎ
 0x0720
 0x079f
 0x0769
 0x0737
 0x072d
 0x0737
 0x0729
 0x070a
 0x06e3

We also want to visualize an elevation profile along this line, so we compute the distances along the line manually. Using GO.distance, we can get the distance between successive points along the line. Then, we add 0.0 to the beginning of the array of distances, and sum along that array to get the cumulative distance along the line.

# Compute distances between successive points along the line
distances_between_points = GO.distance.(zion_transect_pnt.geometry[1:end-1], zion_transect_pnt.geometry[2:end])
# Compute cumulative distances along the line
cumulative_distances = cumsum([0.0; distances_between_points])
# Assign the distances to a new column in the data frame
zion_transect_pnt[!, "distance"] = cumulative_distances
zion_transect_pnt
251ร—3 DataFrame
226 rows omitted
Row geometry elevation distance
Tupleโ€ฆ UInt16 Float64
1 (-113.2, 37.45) 2001 0.0
2 (-113.199, 37.449) 2033 0.00156405
3 (-113.198, 37.448) 2020 0.00312808
4 (-113.196, 37.447) 1989 0.00469209
5 (-113.195, 37.446) 1902 0.00625609
6 (-113.194, 37.445) 1869 0.00782007
7 (-113.193, 37.444) 1790 0.00938404
8 (-113.192, 37.443) 1767 0.010948
9 (-113.19, 37.442) 1747 0.0125119
10 (-113.189, 37.441) 1722 0.0140758
11 (-113.188, 37.44) 1770 0.0156397
12 (-113.187, 37.439) 1888 0.0172036
13 (-113.186, 37.438) 1850 0.0187675
โ‹ฎ โ‹ฎ โ‹ฎ โ‹ฎ
240 (-112.913, 37.211) 1735 0.373351
241 (-112.912, 37.21) 1767 0.374911
242 (-112.911, 37.209) 1778 0.376472
243 (-112.91, 37.208) 1824 0.378032
244 (-112.908, 37.207) 1951 0.379592
245 (-112.907, 37.206) 1897 0.381152
246 (-112.906, 37.205) 1847 0.382712
247 (-112.905, 37.204) 1837 0.384273
248 (-112.904, 37.203) 1847 0.385833
249 (-112.902, 37.202) 1833 0.387393
250 (-112.901, 37.201) 1802 0.388953
251 (-112.9, 37.2) 1763 0.390513

The information in zion_transect_pnt, namely the "dist" and "elev" attributes, can now be used to draw an elevation profile, as illustrated in Figure 5.3.

# Raster and a line transect
fig, ax, plt = plot(src_srtm)
lines!(ax, zion_transect; color = :black)
poly!(ax, zion.geom; color = :transparent, strokecolor = :white, strokewidth = 0.75)
display(fig)
# Elevation profile
fig, ax, plt = lines(
    zion_transect_pnt.distance, 
    zion_transect_pnt.elevation;
    axis = (;
        xlabel = "Distance (m)",
        ylabel = "Elevation (m)",
    )
)
display(fig)
(a) Raster and a line transect
(b) Extracted elevation profile
CairoMakie.Screen{PDF}
(c)
Figure 5.3: Extracting a raster values profile to line

5.3.3 Extraction to polygons

The final type of geographic vector object that can be used for raster extraction is polygons.
Like lines, polygons tend to return many raster values per vector geometry. For continuous rasters (Figure 5.4 (a)), we typically want to generate summary statistics for raster values per polygon, for example to characterize a single region or to compare many regions. The generation of raster summary statistics, by polygons, is demonstrated in the code below using Rasters.zonal, which creates a list of summary statistics for โ€œzonesโ€ defined by geometry. In this case, a vector of length 1 is returned, since there is just one polygon in the DataFrame.

using Statistics # for `mean`
rmean = Rasters.zonal(mean, src_srtm; of = zion)
rmin = Rasters.zonal(minimum, src_srtm; of = zion)
rmax = Rasters.zonal(maximum, src_srtm; of = zion)
result = (rmean, rmin, rmax)
(Union{Missing, Float64}[1818.211830154405], Union{Missing, UInt16}[0x0462], Union{Missing, UInt16}[0x0a65])

Rasters.zonal accepts any function that works on iterables (collections of elements that can be looped over) and returns a single value, like mean, minimum, maximum, std, median, mode, sum, prod, etc.

# This is about the worst possible implementation of a mean function,
# please don't use it in real life!
function my_mean(iterable)
    result = 0.0
    count = 0
    for value in iterable
        result += value
        count += 1
    end
    result /= count
    return result
end

Rasters.zonal(my_mean, src_srtm; of = zion)
1-element Vector{Union{Missing, Float64}}:
 1818.211830154405

Sometimes, in order to be most efficient with raster access, we might want to pass multiple functions to a single Rasters.zonal call. We can do this by passing a function that returns a tuple of values to Rasters.zonal.

For example, consider:

using Statistics # for `mean`
Rasters.zonal(x -> (mean(x), minimum(x), maximum(x)), src_srtm; of = zion)
1-element Vector{Union{Missing, Tuple{Float64, UInt16, UInt16}}}:
 (1818.211830154405, 0x0462, 0x0a65)

We then transform the result to a DataFrame, which makes it easier to track and handle geometry attributes:

DataFrame(mean = rmean, min = rmin, max = rmax)
1ร—3 DataFrame
Row mean min max
Float64? UInt16? UInt16?
1 1818.21 1122 2661

Because there is only one polygon in the example, single-element vectors are returned. However, if zion was composed of more than one polygon, we would accordingly get more elements in the returned vectors. The result provides useful summaries, for example that the maximum height in the park is 2661 \(m\) above sea level.

To count occurrences of categorical raster values within polygons (Figure 5.4 (b)), we can use masking (Section 5.2) combined with StatsBase.countmap, as follows.

out_image = Rasters.mask(src_nlcd; with = GO.reproject(zion; target_crs = GI.crs(src_nlcd)))
using StatsBase
counts = StatsBase.countmap(out_image)
Dict{UInt8, Int64} with 8 entries:
  0x05 => 203701
  0x04 => 298299
  0x06 => 235
  0x07 => 62
  0x02 => 4205
  0xff => 852741
  0x08 => 679
  0x03 => 98285

According to the result, for example, the value 2 (โ€œDevelopedโ€ class) appears in 4205 pixels within the Zion polygon.

Figure 5.4 illustrates the two types of raster extraction to polygons described above.

# Continuous raster
fig, ax, plt = plot(src_srtm)
poly!(ax, zion.geom; color = :transparent, strokecolor = :black, strokewidth = 0.75)
display(fig)
# Categorical raster
fig, ax, plt = plot(src_nlcd; colormap = cgrad(:Set3; categorical = true), source = GI.crs(src_nlcd), axis = (; type = GeoAxis, dest = GI.crs(src_nlcd)))
poly!(ax, zion.geom; source = GI.crs(zion.geom[1]), color = :transparent, strokecolor = :black, strokewidth = 0.75)
cm = Colorbar(fig[1, 2], plt)
ax.xgridvisible = false
ax.ygridvisible = false
display(fig) # TODO: make GeoMakie better on small extents
(a) Continuous raster
(b) Categorical raster
CairoMakie.Screen{PDF}
(c)
Figure 5.4: Sample data used for continuous and categorical raster extraction to a polygon

5.4 Rasterization

Rasterization is the conversion of vector objects into their representation in raster objects. Usually, the output raster is used for quantitative analysis (e.g., analysis of terrain) or modeling. As we saw in Chapter 1, the raster data model has some characteristics that make it conducive to certain methods. Furthermore, the process of rasterization can help simplify datasets because the resulting values all have the same spatial resolution: rasterization can be seen as a special type of geographic data aggregation.

Rasters.jl provides the Rasters.rasterize function for rasterizing vector data. To make this happen, we need to have some definition for a โ€œtemplateโ€ grid, i.e., the โ€œtemplateโ€ raster defining the extent, resolution and CRS of the output. We can also pass a pre-existing raster, in which case Rasters uses the existing grid definition.

As for the vector geometries and their associated values, the Rasters.rasterize function can take input in multiple ways: - ; data::FeatureCollection, fill::Symbol means that the values of the column passed to fill will be used to fill the raster. - ; data::Vector{Geometry}, fill::Vector means that the values passed to fill will be associated with each geometry in data and used to fill the raster. - ; data::Any, fill::Function calls fill with the current value

Furthermore, we define how to handle multiple values burned into the same pixel, in the first argument called reducer.
By default, this is last, meaning that the last polygon to be rasterized takes precedence.
However, we can pass any function that takes in an iterable and returns a single value, like mean, minimum, maximum, std, median, mode, sum, prod, etc.
Many of these may not make sense to use but they are useful to know of.

Furthermore, we define how to deal with multiple values burned into the same pixel, in the first argument called reducer. By default, this is last, meaning that the last polygon to be rasterized takes precedence. However, we can pass any function that takes in an iterable and returns a single value, like mean, minimum, maximum, std, median, mode, sum, prod, etc. Many of these may not make sense to use but they are useful to know of.

Finally, we can set the fill value, which is the value that โ€œunaffectedโ€ pixels get, with fill=0 being the default.

How the Rasters.rasterize function works with all of these various parameters will be made clear in the next examples.

The spatial resolution (area of each โ€œpixelโ€) of the โ€œtemplateโ€ raster has a major impact on the results: if it is too low (cell size is too large), the result may miss the full geographic variability of the vector data; if it is too high, computational times may be excessive. There are no simple rules to follow when deciding an appropriate geographic resolution, which is heavily dependent on the intended use of the results. Often the target resolution is imposed on the user, for example when the output of rasterization needs to be aligned to an existing raster.

Depending on the input data, rasterization typically takes one of two forms which we demonstrate next:

  • in point rasterization (Section 5.4.1), we typically choose how to treat multiple points: either to summarize presence/absence, point count, or summed attribute values (Figure 5.5)
  • in line and polygon rasterization (Section 5.4.2), there are typically no such โ€œoverlapsโ€ and we simply โ€œburnโ€ attribute values, or fixed values, into pixels coinciding with the given geometries (Figure 5.6)

5.4.1 Rasterizing points

To demonstrate point rasterization, we will prepare a โ€œtemplateโ€ raster that has the same extent and CRS as the input vector data cycle_hire_osm_projected (a dataset on cycle hire points in London, illustrated in Figure 5.5 (a)) and a spatial resolution of 1000 \(m\). To do that, we first take our point layer and transform it to a projected CRS.

cycle_hire_osm_projected = GO.reproject(cycle_hire_osm; target_crs = GFT.EPSG(27700))
540ร—6 DataFrame
515 rows omitted
Row geom osm_id name capacity cyclestreets_id description
Tupleโ€ฆ String String? Float64? String? String?
1 (5.32354e5, 1.82858e5) 108539 Windsor Terrace 14.0 missing missing
2 (5.29848e5, 1.83337e5) 598093293 Pancras Road, King's Cross missing missing missing
3 (5.30636e5, 182609.0) 772536185 Clerkenwell, Ampton Street 11.0 missing missing
4 (5.3254e5, 1.82496e5) 772541878 missing missing missing missing
5 (5.30432e5, 1.82907e5) 781506147 missing missing missing missing
6 (5.31639e5, 182485.0) 783824668 Finsbury Library, EC1 missing missing missing
7 (5.31047e5, 1.82548e5) 783824677 Margery Street, WC1 missing missing missing
8 (5.32621e5, 1.81945e5) 787313764 Chiswell Street, EC1 20.0 missing missing
9 (5.27892e5, 1.81374e5) 814136668 George Place Mews 6.0 missing missing
10 (5.30399e5, 1.81206e5) 820945916 Drury Lane 17.0 missing missing
11 (5.26267e5, 1.78812e5) 823375879 Gloucester Road Station missing missing missing
12 (5.32358e5, 1.79686e5) 826142422 The Borough, Borough High Street 11.0 missing missing
13 (5.32179e5, 1.81593e5) 826142862 missing missing missing missing
โ‹ฎ โ‹ฎ โ‹ฎ โ‹ฎ โ‹ฎ โ‹ฎ โ‹ฎ
529 (5.38051e5, 1.84316e5) 4692548549 Stratford Aquatics Centre missing missing missing
530 (5.38235e5, 1.84568e5) 4692548550 missing missing missing missing
531 (5.37978e5, 1.84043e5) 4692548551 missing missing missing missing
532 (5.35957e5, 1.82213e5) 4692553573 missing missing missing missing
533 (5.29512e5, 1.77578e5) 4715048896 Riverlight North missing missing missing
534 (535867.0, 1.81696e5) 4872955240 missing missing missing missing
535 (535870.0, 1.81699e5) 4872955241 missing missing missing missing
536 (5.35854e5, 1.81701e5) 4872955242 missing missing missing missing
537 (5.35867e5, 1.81716e5) 4872955243 missing missing missing missing
538 (5.32531e5, 1.78805e5) 5121513755 missing 5.0 missing missing
539 (5.38723e5, 1.80836e5) 5188912370 missing missing missing missing
540 (5.38413e5, 1.80774e5) 5188912371 missing missing missing missing

We can then use the Rasters.rasterize function to rasterize the points.

This isnโ€™t a great way to get an extent, but needs must. Currently we get the extent of cycle_hire_osm_projected by GI.extent(GI.LineString(cycle_hire_osm_projected.geom)).

Track https://github.com/geocompx/geocompjl/issues/5 to see if thereโ€™s a better way to get an extent from a vector of geometries.

As mentioned above, point rasterization can be a very flexible operation: the results depend not only on the nature of the template raster, but also on the the pixel โ€œactivationโ€ method, namely the way we deal with multiple points matching the same pixel.

To illustrate this flexibility, we will try three different approaches to point rasterization (Figure 5.5 (b)-(d)). First, we create a raster representing the presence or absence of cycle hire points (known as presence/absence rasters). In this case, we transfer the value of 1 to all pixels where at least one point falls in. In the Rasters.jl framework, we use the Rasters.rasterize function, as described above. In this first example, we want to write the value 1 where the points are present, and 0 otherwise.

ch_raster1 = Rasters.rasterize(
    last,                     # reducer
    cycle_hire_osm_projected; # data
    fill = 1, 
    size = (1000, 1000) # specify size in "pixels"
)
โ•ญโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ•ฎ
โ”‚ 1000ร—1000 Raster{Union{Missing, Int64},2} last โ”‚
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”ดโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ dims โ”
  โ†“ X Projected{Float64} 523038.61452275474:15.684723604717:538707.653403867 ForwardOrdered Regular Intervals{Start},
  โ†’ Y Projected{Float64} 174934.65727249475:10.03675127048517:184961.37179170942 ForwardOrdered Regular Intervals{Start}
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ metadata โ”ค
  Metadata of Dict{Any, Any}()
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ raster โ”ค
  extent: Extent(X = (523038.61452275474, 538723.3381274717), Y = (174934.65727249475, 184971.40854297992))
  missingval: missing
โ””โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”˜
      โ†“ โ†’        1.74935e5  1.74945e5  โ€ฆ  1.84941e5  1.84951e5  1.84961e5
      5.23039e5   missing    missing       missing    missing    missing
      5.23054e5   missing    missing       missing    missing    missing
 523070.0         missing    missing       missing    missing    missing
      โ‹ฎ                                โ‹ฑ                        โ‹ฎ
      5.38676e5   missing    missing       missing    missing    missing
 538692.0         missing    missing       missing    missing    missing
      5.38708e5   missing    missing   โ€ฆ   missing    missing    missing

In our second variant of point rasterization, we count the number of bike hire stations. To do that, we use the fixed value of 1 (same as in the last example), but this time combined with the reducer=sum argument. That way, multiple values burned into the same pixel are summed, rather than replaced keeping last (which is the default). The new output, ch_raster2, shows the number of cycle hire points in each grid cell.

ch_raster2 = Rasters.rasterize(
    sum,                     # reducer
    cycle_hire_osm_projected; # data
    fill = 1, 
    size = (1000, 1000) # specify size in "pixels"
)
โ•ญโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ•ฎ
โ”‚ 1000ร—1000 Raster{Union{Missing, Int64},2} sum โ”‚
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”ดโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ dims โ”
  โ†“ X Projected{Float64} 523038.61452275474:15.684723604717:538707.653403867 ForwardOrdered Regular Intervals{Start},
  โ†’ Y Projected{Float64} 174934.65727249475:10.03675127048517:184961.37179170942 ForwardOrdered Regular Intervals{Start}
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ metadata โ”ค
  Metadata of Dict{Any, Any}()
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ raster โ”ค
  extent: Extent(X = (523038.61452275474, 538723.3381274717), Y = (174934.65727249475, 184971.40854297992))
  missingval: missing
โ””โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”˜
      โ†“ โ†’        1.74935e5  1.74945e5  โ€ฆ  1.84941e5  1.84951e5  1.84961e5
      5.23039e5   missing    missing       missing    missing    missing
      5.23054e5   missing    missing       missing    missing    missing
 523070.0         missing    missing       missing    missing    missing
      โ‹ฎ                                โ‹ฑ                        โ‹ฎ
      5.38676e5   missing    missing       missing    missing    missing
 538692.0         missing    missing       missing    missing    missing
      5.38708e5   missing    missing   โ€ฆ   missing    missing    missing

The cycle hire locations have different numbers of bicycles described by the capacity variable, raising the question, what is the capacity in each grid cell? To calculate that, in our third point rasterization variant we sum the field ('capacity') rather than the fixed values of 1.

This is extremely simple to run, but we will show how to do this two ways: first, by passing the column name in the feature collection to fill.

ch_raster3 = Rasters.rasterize(
    sum,                     # reducer
    cycle_hire_osm_projected; # data
    fill = :capacity, 
    size = (1000, 1000) # specify size in "pixels"
)
โ•ญโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ•ฎ
โ”‚ 1000ร—1000 Raster{Union{Missing, Float64},2} capacity โ”‚
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”ดโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ dims โ”
  โ†“ X Projected{Float64} 523038.61452275474:15.684723604717:538707.653403867 ForwardOrdered Regular Intervals{Start},
  โ†’ Y Projected{Float64} 174934.65727249475:10.03675127048517:184961.37179170942 ForwardOrdered Regular Intervals{Start}
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ metadata โ”ค
  Metadata of Dict{Any, Any}()
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ raster โ”ค
  extent: Extent(X = (523038.61452275474, 538723.3381274717), Y = (174934.65727249475, 184971.40854297992))
  missingval: missing
โ””โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”˜
      โ†“ โ†’        1.74935e5  1.74945e5  โ€ฆ  1.84941e5  1.84951e5  1.84961e5
      5.23039e5   missing    missing       missing    missing    missing
      5.23054e5   missing    missing       missing    missing    missing
 523070.0         missing    missing       missing    missing    missing
      โ‹ฎ                                โ‹ฑ                        โ‹ฎ
      5.38676e5   missing    missing       missing    missing    missing
 538692.0         missing    missing       missing    missing    missing
      5.38708e5   missing    missing   โ€ฆ   missing    missing    missing

Second, by passing the vectors of geometries and values separately.

ch_raster3 = Rasters.rasterize(
    sum,                     # reducer
    cycle_hire_osm_projected.geom; # data
    fill = cycle_hire_osm_projected.capacity, 
    crs = GI.crs(cycle_hire_osm_projected),
    size = (1000, 1000) # specify size in "pixels"
)
โ•ญโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ•ฎ
โ”‚ 1000ร—1000 Raster{Union{Missing, Float64},2} sum โ”‚
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”ดโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ dims โ”
  โ†“ X Projected{Float64} 523038.61452275474:15.684723604717:538707.653403867 ForwardOrdered Regular Intervals{Start},
  โ†’ Y Projected{Float64} 174934.65727249475:10.03675127048517:184961.37179170942 ForwardOrdered Regular Intervals{Start}
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ metadata โ”ค
  Metadata of Dict{Any, Any}()
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ raster โ”ค
  extent: Extent(X = (523038.61452275474, 538723.3381274717), Y = (174934.65727249475, 184971.40854297992))
  missingval: missing
  crs: EPSG:27700
โ””โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”˜
      โ†“ โ†’        1.74935e5  1.74945e5  โ€ฆ  1.84941e5  1.84951e5  1.84961e5
      5.23039e5   missing    missing       missing    missing    missing
      5.23054e5   missing    missing       missing    missing    missing
 523070.0         missing    missing       missing    missing    missing
      โ‹ฎ                                โ‹ฑ                        โ‹ฎ
      5.38676e5   missing    missing       missing    missing    missing
 538692.0         missing    missing       missing    missing    missing
      5.38708e5   missing    missing   โ€ฆ   missing    missing    missing

The input point layer cycle_hire_osm_projected and the three variants of rasterizing it ch_raster1, ch_raster2, and ch_raster3 are shown in Figure 5.5.

# Input points
nonmissing_df = dropmissing(cycle_hire_osm_projected, [:capacity, :geom])
f, a, p = scatter(nonmissing_df.geom; color = nonmissing_df.capacity)
Colorbar(f[1, 2], p)
display(f)
# Presence/Absence
plot(ch_raster1) |> display
# Point counts
plot(ch_raster2) |> display
# Summed attribute values
plot(ch_raster3)
(a) Input points
(b) Presence/Absence
(c) Point counts
(d) Summed attribute values
Figure 5.5: Original data and three variants of point rasterization

5.4.2 Rasterizing lines and polygons

Another dataset based on Californiaโ€™s polygons and borders (created below) illustrates rasterization of lines. There are three preliminary steps. First, we subset the California polygon.

california = us_states[ us_states[!, "NAME"] .== "California", :]
1ร—7 DataFrame
Row geom GEOID NAME REGION AREA total_pop_10 total_pop_15
IGeometrโ€ฆ String String String Float64 Float64 Float64
1 Geometry: wkbMultiPolygon 06 California West 4.09747e5 3.66373e7 3.84215e7

Second, we obtain the borders of the polygon as a `โ€˜MultiLineStringโ€™

california_geom = only(california.geom)
california_borders = GI.MultiLineString(GI.LineString.(GI.getexterior.(GI.getgeom(california_geom))); crs = GI.crs(california_geom)) # TODO: make this a lot better....
GeoInterface.Wrappers.MultiLineString{false, false, Vector{GeoInterface.Wrappers.LineString{false, false, ArchGDAL.IGeometry{ArchGDAL.wkbLineString}, Nothing, Nothing}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LineString{false, false, ArchGDAL.IGeometry{ArchGDAL.wkbLineString}, Nothing, Nothing}[GeoInterface.Wrappers.LineString{false, false, ArchGDAL.IGeometry{ArchGDAL.wkbLineString}, Nothing, Nothing}(Geometry: LINEARRING (-118.603375 33.478098,-118.368301 33.4 ... 8098), nothing, nothing), GeoInterface.Wrappers.LineString{false, false, ArchGDAL.IGeometry{ArchGDAL.wkbLineString}, Nothing, Nothing}(Geometry: LINEARRING (-118.605534 33.030999,-118.369984 32.8 ... 0999), nothing, nothing), GeoInterface.Wrappers.LineString{false, false, ArchGDAL.IGeometry{ArchGDAL.wkbLineString}, Nothing, Nothing}(Geometry: LINEARRING (-119.919155 34.07728,-119.755521 34.05 ... 7728), nothing, nothing), GeoInterface.Wrappers.LineString{false, false, ArchGDAL.IGeometry{ArchGDAL.wkbLineString}, Nothing, Nothing}(Geometry: LINEARRING (-120.247393 34.001911,-120.073609 34.0 ... 1911), nothing, nothing), GeoInterface.Wrappers.LineString{false, false, ArchGDAL.IGeometry{ArchGDAL.wkbLineString}, Nothing, Nothing}(Geometry: LINEARRING (-124.211605 41.99846,-123.347562 41.99 ... 9846), nothing, nothing)], nothing, WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"NAD83\",DATUM[\"North_American_Datum_1983\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"6269\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4269\"]]"))

Finally, we rasterize california_borders on a grid with resolution of 0.5 degrees per pixel.

california_raster1 = Rasters.rasterize(
    last,
    california_borders;
    fill = 1,
    res = 0.5, # degrees - this is in units of GI.crs(california_borders)
    boundary = :touches,
)
โ•ญโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ•ฎ
โ”‚ 21ร—19 Raster{Union{Missing, Int64},2} last โ”‚
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”ดโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ dims โ”
  โ†“ X Projected{Float64} -124.409591:0.5:-114.409591 ForwardOrdered Regular Intervals{Start},
  โ†’ Y Projected{Float64} 32.534156:0.5:41.534156 ForwardOrdered Regular Intervals{Start}
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ metadata โ”ค
  Metadata of Dict{Any, Any}()
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ raster โ”ค
  extent: Extent(X = (-124.409591, -113.909591), Y = (32.534156, 42.034156))
  missingval: missing
โ””โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”˜
    โ†“ โ†’   32.5342    33.0342    33.5342    โ€ฆ  40.5342    41.0342    41.5342
 -124.41    missing    missing    missing      1          1          1
 -123.91    missing    missing    missing       missing    missing   1
 -123.41    missing    missing    missing       missing    missing   1
    โ‹ฎ                                      โ‹ฑ                        
 -115.41   1           missing    missing       missing    missing    missing
 -114.91   1          1          1         โ€ฆ    missing    missing    missing
 -114.41    missing    missing    missing       missing    missing    missing

Compare it to a polygon rasterization, with all_touched=False (the default), which selects only raster cells whose centroids are inside the selector polygon, as illustrated in Figure 5.6 (right).

california_raster2 = Rasters.rasterize(
    last, 
    california;
    geometrycolumn = :geom,
    fill = 1,
    res = 0.5,
    boundary = :center,
)
โ•ญโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ•ฎ
โ”‚ 21ร—19 Raster{Union{Missing, Int64},2} last โ”‚
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”ดโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ dims โ”
  โ†“ X Projected{Float64} -124.409591:0.5:-114.409591 ForwardOrdered Regular Intervals{Start},
  โ†’ Y Projected{Float64} 32.534156:0.5:41.534156 ForwardOrdered Regular Intervals{Start}
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ metadata โ”ค
  Metadata of Dict{Any, Any} with 1 entry:
  "missed_geometries" => Bool[0]
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ raster โ”ค
  extent: Extent(X = (-124.409591, -113.909591), Y = (32.534156, 42.034156))
  missingval: missing
โ””โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”˜
    โ†“ โ†’   32.5342    33.0342    33.5342    โ€ฆ  40.5342    41.0342    41.5342
 -124.41    missing    missing    missing      1           missing   1
 -123.91    missing    missing    missing      1          1          1
    โ‹ฎ                                      โ‹ฑ                        
 -115.41   1          1          1              missing    missing    missing
 -114.91   1           missing   1         โ€ฆ    missing    missing    missing
 -114.41    missing    missing    missing       missing    missing    missing

To illustrate which raster pixels are actually selected as part of rasterization, we also show them as points. This also requires the following code section to calculate the points, which we explain in Section 5.5.

dp = DimPoints(california_raster1)
โ•ญโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ•ฎ
โ”‚ 21ร—19 DimPoints{Tuple{Float64, Float64},2} โ”‚
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”ดโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ dims โ”
  โ†“ X Projected{Float64} -124.409591:0.5:-114.409591 ForwardOrdered Regular Intervals{Start},
  โ†’ Y Projected{Float64} 32.534156:0.5:41.534156 ForwardOrdered Regular Intervals{Start}
โ””โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”˜
    โ†“ โ†’   32.5342               33.0342               โ€ฆ  41.5342
 -124.41    (-124.41, 32.5342)    (-124.41, 33.0342)       (-124.41, 41.5342)
 -123.91    (-123.91, 32.5342)    (-123.91, 33.0342)       (-123.91, 41.5342)
 -123.41    (-123.41, 32.5342)    (-123.41, 33.0342)       (-123.41, 41.5342)
 -122.91    (-122.91, 32.5342)    (-122.91, 33.0342)       (-122.91, 41.5342)
    โ‹ฎ                                                 โ‹ฑ  
 -116.41    (-116.41, 32.5342)    (-116.41, 33.0342)       (-116.41, 41.5342)
 -115.91    (-115.91, 32.5342)    (-115.91, 33.0342)       (-115.91, 41.5342)
 -115.41    (-115.41, 32.5342)    (-115.41, 33.0342)       (-115.41, 41.5342)
 -114.91    (-114.91, 32.5342)    (-114.91, 33.0342)  โ€ฆ    (-114.91, 41.5342)
 -114.41    (-114.41, 32.5342)    (-114.41, 33.0342)       (-114.41, 41.5342)

DimensionalData.jl (which underpins Rasters.jl) provides easy ways to get โ€œlookupsโ€, i.e, axis index values, from a raster.

Note that these lookups may encode intervals, points, or anything in between - so you should use shiftlocus or set to get the actual point values! DimPoints does this for you.

But if you want to see how this can be done automatically, here you go.

[(x, y) for x in dims(california_raster1, X), y in dims(california_raster1, Y)]
21ร—19 Matrix{Tuple{Float64, Float64}}:
 (-124.41, 32.5342)  (-124.41, 33.0342)  โ€ฆ  (-124.41, 41.5342)
 (-123.91, 32.5342)  (-123.91, 33.0342)     (-123.91, 41.5342)
 (-123.41, 32.5342)  (-123.41, 33.0342)     (-123.41, 41.5342)
 (-122.91, 32.5342)  (-122.91, 33.0342)     (-122.91, 41.5342)
 (-122.41, 32.5342)  (-122.41, 33.0342)     (-122.41, 41.5342)
 (-121.91, 32.5342)  (-121.91, 33.0342)  โ€ฆ  (-121.91, 41.5342)
 (-121.41, 32.5342)  (-121.41, 33.0342)     (-121.41, 41.5342)
 (-120.91, 32.5342)  (-120.91, 33.0342)     (-120.91, 41.5342)
 (-120.41, 32.5342)  (-120.41, 33.0342)     (-120.41, 41.5342)
 (-119.91, 32.5342)  (-119.91, 33.0342)     (-119.91, 41.5342)
 โ‹ฎ                                       โ‹ฑ  
 (-118.41, 32.5342)  (-118.41, 33.0342)     (-118.41, 41.5342)
 (-117.91, 32.5342)  (-117.91, 33.0342)     (-117.91, 41.5342)
 (-117.41, 32.5342)  (-117.41, 33.0342)     (-117.41, 41.5342)
 (-116.91, 32.5342)  (-116.91, 33.0342)  โ€ฆ  (-116.91, 41.5342)
 (-116.41, 32.5342)  (-116.41, 33.0342)     (-116.41, 41.5342)
 (-115.91, 32.5342)  (-115.91, 33.0342)     (-115.91, 41.5342)
 (-115.41, 32.5342)  (-115.41, 33.0342)     (-115.41, 41.5342)
 (-114.91, 32.5342)  (-114.91, 33.0342)     (-114.91, 41.5342)
 (-114.41, 32.5342)  (-114.41, 33.0342)  โ€ฆ  (-114.41, 41.5342)

You can see that this encodes the same values as dp in the other tab.

TODO: firm up the description here and add links.

Figure 5.6 shows the input vector layer, the rasterization results, and the points pnt.

# Line rasterization
fig, ax, plt = plot(california_raster1; colormap = cgrad(:Set3; categorical = true))
lines!(ax, california_borders; color = :darkgrey, linewidth = 1)
scatter!(ax, vec(dp); markersize = 3, color = :black)
display(fig)
# Polygon rasterization
fig, ax, plt = plot(california_raster2; colormap = cgrad(:Set3; categorical = true))
lines!(ax, california_borders; color = :darkgrey, linewidth = 1)
scatter!(ax, vec(dp); markersize = 3, color = :black)
fig
(a) Line rasterization w/ boundary=:touches
(b) Polygon rasterization w/ boundary=:center
Figure 5.6: Examples of line and polygon rasterization

5.5 Spatial vectorization

Spatial vectorization is the counterpart of rasterization (Section 5.4). It involves converting spatially continuous raster data into spatially discrete vector data such as points, lines or polygons. There are three standard methods to convert a raster to a vector layer, which we cover next:

  • Raster to polygons (Section 5.5.1)โ€”converting raster cells to rectangular polygons, representing pixel areas
  • Raster to points (Section 5.5.2)โ€”converting raster cells to points, representing pixel centroids
  • Raster to contours (Section 5.5.3)

Let us demonstrate all three in the given order.

5.5.1 Raster to polygons

Rasters.jl does not currently have a function to convert a raster to a feature collection with one polygon per pixel or cell. This is a similar situation in Python with rasterio.

GeometryOps.jl offers a polygonize function that returns a feature collection of polygons, where each feature has a value property that encodes the value of all pixels within that polygon. Each polygon contains pixels with the same value.

fc = GO.polygonize(src_grain)
GeoInterface.Wrappers.FeatureCollection{Vector{GeoInterface.Wrappers.Feature{@NamedTuple{geometry::GeoInterface.Wrappers.MultiPolygon{false, false, Vector{GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}, value::UInt8}, WellKnownText{GeoFormatTypes.CRS}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}}}, WellKnownText{GeoFormatTypes.CRS}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}}(GeoInterface.Wrappers.Feature{@NamedTuple{geometry::GeoInterface.Wrappers.MultiPolygon{false, false, Vector{GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}, value::UInt8}, WellKnownText{GeoFormatTypes.CRS}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}}[GeoInterface.Wrappers.Feature{@NamedTuple{geometry::GeoInterface.Wrappers.MultiPolygon{false, false, Vector{GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}, value::UInt8}, WellKnownText{GeoFormatTypes.CRS}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}}((geometry = GeoInterface.Wrappers.MultiPolygon{false, false, Vector{GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}([(1.5, 1.5), (1.5, 3.5), (2.5, 3.5), (2.5, 4.5), (0.5, 4.5), (0.5, 1.5), (1.5, 1.5)], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 2.5), Y = (1.5, 4.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 2.5), Y = (1.5, 4.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]")), GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}([(1.5, 1.5), (1.5, 0.5), (2.5, 0.5), (2.5, 1.5), (1.5, 1.5)], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (1.5, 2.5), Y = (0.5, 1.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (1.5, 2.5), Y = (0.5, 1.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]")), GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}([(5.5, 5.5), (5.5, 6.5), (4.5, 6.5), (4.5, 5.5), (5.5, 5.5)], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (4.5, 5.5), Y = (5.5, 6.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (4.5, 5.5), Y = (5.5, 6.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]")), GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}([(4.5, 3.5), (3.5, 3.5), (3.5, 2.5), (2.5, 2.5), (2.5, 1.5), (4.5, 1.5), (4.5, 2.5), (5.5, 2.5), (5.5, 3.5), (4.5, 3.5)], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (2.5, 5.5), Y = (1.5, 3.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (2.5, 5.5), Y = (1.5, 3.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 5.5), Y = (0.5, 6.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]")), value = 0x00), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"), Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 5.5), Y = (0.5, 6.5)))), GeoInterface.Wrappers.Feature{@NamedTuple{geometry::GeoInterface.Wrappers.MultiPolygon{false, false, Vector{GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}, value::UInt8}, WellKnownText{GeoFormatTypes.CRS}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}}((geometry = GeoInterface.Wrappers.MultiPolygon{false, false, Vector{GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}([(0.5, 4.5), (2.5, 4.5), (2.5, 3.5), (6.5, 3.5), (6.5, 5.5), (4.5, 5.5), (4.5, 4.5), (3.5, 4.5), (3.5, 5.5), (2.5, 5.5), (2.5, 6.5), (1.5, 6.5), (1.5, 5.5), (0.5, 5.5), (0.5, 4.5)], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 6.5), Y = (3.5, 6.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 6.5), Y = (3.5, 6.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]")), GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}([(1.5, 1.5), (0.5, 1.5), (0.5, 0.5), (1.5, 0.5), (1.5, 1.5)], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 1.5), Y = (0.5, 1.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 1.5), Y = (0.5, 1.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]")), GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}([(5.5, 2.5), (5.5, 1.5), (6.5, 1.5), (6.5, 2.5), (5.5, 2.5)], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (5.5, 6.5), Y = (1.5, 2.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (5.5, 6.5), Y = (1.5, 2.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]")), GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}([(2.5, 0.5), (3.5, 0.5), (3.5, 1.5), (2.5, 1.5), (2.5, 0.5)], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (2.5, 3.5), Y = (0.5, 1.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (2.5, 3.5), Y = (0.5, 1.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 6.5), Y = (0.5, 6.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]")), value = 0x01), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"), Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 6.5), Y = (0.5, 6.5)))), GeoInterface.Wrappers.Feature{@NamedTuple{geometry::GeoInterface.Wrappers.MultiPolygon{false, false, Vector{GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}, value::UInt8}, WellKnownText{GeoFormatTypes.CRS}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}}((geometry = GeoInterface.Wrappers.MultiPolygon{false, false, Vector{GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}([(1.5, 1.5), (2.5, 1.5), (2.5, 2.5), (3.5, 2.5), (3.5, 3.5), (1.5, 3.5), (1.5, 1.5)], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (1.5, 3.5), Y = (1.5, 3.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (1.5, 3.5), Y = (1.5, 3.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]")), GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}([(5.5, 5.5), (6.5, 5.5), (6.5, 6.5), (5.5, 6.5), (5.5, 5.5)], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (5.5, 6.5), Y = (5.5, 6.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (5.5, 6.5), Y = (5.5, 6.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]")), GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}([(4.5, 6.5), (2.5, 6.5), (2.5, 5.5), (3.5, 5.5), (3.5, 4.5), (4.5, 4.5), (4.5, 6.5)], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (2.5, 4.5), Y = (4.5, 6.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (2.5, 4.5), Y = (4.5, 6.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]")), GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}([(5.5, 0.5), (6.5, 0.5), (6.5, 1.5), (5.5, 1.5), (5.5, 2.5), (4.5, 2.5), (4.5, 1.5), (3.5, 1.5), (3.5, 0.5), (5.5, 0.5)], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (3.5, 6.5), Y = (0.5, 2.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (3.5, 6.5), Y = (0.5, 2.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]")), GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}([(6.5, 3.5), (5.5, 3.5), (5.5, 2.5), (6.5, 2.5), (6.5, 3.5)], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (5.5, 6.5), Y = (2.5, 3.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (5.5, 6.5), Y = (2.5, 3.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]")), GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}([(1.5, 6.5), (0.5, 6.5), (0.5, 5.5), (1.5, 5.5), (1.5, 6.5)], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 1.5), Y = (5.5, 6.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 1.5), Y = (5.5, 6.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 6.5), Y = (0.5, 6.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]")), value = 0x02), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"), Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 6.5), Y = (0.5, 6.5))))], WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"), Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 6.5), Y = (0.5, 6.5))))

We can convert this โ€œfeature collectionโ€ to a DataFrame as follows. That makes it a lot easier to work with.

df = DataFrame([GI.properties(f) for f in GI.getfeature(fc)])
df.geometry = [GI.geometry(f) for f in GI.getfeature(fc)]
df
3ร—2 DataFrame
Row value geometry
UInt8 MultiPolโ€ฆ
1 0 MultiPolygon{false, false, Vector{Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}([(1.5, 1.5), (1.5, 3.5), (2.5, 3.5), (2.5, 4.5), (0.5, 4.5), (0.5, 1.5), (1.5, 1.5)], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 2.5), Y = (1.5, 4.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 2.5), Y = (1.5, 4.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]")), Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}([(1.5, 1.5), (1.5, 0.5), (2.5, 0.5), (2.5, 1.5), (1.5, 1.5)], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (1.5, 2.5), Y = (0.5, 1.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (1.5, 2.5), Y = (0.5, 1.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]")), Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}([(5.5, 5.5), (5.5, 6.5), (4.5, 6.5), (4.5, 5.5), (5.5, 5.5)], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (4.5, 5.5), Y = (5.5, 6.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (4.5, 5.5), Y = (5.5, 6.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]")), Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}([(4.5, 3.5), (3.5, 3.5), (3.5, 2.5), (2.5, 2.5), (2.5, 1.5), (4.5, 1.5), (4.5, 2.5), (5.5, 2.5), (5.5, 3.5), (4.5, 3.5)], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (2.5, 5.5), Y = (1.5, 3.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (2.5, 5.5), Y = (1.5, 3.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 5.5), Y = (0.5, 6.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))
2 1 MultiPolygon{false, false, Vector{Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}([(0.5, 4.5), (2.5, 4.5), (2.5, 3.5), (6.5, 3.5), (6.5, 5.5), (4.5, 5.5), (4.5, 4.5), (3.5, 4.5), (3.5, 5.5), (2.5, 5.5), (2.5, 6.5), (1.5, 6.5), (1.5, 5.5), (0.5, 5.5), (0.5, 4.5)], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 6.5), Y = (3.5, 6.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 6.5), Y = (3.5, 6.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]")), Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}([(1.5, 1.5), (0.5, 1.5), (0.5, 0.5), (1.5, 0.5), (1.5, 1.5)], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 1.5), Y = (0.5, 1.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 1.5), Y = (0.5, 1.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]")), Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}([(5.5, 2.5), (5.5, 1.5), (6.5, 1.5), (6.5, 2.5), (5.5, 2.5)], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (5.5, 6.5), Y = (1.5, 2.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (5.5, 6.5), Y = (1.5, 2.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]")), Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}([(2.5, 0.5), (3.5, 0.5), (3.5, 1.5), (2.5, 1.5), (2.5, 0.5)], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (2.5, 3.5), Y = (0.5, 1.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (2.5, 3.5), Y = (0.5, 1.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 6.5), Y = (0.5, 6.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))
3 2 MultiPolygon{false, false, Vector{Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}([(1.5, 1.5), (2.5, 1.5), (2.5, 2.5), (3.5, 2.5), (3.5, 3.5), (1.5, 3.5), (1.5, 1.5)], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (1.5, 3.5), Y = (1.5, 3.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (1.5, 3.5), Y = (1.5, 3.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]")), Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}([(5.5, 5.5), (6.5, 5.5), (6.5, 6.5), (5.5, 6.5), (5.5, 5.5)], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (5.5, 6.5), Y = (5.5, 6.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (5.5, 6.5), Y = (5.5, 6.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]")), Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}([(4.5, 6.5), (2.5, 6.5), (2.5, 5.5), (3.5, 5.5), (3.5, 4.5), (4.5, 4.5), (4.5, 6.5)], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (2.5, 4.5), Y = (4.5, 6.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (2.5, 4.5), Y = (4.5, 6.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]")), Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}([(5.5, 0.5), (6.5, 0.5), (6.5, 1.5), (5.5, 1.5), (5.5, 2.5), (4.5, 2.5), (4.5, 1.5), (3.5, 1.5), (3.5, 0.5), (5.5, 0.5)], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (3.5, 6.5), Y = (0.5, 2.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (3.5, 6.5), Y = (0.5, 2.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]")), Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}([(6.5, 3.5), (5.5, 3.5), (5.5, 2.5), (6.5, 2.5), (6.5, 3.5)], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (5.5, 6.5), Y = (2.5, 3.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (5.5, 6.5), Y = (2.5, 3.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]")), Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}([(1.5, 6.5), (0.5, 6.5), (0.5, 5.5), (1.5, 5.5), (1.5, 6.5)], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 1.5), Y = (5.5, 6.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 1.5), Y = (5.5, 6.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 6.5), Y = (0.5, 6.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))

The polygon layer df is shown in Figure 5.7.

f, a, p = poly(df.geometry; color = df.value, strokecolor = :black, strokewidth = 0.75)
Colorbar(f[1, 2], p)
f
Figure 5.7: grain.tif converted to a polygon layer

As highlighted using edgecolor='black', neighboring pixels sharing the same raster value are dissolved into larger polygons.

One suggestion is to add unique values between 0 and 0.9999 to all pixels, convert to polygons, and then get back to the original values using floor.

5.5.2 Raster to points

To transform a raster to points, Rasters.jl provides the Rasters.DimTable constructor, which converts a raster into a lazy, table-like form. This can be converted directly to a DataFrame, or operated on independently.

dt = DimTable(Raster("output/elev.tif"))
DimensionalData.DimTable with 36 rows, 3 columns, and schema:
 :X       Float64
 :Y       Float64
 :layer1  UInt8

Notice that this has three columns, :X, :Y, and :layer1, corresponding to the pixel centroids and elevation values. But what if we want to treat the X and Y dimensionas as point geometries?

DimTable has a mergedims keyword argument for this, which allows us to merge the X and Y dimensions into a single dimension.

dt = DimTable(Raster("output/elev.tif"), mergedims = (X, Y))
DimensionalData.DimTable with 36 rows, 2 columns, and schema:
 :XY      Tuple{Float64, Float64}
 :layer1  UInt8

This has created a DimTable with a column :XY, which contains the pixel centroids as point-like objects. We can convert this to a DataFrame, set some metadata to indicate that geometry is in :XY, and plot the result.

df = DataFrame(dt)
DataFrames.metadata!(df, "GEOINTERFACE:geometrycolumns", (:XY,); style = :note)
df
36ร—2 DataFrame
11 rows omitted
Row XY layer1
Tupleโ€ฆ UInt8
1 (-1.5, 1.0) 1
2 (-1.0, 1.0) 2
3 (-0.5, 1.0) 3
4 (-1.11022e-16, 1.0) 4
5 (0.5, 1.0) 5
6 (1.0, 1.0) 6
7 (-1.5, 0.5) 7
8 (-1.0, 0.5) 8
9 (-0.5, 0.5) 9
10 (-1.11022e-16, 0.5) 10
11 (0.5, 0.5) 11
12 (1.0, 0.5) 12
13 (-1.5, -1.11022e-16) 13
โ‹ฎ โ‹ฎ โ‹ฎ
25 (-1.5, -1.0) 25
26 (-1.0, -1.0) 26
27 (-0.5, -1.0) 27
28 (-1.11022e-16, -1.0) 28
29 (0.5, -1.0) 29
30 (1.0, -1.0) 30
31 (-1.5, -1.5) 31
32 (-1.0, -1.5) 32
33 (-0.5, -1.5) 33
34 (-1.11022e-16, -1.5) 34
35 (0.5, -1.5) 35
36 (1.0, -1.5) 36
scatter(df.XY; color = df.layer1)

We can even save this to a file trivially easily:

GeoDataFrames.write("output/elev.gpkg", df)
GeoDataFrames.read("output/elev.gpkg")
36ร—2 DataFrame
11 rows omitted
Row XY layer1
IGeometrโ€ฆ Int16
1 Geometry: wkbPoint 1
2 Geometry: wkbPoint 2
3 Geometry: wkbPoint 3
4 Geometry: wkbPoint 4
5 Geometry: wkbPoint 5
6 Geometry: wkbPoint 6
7 Geometry: wkbPoint 7
8 Geometry: wkbPoint 8
9 Geometry: wkbPoint 9
10 Geometry: wkbPoint 10
11 Geometry: wkbPoint 11
12 Geometry: wkbPoint 12
13 Geometry: wkbPoint 13
โ‹ฎ โ‹ฎ โ‹ฎ
25 Geometry: wkbPoint 25
26 Geometry: wkbPoint 26
27 Geometry: wkbPoint 27
28 Geometry: wkbPoint 28
29 Geometry: wkbPoint 29
30 Geometry: wkbPoint 30
31 Geometry: wkbPoint 31
32 Geometry: wkbPoint 32
33 Geometry: wkbPoint 33
34 Geometry: wkbPoint 34
35 Geometry: wkbPoint 35
36 Geometry: wkbPoint 36

Figure 5.8 shows the input raster and the resulting point layer.

# Input raster
fig, ax, plt = plot(src_elev)
scatter!(ax, df.XY; color = df.layer1)
display(fig)
# Points
fig, ax, plt = plot(src_elev; alpha = 0.1)
scatter!(ax, df.XY; color = df.layer1, strokecolor = :black, strokewidth = 1)
fig
(a) Input raster
(b) Points
Figure 5.8: Raster and point representation of elev.tif

TODO: nodata pixels

5.5.3 Raster to contours

Another common type of spatial vectorization is the creation of contour lines, representing lines of continuous height or temperatures (isotherms), for example. We will use a real-world digital elevation model (DEM) because the artificial raster elev.tif produces parallel lines (task for the reader: verify this and explain why this happens). Plotting contour lines is straightforward, using the contour or contourf functions in Makie.

f, ax, plt = contour(src_dem; levels = LinRange(0, 1200, 50), color = :black)

TODO: gdal_contour (via ArchGDAL??)

It would be good to show how to use the provided GDAL executables thoughโ€ฆ

5.6 Distance to nearest geometry

Calculating a raster of distances to the nearest geometry is an example of a โ€œglobalโ€ raster operation (?sec-global-operations-and-distances). To demonstrate it, suppose that we need to calculate a raster representing the distance to the nearest coast in New Zealand. This example also wraps many of the concepts introduced in this chapter and in previous chapters, such as raster aggregation (?sec-raster-agg-disagg), raster conversion to points (Section 5.5.2), and rasterizing points (Section 5.4.1).

For the coastline, we will dissolve the New Zealand administrative division polygon layer and โ€œextractโ€ the boundary as a 'MultiLineString' geometry.

using LibGEOS
coastline_linestrings = GI.getexterior.(GI.getgeom(LibGEOS.unaryUnion(GI.GeometryCollection(nz.geom)))) .|> x -> GI.LineString(collect(GI.getpoint(x)))
coastline = GI.MultiLineString(coastline_linestrings)
coastline = GO.reproject(coastline; target_crs = GI.crs(src_nz_elev), source_crs = GI.crs(nz))
GeoInterface.Wrappers.MultiLineString{false, false, Vector{GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}([(1.2299979013e6, 4.7934322369e6), (1.2216031693e6, 4.7943017685e6), (1.2177895399e6, 4.7904898277e6), (1.2311152067e6, 4.787137924e6), (1.2361839211e6, 4.7843314255e6), (1.2325099445e6, 4.7718700237e6), (1.2173644933e6, 4.7681066692e6), (1.2080683451e6, 4.7605677875e6), (1.199129799e6, 4.7634183939e6), (1.1938565725e6, 4.7574957427e6)  โ€ฆ  (1.1964031733e6, 4.7848147499e6), (1.2011220591e6, 4.7884807601e6), (1.2003610396e6, 4.7969222422e6), (1.1955720567e6, 4.8084895651e6), (1.1989855975e6, 4.8153873633e6), (1.2077044806e6, 4.8171297889e6), (1.2161546733e6, 4.8132662203e6), (1.2184777781e6, 4.8064984488e6), (1.2297297347e6, 4.7987186067e6), (1.2299979013e6, 4.7934322369e6)], nothing, WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "PROJCS[\"unknown\",GEOGCS[\"unknown\",DATUM[\"Unknown based on GRS80 ellipsoid\",SPHEROID[\"GRS 1980\",6378137,298.257222101004,AUTHORITY[\"EPSG\",\"7019\"]],TOWGS84[0,0,0,0,0,0,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",173],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",1600000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]")), GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}([(1.1010682698e6, 4.9169506843e6), (1.0962776826e6, 4.9262213416e6), (1.1070468733e6, 4.9308841834e6), (1.1101943531e6, 4.9284267394e6), (1.1114045691e6, 4.9201779986e6), (1.1088006957e6, 4.914502166e6), (1.1010682698e6, 4.9169506843e6)], nothing, WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "PROJCS[\"unknown\",GEOGCS[\"unknown\",DATUM[\"Unknown based on GRS80 ellipsoid\",SPHEROID[\"GRS 1980\",6378137,298.257222101004,AUTHORITY[\"EPSG\",\"7019\"]],TOWGS84[0,0,0,0,0,0,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",173],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",1600000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]")), GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}([(1.2822947146e6, 4.8243903885e6), (1.2777592139e6, 4.8327835132e6), (1.2629460033e6, 4.8352306146e6), (1.2563895979e6, 4.832365038e6), (1.2436815529e6, 4.8358664271e6), (1.237141582e6, 4.8332854896e6), (1.2365580519e6, 4.8400088904e6), (1.2301828929e6, 4.8518834763e6), (1.2202312165e6, 4.8564065645e6), (1.2167650225e6, 4.8512224604e6)  โ€ฆ  (1.3536635084e6, 4.861558587e6), (1.3540190062e6, 4.8508162544e6), (1.3384831555e6, 4.837432841e6), (1.3299957273e6, 4.8375532567e6), (1.3247154647e6, 4.8316347032e6), (1.3151459743e6, 4.8300668586e6), (1.310825215e6, 4.826096002e6), (1.3025458256e6, 4.8269904113e6), (1.2945187261e6, 4.8232725043e6), (1.2822947146e6, 4.8243903885e6)], nothing, WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "PROJCS[\"unknown\",GEOGCS[\"unknown\",DATUM[\"Unknown based on GRS80 ellipsoid\",SPHEROID[\"GRS 1980\",6378137,298.257222101004,AUTHORITY[\"EPSG\",\"7019\"]],TOWGS84[0,0,0,0,0,0,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",173],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",1600000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]")), GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}([(1.1253763259e6, 4.982869816e6), (1.1248143323e6, 4.976641044e6), (1.1282897921e6, 4.9707170423e6), (1.1254201958e6, 4.9665221375e6), (1.1188173258e6, 4.9720125867e6), (1.1253763259e6, 4.982869816e6)], nothing, WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "PROJCS[\"unknown\",GEOGCS[\"unknown\",DATUM[\"Unknown based on GRS80 ellipsoid\",SPHEROID[\"GRS 1980\",6378137,298.257222101004,AUTHORITY[\"EPSG\",\"7019\"]],TOWGS84[0,0,0,0,0,0,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",173],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",1600000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]")), GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}([(1.7015121629e6, 5.9962048743e6), (1.6940719249e6, 5.9966701668e6), (1.6984728948e6, 5.9883324659e6), (1.7037685647e6, 5.9848085328e6), (1.7060307416e6, 5.9751359042e6), (1.7018195613e6, 5.9719674757e6), (1.6933144403e6, 5.9719708739e6), (1.6883953701e6, 5.9896812008e6), (1.6779656958e6, 6.0055889233e6), (1.662101242e6, 6.0277652897e6)  โ€ฆ  (1.725607508e6, 5.945892267e6), (1.7281682052e6, 5.9584185318e6), (1.7268710844e6, 5.9703482785e6), (1.7127870833e6, 5.9730682541e6), (1.7128762566e6, 5.9791634373e6), (1.7232886065e6, 5.9797407952e6), (1.7331650316e6, 5.9897140167e6), (1.7201965583e6, 5.9800776907e6), (1.709110421e6, 5.9866719829e6), (1.7015121629e6, 5.9962048743e6)], nothing, WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "PROJCS[\"unknown\",GEOGCS[\"unknown\",DATUM[\"Unknown based on GRS80 ellipsoid\",SPHEROID[\"GRS 1980\",6378137,298.257222101004,AUTHORITY[\"EPSG\",\"7019\"]],TOWGS84[0,0,0,0,0,0,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",173],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",1600000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]")), GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}([(1.6650793317e6, 5.4762796377e6), (1.6679650479e6, 5.4828286371e6), (1.6790721723e6, 5.4888496631e6), (1.6808040605e6, 5.4949012903e6), (1.6810365469e6, 5.4869445483e6), (1.6763032282e6, 5.4766970825e6), (1.6685295209e6, 5.4686427381e6), (1.6650793317e6, 5.4762796377e6)], nothing, WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "PROJCS[\"unknown\",GEOGCS[\"unknown\",DATUM[\"Unknown based on GRS80 ellipsoid\",SPHEROID[\"GRS 1980\",6378137,298.257222101004,AUTHORITY[\"EPSG\",\"7019\"]],TOWGS84[0,0,0,0,0,0,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",173],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",1600000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]")), GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}([(1.7939818643e6, 5.9319178426e6), (1.7912570345e6, 5.9199068808e6), (1.7842798937e6, 5.9215788138e6), (1.7818307229e6, 5.9258456737e6), (1.7939818643e6, 5.9319178426e6)], nothing, WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "PROJCS[\"unknown\",GEOGCS[\"unknown\",DATUM[\"Unknown based on GRS80 ellipsoid\",SPHEROID[\"GRS 1980\",6378137,298.257222101004,AUTHORITY[\"EPSG\",\"7019\"]],TOWGS84[0,0,0,0,0,0,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",173],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",1600000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]")), GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}([(1.8266292121e6, 5.9745333354e6), (1.8188737655e6, 5.9791838569e6), (1.8185878044e6, 5.9830890239e6), (1.8112269977e6, 5.9881418267e6), (1.8126655001e6, 5.9932466009e6), (1.8093053687e6, 5.9982083619e6), (1.8108957047e6, 6.0059053369e6), (1.8159541251e6, 6.0078779063e6), (1.8180162185e6, 5.9990038299e6), (1.823326993e6, 5.9918794741e6), (1.8227997093e6, 5.9860917836e6), (1.8283790056e6, 5.9803071686e6), (1.8266292121e6, 5.9745333354e6)], nothing, WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "PROJCS[\"unknown\",GEOGCS[\"unknown\",DATUM[\"Unknown based on GRS80 ellipsoid\",SPHEROID[\"GRS 1980\",6378137,298.257222101004,AUTHORITY[\"EPSG\",\"7019\"]],TOWGS84[0,0,0,0,0,0,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",173],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",1600000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]"))], nothing, WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "PROJCS[\"unknown\",GEOGCS[\"unknown\",DATUM[\"Unknown based on GRS80 ellipsoid\",SPHEROID[\"GRS 1980\",6378137,298.257222101004,AUTHORITY[\"EPSG\",\"7019\"]],TOWGS84[0,0,0,0,0,0,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",173],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",1600000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]"))

For a โ€œtemplateโ€ raster, we will aggregate the New Zealand DEM, in the nz_elev.tif file, to 5 times coarser resolution. The code section below follows the aggeregation example in ?sec-raster-agg-disagg.

factor = 2/10
r = Rasters.resample(src_nz_elev; size = round.(Int, size(src_nz_elev) .* factor), method = :average)
โ•ญโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ•ฎ
โ”‚ 223ร—290 Raster{Float32,2} โ”‚
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”ดโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ dims โ”
  โ†“ X Projected{Float64} LinRange{Float64}(995456.5461976258, 2.105456546197626e6, 223) ForwardOrdered Regular Intervals{Start},
  โ†’ Y Projected{Float64} LinRange{Float64}(6.186960927303582e6, 4.741960927303582e6, 290) ReverseOrdered Regular Intervals{Start}
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ metadata โ”ค
  Metadata{Rasters.GDALsource} of Dict{String, Any} with 4 entries:
  "units"    => ""
  "offset"   => 0.0
  "filepath" => "/vsimem/tmp"
  "scale"    => 1.0
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ raster โ”ค
  extent: Extent(X = (995456.5461976258, 2.110456546197626e6), Y = (4.741960927303582e6, 6.191960927303582e6))
  missingval: NaN32
  crs: PROJCS["unknown",GEOGCS["unknown",DATUM["Unknown based on GRS80 ellipsoid",SPHEROID["GRS 1980",6378137,298.257222101004,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",173],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",1600000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]
โ””โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”˜
 โ†“ โ†’          6.18696e6    6.18196e6  โ€ฆ    4.74696e6    4.74196e6
 9.95457e5  NaN          NaN             NaN          NaN
 โ‹ฎ                                    โ‹ฑ                 โ‹ฎ
 2.10546e6  NaN          NaN             NaN          NaN

The resulting raster r and the lines layer coastline are plotted in Figure 5.9. Note that the raster values are average elevations based on \(5 \times 5\) pixels, but this is irrelevant for the subsequent calculation; the raster is going to be used as a template, and all of its values will be replaced with distances to coastline (Figure 5.10).

fig, ax, plt = plot(r)
lines!(ax, coastline; color = :red)
fig
Figure 5.9: Template to calculate distance to nearest geometry (coastlines, in red)

To calculate the actual distances, we must convert each pixel to a vector (point) geometry. For this purpose, we use the technique demonstrated in Section 5.5.2, but simply select the pixels that are not missing.

dp = DimPoints(r)
nonmissing_points = dp[r .=== missingval(r)]
53316-element Vector{Tuple{Float64, Float64}}:
 (995456.5461976258, 6.186960927303582e6)
 (1.0004565461976258e6, 6.186960927303582e6)
 (1.0054565461976257e6, 6.186960927303582e6)
 (1.0104565461976259e6, 6.186960927303582e6)
 (1.0154565461976258e6, 6.186960927303582e6)
 (1.0204565461976258e6, 6.186960927303582e6)
 (1.0254565461976258e6, 6.186960927303582e6)
 (1.0304565461976257e6, 6.186960927303582e6)
 (1.0354565461976259e6, 6.186960927303582e6)
 (1.0404565461976258e6, 6.186960927303582e6)
 โ‹ฎ
 (2.0654565461976258e6, 4.741960927303582e6)
 (2.0704565461976258e6, 4.741960927303582e6)
 (2.0754565461976258e6, 4.741960927303582e6)
 (2.0804565461976258e6, 4.741960927303582e6)
 (2.0854565461976258e6, 4.741960927303582e6)
 (2.0904565461976258e6, 4.741960927303582e6)
 (2.0954565461976258e6, 4.741960927303582e6)
 (2.100456546197626e6, 4.741960927303582e6)
 (2.105456546197626e6, 4.741960927303582e6)

The result is a vector of 2-tuples, which are recognized as GeoInterface point geometries.

We can compute the Cartesian distance from each point to the nearest line in the coastline multilinestring using the distance method from GeometryOps.

distances = GO.distance.((coastline,), nonmissing_points)
53316-element Vector{Float64}:
 572764.2876244619
 567764.3239539921
 562764.3609290748
 557764.3985670706
 552764.436885969
 547764.4759044152
 542764.5156417422
 537764.5561180016
 532764.5973539978
 527764.639371324
      โ‹ฎ
 610809.9189554364
 614562.8629679172
 618333.460184955
 622121.3896252689
 625926.3366197749
 629747.9927075815
 633586.0555316888
 637440.2287345697
 641310.2218537896

Finally, we rasterize (see Section 5.4.1) the distances into our raster template.

img = Rasters.rasterize(
    last,
    nonmissing_points;
    to = r,
    fill = distances,
)
โ•ญโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ•ฎ
โ”‚ 223ร—290 Raster{Union{Missing, Float64},2} last โ”‚
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”ดโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ dims โ”
  โ†“ X Projected{Float64} LinRange{Float64}(995456.5461976258, 2.105456546197626e6, 223) ForwardOrdered Regular Intervals{Start},
  โ†’ Y Projected{Float64} LinRange{Float64}(6.186960927303582e6, 4.741960927303582e6, 290) ReverseOrdered Regular Intervals{Start}
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ metadata โ”ค
  Metadata of Dict{Any, Any}()
โ”œโ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ raster โ”ค
  extent: Extent(X = (995456.5461976258, 2.110456546197626e6), Y = (4.741960927303582e6, 6.191960927303582e6))
  missingval: missing
  crs: PROJCS["unknown",GEOGCS["unknown",DATUM["Unknown based on GRS80 ellipsoid",SPHEROID["GRS 1980",6378137,298.257222101004,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",173],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",1600000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]
โ””โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”˜
 โ†“ โ†’             6.18696e6  6.18196e6  โ€ฆ       4.74696e6       4.74196e6
 9.95457e5       5.72764e5  5.72767e5     164188.0             1.67786e5
 1.00046e6       5.67764e5  5.67767e5          1.60717e5  164391.0
 1.00546e6       5.62764e5  5.62767e5          1.57329e5  161080.0
 โ‹ฎ                                     โ‹ฑ                       โ‹ฎ
 2.09546e6       3.31952e5  3.29282e5          6.30403e5       6.33586e5
 2.10046e6       3.36173e5  3.33536e5          6.34276e5       6.3744e5
 2.10546e6  340415.0        3.37811e5          6.38166e5       6.4131e5

The final result, a raster of distances to the nearest coastline, is shown in Figure 5.10.

fig, ax, plt = plot(img)
lines!(ax, coastline; color = :red)
Colorbar(fig[1, 2], plt; label = "Distance to coastline (m)")
fig
Figure 5.10: Distance to nearest coastline in New Zealand