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

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{Union{Missing, Float32}, 2} ┐
├──────────────────────────────────────────────┴───────────────────────── dims ┐
  ↓ X Projected{Float64} 995456.5461976258:1000.0:2.109456546197626e6 ForwardOrdered Regular Intervals{Start},
  → Y Projected{Float64} 6.190960927303582e6:-1000.0:4.741960927303582e6 ReverseOrdered Regular Intervals{Start}
├──────────────────────────────────────────────────────────────────── metadata ┤
  Metadata{Rasters.GDALsource} of Dict{String, Any} with 1 entry:
  "filepath" => "data/nz_elev.tif"
├────────────────────────────────────────────────────────────────────── raster ┤
  missingval: missing
  extent: Extent(X = (995456.5461976258, 2.110456546197626e6), Y = (4.741960927303582e6, 6.191960927303582e6))
  crs: PROJCS["unknown",GEOGCS["unknown",DATUM["Unknown based on GRS80 elli...
└──────────────────────────────────────────────────────────────────────────────┘
 ↓ →        6.19096e6  6.18996e6  …  4.74396e6  4.74296e6  4.74196e6
 9.95457e5   missing    missing       missing    missing    missing
 9.96457e5   missing    missing       missing    missing    missing
 ⋮                                ⋱                        ⋮
 2.10846e6   missing    missing       missing    missing    missing
 2.10946e6   missing    missing   …   missing    missing    missing

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), always_xy=true)

1×12 DataFrame

Row	geometry	UNIT_CODE	GIS_Notes	UNIT_NAME	DATE_EDIT	STATE	REGION	GNIS_ID	UNIT_TYPE	CREATED_BY	METADATA	PARKNAME
	IGeometr…	String	String	String	DateTime	String	String	String	String	String	String	String
1	Geometry: wkbPolygon	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.

DataFrame

out_image_mask = Rasters.mask(src_srtm; with = zion)

┌ 465×457 Raster{Union{Missing, UInt16}, 2} ┐
├───────────────────────────────────────────┴──────────────────────────── dims ┐
  ↓ X Projected{Float64} -113.23958321278403:0.000833333333277784:-112.85291654614313 ForwardOrdered Regular Intervals{Start},
  → Y Projected{Float64} 37.51208342983253:-0.0008333333332777888:37.13208342985786 ReverseOrdered Regular Intervals{Start}
├────────────────────────────────────────────────────────────────────── raster ┤
  missingval: missing
  extent: Extent(X = (-113.23958321278403, -112.85208321280986), Y = (37.13208342985786, 37.512916763165805))
  crs: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.25722...
└──────────────────────────────────────────────────────────────────────────────┘
    ↓ →    37.5121    37.5113    …  37.1338    37.1329    37.1321
 -113.24     missing    missing       missing    missing    missing
 -113.239    missing    missing       missing    missing    missing
 -113.238    missing    missing       missing    missing    missing
 -113.237    missing    missing       missing    missing    missing
    ⋮                            ⋱   ⋮                    
 -112.855    missing    missing       missing    missing    missing
 -112.855    missing    missing       missing    missing    missing
 -112.854    missing    missing       missing    missing    missing
 -112.853    missing    missing  …    missing    missing    missing

Single geometry

masker = zion.geometry[1]

Geometry: POLYGON ((-113.084386842835 37.1713038811739,-113. ... 739))

Rasters.mask(src_srtm; with = masker)

┌ 465×457 Raster{Union{Missing, UInt16}, 2} ┐
├───────────────────────────────────────────┴──────────────────────────── dims ┐
  ↓ X Projected{Float64} -113.23958321278403:0.000833333333277784:-112.85291654614313 ForwardOrdered Regular Intervals{Start},
  → Y Projected{Float64} 37.51208342983253:-0.0008333333332777888:37.13208342985786 ReverseOrdered Regular Intervals{Start}
├────────────────────────────────────────────────────────────────────── raster ┤
  missingval: missing
  extent: Extent(X = (-113.23958321278403, -112.85208321280986), Y = (37.13208342985786, 37.512916763165805))
  crs: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.25722...
└──────────────────────────────────────────────────────────────────────────────┘
    ↓ →    37.5121    37.5113    …  37.1338    37.1329    37.1321
 -113.24     missing    missing       missing    missing    missing
 -113.239    missing    missing       missing    missing    missing
 -113.238    missing    missing       missing    missing    missing
 -113.237    missing    missing       missing    missing    missing
    ⋮                            ⋱   ⋮                    
 -112.855    missing    missing       missing    missing    missing
 -112.855    missing    missing       missing    missing    missing
 -112.854    missing    missing       missing    missing    missing
 -112.853    missing    missing  …    missing    missing    missing

Vector of geometries

masker = zion.geometry

1-element GeoDataFrames.GeometryVector{ArchGDAL.IGeometry{ArchGDAL.wkbPolygon}}:
 Geometry: POLYGON ((-113.084386842835 37.1713038811739,-113. ... 739))

Rasters.mask(src_srtm; with = masker)

┌ 465×457 Raster{Union{Missing, UInt16}, 2} ┐
├───────────────────────────────────────────┴──────────────────────────── dims ┐
  ↓ X Projected{Float64} -113.23958321278403:0.000833333333277784:-112.85291654614313 ForwardOrdered Regular Intervals{Start},
  → Y Projected{Float64} 37.51208342983253:-0.0008333333332777888:37.13208342985786 ReverseOrdered Regular Intervals{Start}
├────────────────────────────────────────────────────────────────────── raster ┤
  missingval: missing
  extent: Extent(X = (-113.23958321278403, -112.85208321280986), Y = (37.13208342985786, 37.512916763165805))
  crs: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.25722...
└──────────────────────────────────────────────────────────────────────────────┘
    ↓ →    37.5121    37.5113    …  37.1338    37.1329    37.1321
 -113.24     missing    missing       missing    missing    missing
 -113.239    missing    missing       missing    missing    missing
 -113.238    missing    missing       missing    missing    missing
 -113.237    missing    missing       missing    missing    missing
    ⋮                            ⋱   ⋮                    
 -112.855    missing    missing       missing    missing    missing
 -112.855    missing    missing       missing    missing    missing
 -112.854    missing    missing       missing    missing    missing
 -112.853    missing    missing  …    missing    missing    missing

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
)

[ Info: `missingval` set to 65535 on disk

"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.geometry[1], touches = true)

┌ 439×436 Raster{Union{Missing, UInt16}, 2} ┐
├───────────────────────────────────────────┴──────────────────────────── dims ┐
  ↓ X Projected{Float64} -113.22874987945141:0.000833333333277784:-112.86374987947575 ForwardOrdered Regular Intervals{Start},
  → Y Projected{Float64} 37.50375009649975:-0.0008333333332777888:37.14125009652391 ReverseOrdered Regular Intervals{Start}
├──────────────────────────────────────────────────────────────────── metadata ┤
  Metadata{Rasters.GDALsource} of Dict{String, Any} with 1 entry:
  "filepath" => "data/srtm.tif"
├────────────────────────────────────────────────────────────────────── raster ┤
  missingval: missing
  extent: Extent(X = (-113.22874987945141, -112.86291654614247), Y = (37.14125009652391, 37.504583429833026))
  crs: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.25722...
└──────────────────────────────────────────────────────────────────────────────┘
    ↓ →        37.5038      37.5029      37.5021  …      37.1421      37.1413
 -113.229  0x0668       0x065f       0x0659          0x05a0       0x05df
 -113.228  0x0663       0x065e       0x0658          0x059e       0x05e4
    ⋮                                             ⋱       ⋮       
 -112.865  0x09ef       0x0a02       0x0a0b          0x06d9       0x06d4
 -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{Union{Missing, UInt16}, 2} ┐
├───────────────────────────────────────────┴──────────────────────────── dims ┐
  ↓ X Projected{Float64} -113.22874987945141:0.000833333333277784:-112.86374987947575 ForwardOrdered Regular Intervals{Start},
  → Y Projected{Float64} 37.50375009649975:-0.0008333333332777888:37.14125009652391 ReverseOrdered Regular Intervals{Start}
├────────────────────────────────────────────────────────────────────── raster ┤
  missingval: missing
  extent: Extent(X = (-113.22874987945141, -112.86291654614247), Y = (37.14125009652391, 37.504583429833026))
  crs: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.25722...
└──────────────────────────────────────────────────────────────────────────────┘
    ↓ →    37.5038    37.5029    …  37.1429    37.1421    37.1413
 -113.229    missing    missing       missing    missing    missing
 -113.228    missing    missing       missing    missing    missing
 -113.227    missing    missing       missing    missing    missing
 -113.226    missing    missing       missing    missing    missing
    ⋮                            ⋱              ⋮         
 -112.866    missing    missing       missing    missing    missing
 -112.865    missing    missing       missing    missing    missing
 -112.865    missing    missing       missing    missing    missing
 -112.864    missing    missing       missing    missing    missing

and we write it to file using Rasters.write:

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

[ Info: `missingval` set to 65535 on disk

"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.geometry; color = :transparent, strokecolor = :black, strokewidth = 0.75)

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

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

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

display(fig);

Figure 5.1: Raster masking and cropping

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.

• To points (Section 5.3.1) or to lines (Section 5.3.2), via the Rasters.extract function
• To polygons (Section 5.3.3), via the Rasters.zonal function

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.geometry, color=:black, strokecolor=:white, strokewidth = 1);
display(fig);

Figure 5.2: 30 point locations within the Zion National Park, with elevation in the background

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

24 rows omitted

Row
	UInt16?
1	1802
2	2433
3	1886
⋮	⋮
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

24 rows omitted

Row	geometry	elev1
	IGeometr…	UInt16?
1	Geometry: wkbPoint	1802
2	Geometry: wkbPoint	2433
3	Geometry: wkbPoint	1886
⋮	⋮	⋮
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 Section 1.2 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}([[-113.2, 37.45], [-112.9, 37.2]])

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), always_xy=true)

GeoInterface.Wrappers.LineString{false, false}([(305399.67208180577, 4.147066650206682e6), (331380.8917453843, 4.1187500947884847e6)], crs = "EPSG: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}([(305399.67208180577, 4.147066650206682e6), … (249) … , (331380.8917453843, 4.1187500947884847e6)], crs = "EPSG: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), always_xy=true)

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

245 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
⋮	⋮	⋮	⋮
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.geometry; 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

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.

TipYou can pass your own custom function too!

# 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

TipYou can pass multiple functions to a single Rasters.zonal call

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), always_xy=true))
using StatsBase
counts = StatsBase.countmap(out_image)

Dict{Union{Missing, UInt8}, Int64} with 8 entries:
  0x05    => 203701
  0x04    => 298299
  0x06    => 235
  0x07    => 62
  0x02    => 4205
  missing => 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.geometry; 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.geometry; source = GI.crs(zion.geometry[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

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), always_xy=true)

540×6 DataFrame

534 rows omitted

Row	geometry	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
⋮	⋮	⋮	⋮	⋮	⋮	⋮
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 Sampled{Float64} 523038.6145227547:15.684723604718238:538707.6534038682 ForwardOrdered Regular Intervals{Start},
  → Y Sampled{Float64} 174934.65727249545:10.036751270483991:184961.37179170895 ForwardOrdered Regular Intervals{Start}
├────────────────────────────────────────────────────────────────────── raster ┤
  missingval: missing
  extent: Extent(X = (523038.6145227547, 538723.3381274729), Y = (174934.65727249545, 184971.40854297945))
└──────────────────────────────────────────────────────────────────────────────┘
      ↓ →        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.23086e5   missing    missing       missing    missing    missing
      ⋮                                ⋱                        ⋮
      5.38661e5   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 Sampled{Float64} 523038.6145227547:15.684723604718238:538707.6534038682 ForwardOrdered Regular Intervals{Start},
  → Y Sampled{Float64} 174934.65727249545:10.036751270483991:184961.37179170895 ForwardOrdered Regular Intervals{Start}
├────────────────────────────────────────────────────────────────────── raster ┤
  missingval: missing
  extent: Extent(X = (523038.6145227547, 538723.3381274729), Y = (174934.65727249545, 184971.40854297945))
└──────────────────────────────────────────────────────────────────────────────┘
      ↓ →        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.23086e5   missing    missing       missing    missing    missing
      ⋮                                ⋱                        ⋮
      5.38661e5   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 Sampled{Float64} 523038.6145227547:15.684723604718238:538707.6534038682 ForwardOrdered Regular Intervals{Start},
  → Y Sampled{Float64} 174934.65727249545:10.036751270483991:184961.37179170895 ForwardOrdered Regular Intervals{Start}
├────────────────────────────────────────────────────────────────────── raster ┤
  missingval: missing
  extent: Extent(X = (523038.6145227547, 538723.3381274729), Y = (174934.65727249545, 184971.40854297945))
└──────────────────────────────────────────────────────────────────────────────┘
      ↓ →        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.23086e5   missing    missing       missing    missing    missing
      ⋮                                ⋱                        ⋮
      5.38661e5   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.geometry; # 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.6145227547:15.684723604718238:538707.6534038682 ForwardOrdered Regular Intervals{Start},
  → Y Projected{Float64} 174934.65727249545:10.036751270483991:184961.37179170895 ForwardOrdered Regular Intervals{Start}
├────────────────────────────────────────────────────────────────────── raster ┤
  missingval: missing
  extent: Extent(X = (523038.6145227547, 538723.3381274729), Y = (174934.65727249545, 184971.40854297945))
  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.23086e5   missing    missing       missing    missing    missing
      ⋮                                ⋱                        ⋮
      5.38661e5   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, :geometry])
f, a, p = scatter(nonmissing_df.geometry; 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	geometry	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.geometry)
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}([GeoInterface.Wrappers.LineString(Geometry: LINEARRING (-118.603375 33.478098,-118.368301 33.4 ... 8098)), … (3) … , GeoInterface.Wrappers.LineString(Geometry: LINEARRING (-124.211605 41.99846,-123.347562 41.99 ... 9846))], crs = "WellKnownText with CRS mode: 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 Sampled{Float64} -124.409591:0.5:-114.409591 ForwardOrdered Regular Intervals{Start},
  → Y Sampled{Float64} 32.534156:0.5:41.534156 ForwardOrdered Regular Intervals{Start}
├────────────────────────────────────────────────────────────────────── raster ┤
  missingval: missing
  extent: Extent(X = (-124.409591, -113.909591), Y = (32.534156, 42.034156))
└──────────────────────────────────────────────────────────────────────────────┘
    ↓ →   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
 -122.91    missing    missing    missing       missing    missing   1
    ⋮                                      ⋱                        
 -115.91   1           missing    missing       missing    missing    missing
 -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;
    fill = 1,
    res = 0.5,
    boundary = :center,
)

┌ 21×19 Raster{Union{Missing, Int64}, 2} last ┐
├─────────────────────────────────────────────┴────────────────────────── dims ┐
  ↓ X Sampled{Float64} -124.409591:0.5:-114.409591 ForwardOrdered Regular Intervals{Start},
  → Y Sampled{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 ┤
  missingval: missing
  extent: Extent(X = (-124.409591, -113.909591), Y = (32.534156, 42.034156))
└──────────────────────────────────────────────────────────────────────────────┘
    ↓ →   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
    ⋮                                      ⋱                        
 -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.

DimPoints

dp = DimPoints(california_raster1)

┌ 21×19 DimPoints{Tuple{Float64, Float64}, 2} ┐
├─────────────────────────────────────────────┴────────────────────────── dims ┐
  ↓ X Sampled{Float64} -124.409591:0.5:-114.409591 ForwardOrdered Regular Intervals{Start},
  → Y Sampled{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)
    ⋮                                                 ⋱  
 -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)

Constructing from lookups

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 DimArray{Tuple{Float64, Float64}, 2} ┐
├────────────────────────────────────────────┴─────────────────────────── dims ┐
  ↓ X Sampled{Float64} -124.409591:0.5:-114.409591 ForwardOrdered Regular Intervals{Start},
  → Y Sampled{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)
    ⋮                                                 ⋱  
 -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)

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)
display(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)

FeatureCollection([Feature(GeoInterface.Wrappers.MultiPolygon{false, false}([GeoInterface.Wrappers.Polygon([GeoInterface.Wrappers.LinearRing([(1.5, 1.5), … (5) … , (1.5, 1.5)], extent = Extent(X = (0.5, 2.5), Y = (1.5, 4.5)), crs = "WellKnownText with CRS mode: 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 = Extent(X = (0.5, 2.5), Y = (1.5, 4.5)), crs = "WellKnownText with CRS mode: 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) … , GeoInterface.Wrappers.Polygon([GeoInterface.Wrappers.LinearRing([(4.5, 3.5), … (8) … , (4.5, 3.5)], extent = Extent(X = (2.5, 5.5), Y = (1.5, 3.5)), crs = "WellKnownText with CRS mode: 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 = Extent(X = (2.5, 5.5), Y = (1.5, 3.5)), crs = "WellKnownText with CRS mode: 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 = Extent(X = (0.5, 5.5), Y = (0.5, 6.5)), crs = "WellKnownText with CRS mode: 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"]]"), properties = (value = 0x00), extent = Extent(X = (0.5, 5.5), Y = (0.5, 6.5)), crs = "WellKnownText with CRS mode: 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"]]"), Feature(GeoInterface.Wrappers.MultiPolygon{false, false}([GeoInterface.Wrappers.Polygon([GeoInterface.Wrappers.LinearRing([(0.5, 4.5), … (13) … , (0.5, 4.5)], extent = Extent(X = (0.5, 6.5), Y = (3.5, 6.5)), crs = "WellKnownText with CRS mode: 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 = Extent(X = (0.5, 6.5), Y = (3.5, 6.5)), crs = "WellKnownText with CRS mode: 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) … , GeoInterface.Wrappers.Polygon([GeoInterface.Wrappers.LinearRing([(2.5, 0.5), … (3) … , (2.5, 0.5)], extent = Extent(X = (2.5, 3.5), Y = (0.5, 1.5)), crs = "WellKnownText with CRS mode: 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 = Extent(X = (2.5, 3.5), Y = (0.5, 1.5)), crs = "WellKnownText with CRS mode: 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 = Extent(X = (0.5, 6.5), Y = (0.5, 6.5)), crs = "WellKnownText with CRS mode: 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"]]"), properties = (value = 0x01), extent = Extent(X = (0.5, 6.5), Y = (0.5, 6.5)), crs = "WellKnownText with CRS mode: 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"]]"), Feature(GeoInterface.Wrappers.MultiPolygon{false, false}([GeoInterface.Wrappers.Polygon([GeoInterface.Wrappers.LinearRing([(1.5, 1.5), … (5) … , (1.5, 1.5)], extent = Extent(X = (1.5, 3.5), Y = (1.5, 3.5)), crs = "WellKnownText with CRS mode: 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 = Extent(X = (1.5, 3.5), Y = (1.5, 3.5)), crs = "WellKnownText with CRS mode: 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"]]"), … (4) … , GeoInterface.Wrappers.Polygon([GeoInterface.Wrappers.LinearRing([(1.5, 6.5), … (3) … , (1.5, 6.5)], extent = Extent(X = (0.5, 1.5), Y = (5.5, 6.5)), crs = "WellKnownText with CRS mode: 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 = Extent(X = (0.5, 1.5), Y = (5.5, 6.5)), crs = "WellKnownText with CRS mode: 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 = Extent(X = (0.5, 6.5), Y = (0.5, 6.5)), crs = "WellKnownText with CRS mode: 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"]]"), properties = (value = 0x02), extent = Extent(X = (0.5, 6.5), Y = (0.5, 6.5)), crs = "WellKnownText with CRS mode: 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"]]")], crs = "WellKnownText with CRS mode: 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 = Extent(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	GeoInterface.Wrappers.MultiPolygon([GeoInterface.Wrappers.Polygon([GeoInterface.Wrappers.LinearRing([(1.5,1.5),…(5)…,(1.5,1.5)])]),…(2)…,GeoInterface.Wrappers.Polygon([GeoInterface.Wrappers.LinearRing([(4.5,3.5),…(8)…,(4.5,3.5)])])])
2	1	GeoInterface.Wrappers.MultiPolygon([GeoInterface.Wrappers.Polygon([GeoInterface.Wrappers.LinearRing([(0.5,4.5),…(13)…,(0.5,4.5)])]),…(2)…,GeoInterface.Wrappers.Polygon([GeoInterface.Wrappers.LinearRing([(2.5,0.5),…(3)…,(2.5,0.5)])])])
3	2	GeoInterface.Wrappers.MultiPolygon([GeoInterface.Wrappers.Polygon([GeoInterface.Wrappers.LinearRing([(1.5,1.5),…(5)…,(1.5,1.5)])]),…(4)…,GeoInterface.Wrappers.Polygon([GeoInterface.Wrappers.LinearRing([(1.5,6.5),…(3)…,(1.5,6.5)])])])

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

30 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
⋮	⋮	⋮
34	(0.0, -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

30 rows omitted

Row	XY	layer1
	IGeometr…	Int16
1	Geometry: wkbPoint	1
2	Geometry: wkbPoint	2
3	Geometry: wkbPoint	3
⋮	⋮	⋮
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)
display(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.geometry)))) .|> 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), always_xy=true)

GeoInterface.Wrappers.MultiLineString{false, false}([GeoInterface.Wrappers.LineString([(1.2299979013e6, 4.7934322369e6), … (24) … , (1.2299979013e6, 4.7934322369e6)], crs = "WellKnownText with CRS mode: 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) … , GeoInterface.Wrappers.LineString([(1.8266292121e6, 5.9745333354e6), … (11) … , (1.8266292121e6, 5.9745333354e6)], crs = "WellKnownText with CRS mode: 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]]")], crs = "WellKnownText with CRS mode: 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{Union{Missing, Float32}, 2} ┐
├────────────────────────────────────────────┴─────────────────────────── dims ┐
  ↓ X Projected{Float64} 995456.5461976258:5000.0:2.105456546197626e6 ForwardOrdered Regular Intervals{Start},
  → Y Projected{Float64} 6.186960927303582e6:-5000.0:4.741960927303582e6 ReverseOrdered Regular Intervals{Start}
├──────────────────────────────────────────────────────────────────── metadata ┤
  Metadata{Rasters.GDALsource} of Dict{String, Any} with 1 entry:
  "filepath" => "/vsimem/tmp"
├────────────────────────────────────────────────────────────────────── raster ┤
  missingval: missing
  extent: Extent(X = (995456.5461976258, 2.110456546197626e6), Y = (4.741960927303582e6, 6.191960927303582e6))
  crs: PROJCS["unknown",GEOGCS["unknown",DATUM["Unknown based on GRS80 elli...
└──────────────────────────────────────────────────────────────────────────────┘
 ↓ →        6.18696e6  6.18196e6  …  4.75196e6  4.74696e6  4.74196e6
 9.95457e5   missing    missing       missing    missing    missing
 1.00046e6   missing    missing       missing    missing    missing
 ⋮                                ⋱                        ⋮
 2.10046e6   missing    missing       missing    missing    missing
 2.10546e6   missing    missing       missing    missing    missing

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)
display(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.0054565461976258e6, 6.186960927303582e6)
 (1.0104565461976258e6, 6.186960927303582e6)
 (1.0154565461976258e6, 6.186960927303582e6)
 (1.0204565461976258e6, 6.186960927303582e6)
 (1.0254565461976258e6, 6.186960927303582e6)
 (1.0304565461976258e6, 6.186960927303582e6)
 (1.0354565461976258e6, 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.3609290747
 557764.3985670707
 552764.436885969
 547764.4759044152
 542764.5156417422
 537764.5561180015
 532764.5973539979
 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} 995456.5461976258:5000.0:2.105456546197626e6 ForwardOrdered Regular Intervals{Start},
  → Y Projected{Float64} 6.186960927303582e6:-5000.0:4.741960927303582e6 ReverseOrdered Regular Intervals{Start}
├────────────────────────────────────────────────────────────────────── raster ┤
  missingval: missing
  extent: Extent(X = (995456.5461976258, 2.110456546197626e6), Y = (4.741960927303582e6, 6.191960927303582e6))
  crs: PROJCS["unknown",GEOGCS["unknown",DATUM["Unknown based on GRS80 elli...
└──────────────────────────────────────────────────────────────────────────────┘
 ↓ →             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
 1.01046e6       5.57764e5  5.57767e5          1.54029e5  157858.0
 ⋮                                     ⋱                       ⋮
 2.09046e6       3.27753e5  3.25048e5  …       6.26545e5  629748.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)")
display(fig);

Figure 5.10: Distance to nearest coastline in New Zealand