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
6 Raster-vector interactions
Prerequisites
This chapter requires importing the following packages:
set_theme!(
Makie.= (; rasterize = 2),
Heatmap = (; rasterize = 2),
Surface )
It also relies on the following data files:
= Raster("data/srtm.tif")
src_srtm = Raster("data/nlcd.tif")
src_nlcd = Raster("output/grain.tif")
src_grain = Raster("output/elev.tif")
src_elev = Raster("data/dem.tif")
src_dem = GeoDataFrames.read("data/zion.gpkg")
zion = GeoDataFrames.read("data/zion_points.gpkg")
zion_points = GeoDataFrames.read("data/cycle_hire_osm.gpkg")
cycle_hire_osm = GeoDataFrames.read("data/us_states.gpkg")
us_states = GeoDataFrames.read("data/nz.gpkg")
nz = Raster("data/nz_elev.tif") src_nz_elev
โญโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโฎ
โ 1115ร1450 Raster{Union{Missing, 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 1 entry:
"filepath" => "data/nz_elev.tif"
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโ 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.19096e6 6.18996e6 โฆ 4.74396e6 4.74296e6 4.74196e6
9.95457e5 missing missing missing missing missing
9.96457e5 missing missing missing missing missing
โฎ โฑ โฎ
2.10746e6 missing missing missing missing missing
2.10846e6 missing missing missing missing missing
2.10946e6 missing missing โฆ missing missing missing
6.1 Introduction
This chapter focuses on interactions between raster and vector geographic data models, both introduced in Chapter 2. It includes four main techniques:
- Raster cropping and masking using vector objects (Section 6.2)
- Extracting raster values using different types of vector data (Section Section 6.3)
- Raster-vector conversion (Section 6.4 and Section 6.5)
These concepts are demonstrated using data from in previous chapters, to understand their potential real-world applications.
6.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 6.1 (b) and (c), and the following two examples, illustrate the difference between masking and cropping). Both operations reduce object memory use and associated computational resources for subsequent analysis steps, 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 namedsrc_srtm
(see Figure 6.1 (a)) - The
zion.gpkg
vector layer representing the Zion National Park boundaries (aGeoDataFrame
namedzion
)
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 CRS of the raster src_srtm
.
= GO.reproject(zion; target_crs = GI.crs(src_srtm)) zion
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!
= Rasters.mask(src_srtm; with = zion) out_image_mask
โญโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโฎ
โ 465ร457 Raster{Union{Missing, 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}
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโ raster โค
extent: Extent(X = (-113.23958321278403, -112.85208321280986), Y = (37.13208342985786, 37.512916763165805))
missingval: missing
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.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.856 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 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:
= Rasters.mask(src_srtm; with = zion, missingval = 9999) out_image_mask
We can write this masked raster to file with Rasters.write
, as usual:
write("output/srtm_masked.tif", out_image_mask; force = true
Rasters. )
[ 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.
= Rasters.crop(src_srtm; to = zion, touches = true) out_image_crop
โญโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโฎ
โ 439ร436 Raster{Union{Missing, 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 1 entry:
"filepath" => "data/srtm.tif"
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโ raster โค
extent: Extent(X = (-113.22874987945141, -112.86291654614247), Y = (37.14125009652391, 37.504583429833026))
missingval: missing
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
-113.228 0x0663 0x065e 0x0658 0x059e 0x05e4
โฎ โฑ โฎ
-112.865 0x09db 0x09e2 0x09e7 0x06d4 0x06d0
-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:
= Rasters.crop(out_image_mask; to = zion, touches = true) out_image_mask_crop
โญโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโฎ
โ 439ร436 Raster{Union{Missing, 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}
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโ raster โค
extent: Extent(X = (-113.22874987945141, -112.86291654614247), Y = (37.14125009652391, 37.504583429833026))
missingval: missing
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.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.867 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 as usual:
write("output/srtm_masked_cropped.tif", out_image_mask_crop; force = true) Rasters.
[ Info: `missingval` set to 65535 on disk
"output/srtm_masked_cropped.tif"
Figure 6.1 shows the original raster, and the three masking and/or cropping results.
= Figure(size = (600, 600))
fig
= Axis(fig[1, 1]; title = "Original")
ax1 plot!(ax1, src_srtm)
poly!(ax1, zion.geom; color = :transparent, strokecolor = :black, strokewidth = 0.75)
= Axis(fig[1, 2]; title = "Masked")
ax2 plot!(ax2, out_image_mask)
poly!(ax2, zion.geom; color = :transparent, strokecolor = :black, strokewidth = 0.75)
= Axis(fig[2, 1]; title = "Cropped")
ax3 plot!(ax3, out_image_crop)
poly!(ax3, zion.geom; color = :transparent, strokecolor = :black, strokewidth = 0.75)
= Axis(fig[2, 2]; title = "Masked+Cropped")
ax4 plot!(ax4, out_image_mask_crop)
poly!(ax4, zion.geom; color = :transparent, strokecolor = :black, strokewidth = 0.75)
fig
6.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 6.4.
Rasters.jl provides modular raster extraction and statistics functions, and we use this package in the following examples.
- To points (Section 6.3.1) or to lines (Section 6.3.2), via the
Rasters.extract
function - To polygons (Section 6.3.3), via the
Rasters.zonal
function
6.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 6.2).
= plot(src_srtm)
fig, ax, plt scatter!(ax, zion_points.geom, color=:black, strokecolor=:white, strokewidth = 1);
fig
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.
= Rasters.extract(src_srtm, zion_points; geometry = false) .|> first result1
30-element Vector{UInt16}:
0x070a
0x0981
0x075e
0x055a
0x05ac
0x0663
0x0564
0x07f0
0x0726
0x0744
โฎ
0x0711
0x0722
0x060e
0x0707
0x0836
0x0846
0x055c
0x0771
0x0626
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.
Rasters.jl does not yet support interpolation, 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.
1:5] result1[
5-element Vector{UInt16}:
0x070a
0x0981
0x075e
0x055a
0x05ac
To get a GeoDataFrame
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.
"elev1"] = result1
zion_points[!, zion_points
Row | geom | 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
6.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 6.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.
= [[-113.2, 37.45], [-112.9, 37.2]]
coords = GI.LineString(coords) zion_transect
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 due to long climbs.
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 UTM.
= GO.reproject(zion_transect; target_crs = GFT.EPSG(32612), source_crs = GFT.EPSG(4326)) zion_transect_utm
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
.
We first compute the length of the line, and then use this to segmentize the line into approximately 250 points.
= GO.centroid_and_length(zion_transect_utm)
_centroid, linelen = GO.segmentize(zion_transect_utm; max_distance = linelen / 250) zion_transect_line
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 by using GI.getpoint
on the line, and then reproject the points to the CRS of the raster.
= GO.reproject(GI.getpoint(zion_transect_line); target_crs = GI.crs(src_srtm), source_crs = GI.crs(zion_transect_line)) zion_transect_pnt
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 GeoDataFrame
, to accommodate extra attributes), using the point extraction method shown earlier (Section 6.3.1). We also attach the respective distance cutoff points distances
.
= DataFrame(geometry = zion_transect_pnt)
zion_transect_pnt
= Rasters.extract(src_srtm, zion_transect_pnt; geometry = false) .|> first
result "elev"] = result zion_transect_pnt[!,
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
= GO.distance.(zion_transect_pnt.geometry[1:end-1], zion_transect_pnt.geometry[2:end])
distances_between_points # Compute cumulative distances along the line
= cumsum([0.0; distances_between_points])
cumulative_distances "dist"] = cumulative_distances
zion_transect_pnt[!, zion_transect_pnt
Row | geometry | elev | dist |
---|---|---|---|
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 6.3.
# Raster and a line transect
= plot(src_srtm)
fig, ax, plt lines!(ax, zion_transect; color = :black)
poly!(ax, zion.geom; color = :transparent, strokecolor = :white, strokewidth = 0.75)
display(fig)
# Elevation profile
= lines(
fig, ax, plt
zion_transect_pnt.dist,
zion_transect_pnt.elev;= (;
axis = "Distance (m)",
xlabel = "Elevation (m)",
ylabel
)
) fig
6.3.3 Extraction to polygons
The final type of geographic vector object for raster extraction is polygons. Like lines, polygons tend to return many raster values per vector geometry. For continuous rasters (Figure 6.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 (in this case a list of length 1, since there is just one polygon).
Rasters.zonal does not allow passing multiple functions. Would this provide a speedup? Do we care?
using Statistics # for `mean`
= Rasters.zonal(mean, src_srtm; of = zion)
rmean = Rasters.zonal(minimum, src_srtm; of = zion)
rmin = Rasters.zonal(maximum, src_srtm; of = zion)
rmax = (rmean, rmin, rmax) result
(Union{Missing, Float64}[1818.211830154405], Union{Missing, UInt16}[0x0462], Union{Missing, UInt16}[0x0a65])
Itโs straightforward to transform the result to a DataFrame as well:
DataFrame(mean = rmean, min = rmin, max = rmax)
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.
Rasters.zonal
accepts any function that works on iterables and returns a single value, like mean
, minimum
, maximum
, std
, median
, mode
, sum
, prod
, etc.
To count occurrences of categorical raster values within polygons (Figure 6.4 (b)), we can use masking (Section 6.2) combined with StatsBase.countmap
, as follows.
= Rasters.mask(src_nlcd; with = GO.reproject(zion; target_crs = GI.crs(src_nlcd)))
out_image using StatsBase
= StatsBase.countmap(out_image) counts
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 6.4 illustrates the two types of raster extraction to polygons described above.
# Continuous raster
= plot(src_srtm)
fig, ax, plt poly!(ax, zion.geom; color = :transparent, strokecolor = :black, strokewidth = 0.75)
display(fig)
# Categorical raster
= plot(src_nlcd; colormap = cgrad(:Set3; categorical = true), source = GI.crs(src_nlcd), axis = (; type = GeoAxis, dest = GI.crs(src_nlcd)))
fig, ax, plt poly!(ax, zion.geom; source = GI.crs(zion.geom[1]), color = :transparent, strokecolor = :black, strokewidth = 0.75)
= Colorbar(fig[1, 2], plt)
cm = false
ax.xgridvisible = false
ax.ygridvisible # TODO: make GeoMakie better on small extents fig
6.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 2, 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
We will use the ; data::FeatureCollection, fill::Symbol
method to rasterize the zion
polygon we have been using in this chapter.
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 geographic resolution 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 6.4.1), we typically choose how to treat multiple points: either to summarize presence/absence, point count, or summed attribute values (Figure 6.5)
- in line and polygon rasterization (Section 6.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 6.6)
6.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 6.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.
= GO.reproject(cycle_hire_osm; target_crs = GFT.EPSG(27700)) cycle_hire_osm_projected
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 |
โฎ | โฎ | โฎ | โฎ | โฎ | โฎ | โฎ |
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 6.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.
= Rasters.rasterize(
ch_raster1 # reducer
last, # data
cycle_hire_osm_projected; = 1,
fill = (1000, 1000) # specify size in "pixels"
size )
โญโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโฎ
โ 1000ร1000 Raster{Union{Missing, Int64},2} last โ
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโดโโโโโโโโโโโโโโโโโโโโโโโ dims โ
โ X Sampled{Float64} 523038.61452275474:15.684723604717:538707.653403867 ForwardOrdered Regular Intervals{Start},
โ Y Sampled{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.
= Rasters.rasterize(
ch_raster2 # reducer
sum, # data
cycle_hire_osm_projected; = 1,
fill = (1000, 1000) # specify size in "pixels"
size )
โญโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโฎ
โ 1000ร1000 Raster{Union{Missing, Int64},2} sum โ
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโดโโโโโโโโโโโโโโโโโโโโโโโโ dims โ
โ X Sampled{Float64} 523038.61452275474:15.684723604717:538707.653403867 ForwardOrdered Regular Intervals{Start},
โ Y Sampled{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
.
= Rasters.rasterize(
ch_raster3 # reducer
sum, # data
cycle_hire_osm_projected; = :capacity,
fill = (1000, 1000) # specify size in "pixels"
size )
โญโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโฎ
โ 1000ร1000 Raster{Union{Missing, Float64},2} capacity โ
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโดโโโโโโโโโโโโโโโโโ dims โ
โ X Sampled{Float64} 523038.61452275474:15.684723604717:538707.653403867 ForwardOrdered Regular Intervals{Start},
โ Y Sampled{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.
= Rasters.rasterize(
ch_raster3 # reducer
sum, # data
cycle_hire_osm_projected.geom; = cycle_hire_osm_projected.capacity,
fill = GI.crs(cycle_hire_osm_projected),
crs = (1000, 1000) # specify size in "pixels"
size )
โญโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโฎ
โ 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 6.5.
# Input points
= dropmissing(cycle_hire_osm_projected, [:capacity, :geom])
nonmissing_df = scatter(nonmissing_df.geom; color = nonmissing_df.capacity)
f, a, p 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)
6.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.
= us_states[ us_states[!, "NAME"] .== "California", :] california
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โ
= only(california.geom)
california_geom = GI.MultiLineString(GI.LineString.(GI.getexterior.(GI.getgeom(california_geom))); crs = GI.crs(california_geom)) # TODO: make this a lot better.... california_borders
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.
= Rasters.rasterize(
california_raster1
last,
california_borders;= 1,
fill = 0.5, # degrees - this is in units of GI.crs(california_borders)
res = :touches,
boundary )
โญโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโฎ
โ 20ร18 Raster{Union{Missing, Int64},2} last โ
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโดโโโโโโโโโโโโโโโโโโโโโโโโโโโ dims โ
โ X Sampled{Float64} -124.409591:0.5:-114.909591 ForwardOrdered Regular Intervals{Start},
โ Y Sampled{Float64} 32.534156:0.5:41.034156 ForwardOrdered Regular Intervals{Start}
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโ metadata โค
Metadata of Dict{Any, Any}()
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโ raster โค
extent: Extent(X = (-124.409591, -114.409591), Y = (32.534156, 41.534156))
missingval: missing
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโ
โ โ 32.5342 33.0342 33.5342 โฆ 40.0342 40.5342 41.0342
-124.41 missing missing missing 1 1 1
-123.91 missing missing missing missing missing missing
-123.41 missing missing missing missing missing missing
โฎ โฑ
-115.91 1 missing missing missing missing missing
-115.41 1 missing missing missing missing missing
-114.91 1 1 1 โฆ 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 6.6 (right).
= Rasters.rasterize(
california_raster2
last,
california;= :geom,
geometrycolumn = 1,
fill = 0.5,
res = :center,
boundary )
โญโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโฎ
โ 20ร18 Raster{Union{Missing, Int64},2} last โ
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโดโโโโโโโโโโโโโโโโโโโโโโโโโโโ dims โ
โ X Sampled{Float64} -124.409591:0.5:-114.909591 ForwardOrdered Regular Intervals{Start},
โ Y Sampled{Float64} 32.534156:0.5:41.034156 ForwardOrdered Regular Intervals{Start}
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโ metadata โค
Metadata of Dict{Any, Any} with 1 entry:
"missed_geometries" => Bool[0]
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโ raster โค
extent: Extent(X = (-124.409591, -114.409591), Y = (32.534156, 41.534156))
missingval: missing
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโ
โ โ 32.5342 33.0342 33.5342 โฆ 40.0342 40.5342 41.0342
-124.41 missing missing missing 1 1 missing
-123.91 missing missing missing 1 1 1
โฎ โฑ
-115.91 1 1 1 missing missing missing
-115.41 1 1 1 missing missing missing
-114.91 1 missing 1 โฆ 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 6.5.
= DimPoints(california_raster1) dp
โญโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโฎ
โ 20ร18 DimPoints{Tuple{Float64, Float64},2} โ
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโดโโโโโโโโโโโโโโโโโโโโโโโโโโโ dims โ
โ X Sampled{Float64} -124.409591:0.5:-114.909591 ForwardOrdered Regular Intervals{Start},
โ Y Sampled{Float64} 32.534156:0.5:41.034156 ForwardOrdered Regular Intervals{Start}
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโ
โ โ 32.5342 33.0342 โฆ 41.0342
-124.41 (-124.41, 32.5342) (-124.41, 33.0342) (-124.41, 41.0342)
-123.91 (-123.91, 32.5342) (-123.91, 33.0342) (-123.91, 41.0342)
-123.41 (-123.41, 32.5342) (-123.41, 33.0342) (-123.41, 41.0342)
-122.91 (-122.91, 32.5342) (-122.91, 33.0342) (-122.91, 41.0342)
โฎ โฑ
-116.91 (-116.91, 32.5342) (-116.91, 33.0342) (-116.91, 41.0342)
-116.41 (-116.41, 32.5342) (-116.41, 33.0342) (-116.41, 41.0342)
-115.91 (-115.91, 32.5342) (-115.91, 33.0342) (-115.91, 41.0342)
-115.41 (-115.41, 32.5342) (-115.41, 33.0342) (-115.41, 41.0342)
-114.91 (-114.91, 32.5342) (-114.91, 33.0342) โฆ (-114.91, 41.0342)
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.
in dims(california_raster1, X), y in dims(california_raster1, Y)] [(x, y) for x
20ร18 Matrix{Tuple{Float64, Float64}}:
(-124.41, 32.5342) (-124.41, 33.0342) โฆ (-124.41, 41.0342)
(-123.91, 32.5342) (-123.91, 33.0342) (-123.91, 41.0342)
(-123.41, 32.5342) (-123.41, 33.0342) (-123.41, 41.0342)
(-122.91, 32.5342) (-122.91, 33.0342) (-122.91, 41.0342)
(-122.41, 32.5342) (-122.41, 33.0342) (-122.41, 41.0342)
(-121.91, 32.5342) (-121.91, 33.0342) โฆ (-121.91, 41.0342)
(-121.41, 32.5342) (-121.41, 33.0342) (-121.41, 41.0342)
(-120.91, 32.5342) (-120.91, 33.0342) (-120.91, 41.0342)
(-120.41, 32.5342) (-120.41, 33.0342) (-120.41, 41.0342)
(-119.91, 32.5342) (-119.91, 33.0342) (-119.91, 41.0342)
(-119.41, 32.5342) (-119.41, 33.0342) โฆ (-119.41, 41.0342)
(-118.91, 32.5342) (-118.91, 33.0342) (-118.91, 41.0342)
(-118.41, 32.5342) (-118.41, 33.0342) (-118.41, 41.0342)
(-117.91, 32.5342) (-117.91, 33.0342) (-117.91, 41.0342)
(-117.41, 32.5342) (-117.41, 33.0342) (-117.41, 41.0342)
(-116.91, 32.5342) (-116.91, 33.0342) โฆ (-116.91, 41.0342)
(-116.41, 32.5342) (-116.41, 33.0342) (-116.41, 41.0342)
(-115.91, 32.5342) (-115.91, 33.0342) (-115.91, 41.0342)
(-115.41, 32.5342) (-115.41, 33.0342) (-115.41, 41.0342)
(-114.91, 32.5342) (-114.91, 33.0342) (-114.91, 41.0342)
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 6.6 shows the input vector layer, the rasterization results, and the points pnt
.
# Line rasterization
= plot(california_raster1; colormap = cgrad(:Set3; categorical = true))
fig, ax, plt lines!(ax, california_borders; color = :darkgrey, linewidth = 1)
scatter!(ax, vec(dp); markersize = 3, color = :black)
display(fig)
# Polygon rasterization
= plot(california_raster2; colormap = cgrad(:Set3; categorical = true))
fig, ax, plt lines!(ax, california_borders; color = :darkgrey, linewidth = 1)
scatter!(ax, vec(dp); markersize = 3, color = :black)
fig
boundary=:touches
boundary=:center
6.5 Spatial vectorization
Spatial vectorization is the counterpart of rasterization (Section 6.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 6.5.1)โconverting raster cells to rectangular polygons, representing pixel areas
- Raster to points (Section 6.5.2)โconverting raster cells to points, representing pixel centroids
- Raster to contours (Section 6.5.3)
Let us demonstrate all three in the given order.
6.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.
= GO.polygonize(src_grain) fc
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))))
To load this as a DataFrame
, we can take a shortcut and read it as GeoJSON (which is convertable to a DataFrame).
Support for direct conversion from GI.FeatureCollection
to a DataFrame
is in progress. TODO: add a link and track progress.
= DataFrame([GI.properties(f) for f in GI.getfeature(fc)])
df = [GI.geometry(f) for f in GI.getfeature(fc)]
df.geometry df
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 6.7.
= poly(df.geometry; color = df.value, strokecolor = :black, strokewidth = 0.75)
f, a, p Colorbar(f[1, 2], p)
f
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
.
6.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.
= DimTable(Raster("output/elev.tif")) dt
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.
= DimTable(Raster("output/elev.tif"), mergedims = (X, Y)) dt
DimensionalData.DimTable with 36 rows, 2 columns, and schema:
:XY Tuple{Float64, Float64}
:layer1 UInt8
This has created a DimTable
with a single dimension, :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.
= DataFrame(dt)
df metadata!(df, "GEOINTERFACE:geometrycolumns", (:XY,); style = :note)
DataFrames. df
Row | XY | layer1 |
---|---|---|
Tupleโฆ | UInt8 | |
1 | (-1.5, 1.0) | 1 |
2 | (-1.0, 1.0) | 2 |
3 | (-0.5, 1.0) | 3 |
โฎ | โฎ | โฎ |
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:
write("output/elev.gpkg", df)
GeoDataFrames.read("output/elev.gpkg") GeoDataFrames.
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 6.8 shows the input raster and the resulting point layer.
# Input raster
= plot(src_elev)
fig, ax, plt scatter!(ax, df.XY; color = df.layer1)
display(fig)
# Points
= plot(src_elev; alpha = 0.1)
fig, ax, plt scatter!(ax, df.XY; color = df.layer1, strokecolor = :black, strokewidth = 1)
fig
elev.tif
TODO: nodata pixels
6.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.
= contour(src_dem; levels = LinRange(0, 1200, 50), color = :black) f, ax, plt
TODO: gdal_contour (via ArchGDAL??)
It would be good to show how to use the provided GDAL executables thoughโฆ
6.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 6.5.2), and rasterizing points (Section 6.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
= GI.getexterior.(GI.getgeom(LibGEOS.unaryUnion(GI.GeometryCollection(nz.geom)))) .|> x -> GI.LineString(collect(GI.getpoint(x)))
coastline_linestrings = GI.MultiLineString(coastline_linestrings)
coastline = GO.reproject(coastline; target_crs = GI.crs(src_nz_elev), source_crs = GI.crs(nz)) coastline
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.
= 2/10
factor = Rasters.resample(src_nz_elev; size = round.(Int, size(src_nz_elev) .* factor), method = :average) r
โญโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโฎ
โ 223ร290 Raster{Union{Missing, 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 1 entry:
"filepath" => "/vsimem/tmp"
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโ 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.75196e6 4.74696e6 4.74196e6
9.95457e5 missing missing missing missing missing
1.00046e6 missing missing missing missing missing
โฎ โฑ โฎ
2.09546e6 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 6.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 6.10).
= plot(r)
fig, ax, plt lines!(ax, coastline; color = :red)
fig
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 6.5.2, but simply select the pixels that are not missing
.
= DimPoints(r)
dp = dp[r .=== missingval(r)] nonmissing_points
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.
= GO.distance.((coastline,), nonmissing_points) distances
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 6.4.1) the distances into our raster template.
= Rasters.rasterize(
img
last,
nonmissing_points;= r,
to = distances,
fill )
โญโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโฎ
โ 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 6.10.
= plot(img)
fig, ax, plt lines!(ax, coastline; color = :red)
Colorbar(fig[1, 2], plt; label = "Distance to coastline (m)")
fig