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 # plotting6 Raster-vector interactions
Prerequisites
This chapter requires importing the following packages:
Makie.set_theme!(
Heatmap = (; rasterize = 2),
Surface = (; rasterize = 2),
)It also relies on the following data files:
src_srtm = Raster("data/srtm.tif")
src_nlcd = Raster("data/nlcd.tif")
src_grain = Raster("output/grain.tif")
src_elev = Raster("output/elev.tif")
src_dem = Raster("data/dem.tif")
zion = GeoDataFrames.read("data/zion.gpkg")
zion_points = GeoDataFrames.read("data/zion_points.gpkg")
cycle_hire_osm = GeoDataFrames.read("data/cycle_hire_osm.gpkg")
us_states = GeoDataFrames.read("data/us_states.gpkg")
nz = GeoDataFrames.read("data/nz.gpkg")
src_nz_elev = Raster("data/nz_elev.tif")โญโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโฎ
โ 1115ร1450 Raster{Union{Missing, Float32},2} โ
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโดโโโโโโโโโโโโโโโโโโโโโโโโโโ dims โ
โ X Projected{Float64} 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 missing6.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.tifraster 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.gpkgvector layer representing the Zion National Park boundaries (aGeoDataFramenamedzion)
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.
zion = GO.reproject(zion; target_crs = GI.crs(src_srtm))1ร12 DataFrame
| Row | geom | UNIT_CODE | GIS_Notes | UNIT_NAME | DATE_EDIT | STATE | REGION | GNIS_ID | UNIT_TYPE | CREATED_BY | METADATA | PARKNAME |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Polygonโฆ | String | String | String | DateTime | String | String | String | String | String | String | String | |
| 1 | Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{CRS}}}, Nothing, WellKnownText{CRS}}(LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{CRS}}[LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{CRS}}([(-113.084, 37.1713), (-113.086, 37.1716), (-113.091, 37.1722), (-113.091, 37.1744), (-113.091, 37.178), (-113.091, 37.1817), (-113.091, 37.1853), (-113.091, 37.1889), (-113.091, 37.1926), (-113.091, 37.1962) โฆ (-113.06, 37.1713), (-113.062, 37.1713), (-113.064, 37.1713), (-113.069, 37.1713), (-113.071, 37.1713), (-113.073, 37.1713), (-113.078, 37.1713), (-113.08, 37.1713), (-113.082, 37.1713), (-113.084, 37.1713)], nothing, WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], nothing, WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]")) | ZION | Lands - http://landsnet.nps.gov/tractsnet/documents/ZION/Metadata/zion_metadata.xml | Zion National Park | 2017-06-22T00:00:00 | UT | IM | 1455157 | National Park | Lands | https://irma.nps.gov/App/Reference/Profile/2181118#Zion National Monument | Zion |
To mask the image, i.e., convert all pixels which do not intersect with the zion polygon to missing, we use the Rasters.mask function. mask supports any geometry, vector of geometries, feature collection, or table with a geometry column!
out_image_mask = Rasters.mask(src_srtm; with = zion)โญโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโฎ
โ 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
Note that since Julia has a native missing/NODATA value type, we donโt need to specify a NODATA value for the mask function.
However, it can sometimes be useful and more efficient to specify a sentinel value which Rasters treats as missing.
You can do this by specifying the missingval keyword argument, like so:
out_image_mask = Rasters.mask(src_srtm; with = zion, missingval = 9999)We can write this masked raster to file with Rasters.write, as usual:
Rasters.write("output/srtm_masked.tif", out_image_mask; force = true
)[ Info: `missingval` set to 65535 on disk"output/srtm_masked.tif"In Rasters.jl, cropping and masking are distinct operations. Cropping, which reduces the raster extent to the extent of the vector layer, is accomplished with the crop function.
Here, we simply pass the zion feature table to the to keyword argument, which indicates what to crop the raster โtoโ. We also set the touches keyword argument to true, to specify that pixels that partially overlap with the vector layer are included in the output.
out_image_crop = Rasters.crop(src_srtm; to = zion, touches = true)โญโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโฎ
โ 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 0x06d4You can also assemble an extent manually, using Extents.Extent, or extract one using GI.extent.
We can crop our masked raster as well:
out_image_mask_crop = Rasters.crop(out_image_mask; to = zion, touches = true)โญโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโฎ
โ 439ร436 Raster{Union{Missing, UInt16},2} โ
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโดโโโโโโโโโโโโโโโโโโโโโโโโโโโโโ dims โ
โ X Projected{Float64} 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 missingand we write it to file as usual:
Rasters.write("output/srtm_masked_cropped.tif", out_image_mask_crop; force = true)[ Info: `missingval` set to 65535 on disk"output/srtm_masked_cropped.tif"Figure 6.1 shows the original raster, and the three masking and/or cropping results.
fig = Figure(size = (600, 600))
ax1 = Axis(fig[1, 1]; title = "Original")
plot!(ax1, src_srtm)
poly!(ax1, zion.geom; color = :transparent, strokecolor = :black, strokewidth = 0.75)
ax2 = Axis(fig[1, 2]; title = "Masked")
plot!(ax2, out_image_mask)
poly!(ax2, zion.geom; color = :transparent, strokecolor = :black, strokewidth = 0.75)
ax3 = Axis(fig[2, 1]; title = "Cropped")
plot!(ax3, out_image_crop)
poly!(ax3, zion.geom; color = :transparent, strokecolor = :black, strokewidth = 0.75)
ax4 = Axis(fig[2, 2]; title = "Masked+Cropped")
plot!(ax4, out_image_mask_crop)
poly!(ax4, zion.geom; color = :transparent, strokecolor = :black, strokewidth = 0.75)
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.extractfunction - To polygons (Section 6.3.3), via the
Rasters.zonalfunction
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).
fig, ax, plt = plot(src_srtm)
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.
result1 = Rasters.extract(src_srtm, zion_points; geometry = false) .|> first30-element Vector{UInt16}:
0x070a
0x0981
0x075e
0x055a
0x05ac
0x0663
0x0564
0x07f0
0x0726
0x0744
โฎ
0x0711
0x0722
0x060e
0x0707
0x0836
0x0846
0x055c
0x0771
0x0626The first argument is the raster from which to extract values, and the second is the vector object (or collection of objects) according to which to extract the values.
Note
Rasters.jl does not yet support interpolation, so the values extracted are the values of the nearest cell center.
This corresponds to interpolate='nearest' in the Python rasterstats package.
Either way, the resulting object is a vector of raster values, corresponding to zion_points. For example, here are the elevations of the first five points.
result1[1:5]5-element Vector{UInt16}:
0x070a
0x0981
0x075e
0x055a
0x05acTo 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.
zion_points[!, "elev1"] = result1
zion_points30ร2 DataFrame
24 rows omitted
| 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.
coords = [[-113.2, 37.45], [-112.9, 37.2]]
zion_transect = GI.LineString(coords)GeoInterface.Wrappers.LineString{false, false, Vector{Vector{Float64}}, Nothing, Nothing}([[-113.2, 37.45], [-112.9, 37.2]], nothing, nothing)The utility of extracting heights from a linear selector is illustrated by imagining that you are planning a hike. The method demonstrated below provides an โelevation profileโ of the route (the line does not need to be straight), useful for estimating how long it will take 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.
zion_transect_utm = GO.reproject(zion_transect; target_crs = GFT.EPSG(32612), source_crs = GFT.EPSG(4326))GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, EPSG{1}}([(305399.67208180577, 4.147066650206682e6), (331380.8917453843, 4.1187500947884847e6)], nothing, EPSG{1}((32612,)))The printout of the new geometry shows this is still a straight line between two points, only with coordinates in a projected CRS.
Iโve chosen to differ from the Python treatment here - instead of selecting some number of points along the line explicitly, I will segmentize the line and extract the points. This is less precise, but we donโt have the API to do arclength interpolation in GeometryOps yet. Hopefully this will be added soon.
cf. https://github.com/JuliaGeo/GeometryOps.jl/issues/210
Here, we interpolate points along the line using GO.segmentize.
We first compute the length of the line, and then use this to segmentize the line into approximately 250 points.
_centroid, linelen = GO.centroid_and_length(zion_transect_utm)
zion_transect_line = GO.segmentize(zion_transect_utm; max_distance = linelen / 250)GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, EPSG{1}}([(305399.67208180577, 4.147066650206682e6), (305503.59696046007, 4.146953383985009e6), (305607.5218391144, 4.1468401177633363e6), (305711.44671776873, 4.1467268515416635e6), (305815.37159642304, 4.1466135853199908e6), (305919.29647507734, 4.146500319098318e6), (306023.22135373164, 4.1463870528766452e6), (306127.14623238595, 4.1462737866549725e6), (306231.0711110403, 4.1461605204332997e6), (306334.9959896946, 4.1460472542116265e6) โฆ (330445.5678374955, 4.11976949078354e6), (330549.4927161498, 4.119656224561867e6), (330653.41759480414, 4.119542958340194e6), (330757.34247345844, 4.1194296921185213e6), (330861.26735211274, 4.1193164258968486e6), (330965.19223076705, 4.119203159675176e6), (331069.11710942135, 4.119089893453503e6), (331173.0419880757, 4.1189766272318303e6), (331276.96686673, 4.1188633610101575e6), (331380.8917453843, 4.1187500947884847e6)], nothing, EPSG{1}((32612,)))This gives us a collection of 251 points along the line. We can extract the points by using GI.getpoint on the line, and then reproject the points to the CRS of the raster.
zion_transect_pnt = GO.reproject(GI.getpoint(zion_transect_line); target_crs = GI.crs(src_srtm), source_crs = GI.crs(zion_transect_line))251-element Vector{Tuple{Float64, Float64}}:
(-113.2, 37.45)
(-113.1987960226292, 37.449001661297636)
(-113.19759207733796, 37.448003309143594)
(-113.1963881641247, 37.44700494353919)
(-113.1951842829878, 37.44600656448574)
(-113.19398043392566, 37.44500817198453)
(-113.19277661693668, 37.444009766036885)
(-113.19157283201925, 37.4430113466441)
(-113.19036907917175, 37.4420129138075)
(-113.18916535839257, 37.44101446752837)
โฎ
(-112.90956920045984, 37.20801281544639)
(-112.90837293949568, 37.20701125948384)
(-112.90717671022236, 37.20600969038492)
(-112.90598051263827, 37.205008108150906)
(-112.90478434674178, 37.2040065127831)
(-112.90358821253133, 37.20300490428282)
(-112.9023921100053, 37.20200328265134)
(-112.90119603916204, 37.201001647889974)
(-112.90000000000002, 37.20000000000002)Finally, we extract the elevation values for each point in our transect and combine the information with zion_transect_pnt (after โpromotingโ it to a 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.
zion_transect_pnt = DataFrame(geometry = zion_transect_pnt)
result = Rasters.extract(src_srtm, zion_transect_pnt; geometry = false) .|> first
zion_transect_pnt[!, "elev"] = result251-element Vector{UInt16}:
0x07d1
0x07f1
0x07e4
0x07c5
0x076e
0x074d
0x06fe
0x06e7
0x06d3
0x06ba
โฎ
0x0720
0x079f
0x0769
0x0737
0x072d
0x0737
0x0729
0x070a
0x06e3We also want to visualize an elevation profile along this line, so we compute the distances along the line manually. Using GO.distance, we can get the distance between successive points along the line. Then, we add 0.0 to the beginning of the array of distances, and sum along that array to get the cumulative distance along the line.
# Compute distances between successive points along the line
distances_between_points = GO.distance.(zion_transect_pnt.geometry[1:end-1], zion_transect_pnt.geometry[2:end])
# Compute cumulative distances along the line
cumulative_distances = cumsum([0.0; distances_between_points])
zion_transect_pnt[!, "dist"] = cumulative_distances
zion_transect_pnt251ร3 DataFrame
245 rows omitted
| 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
fig, ax, plt = plot(src_srtm)
lines!(ax, zion_transect; color = :black)
poly!(ax, zion.geom; color = :transparent, strokecolor = :white, strokewidth = 0.75)
display(fig)
# Elevation profile
fig, ax, plt = lines(
zion_transect_pnt.dist,
zion_transect_pnt.elev;
axis = (;
xlabel = "Distance (m)",
ylabel = "Elevation (m)",
)
)
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`
rmean = Rasters.zonal(mean, src_srtm; of = zion)
rmin = Rasters.zonal(minimum, src_srtm; of = zion)
rmax = Rasters.zonal(maximum, src_srtm; of = zion)
result = (rmean, rmin, rmax)(Union{Missing, Float64}[1818.211830154405], Union{Missing, UInt16}[0x0462], Union{Missing, UInt16}[0x0a65])Itโs straightforward to transform the result to a DataFrame as well:
DataFrame(mean = rmean, min = rmin, max = rmax)1ร3 DataFrame
| Row | mean | min | max |
|---|---|---|---|
| Float64? | UInt16? | UInt16? | |
| 1 | 1818.21 | 1122 | 2661 |
Because there is only one polygon in the example, single-element vectors are returned. However, if zion was composed of more than one polygon, we would accordingly get more elements in the returned vectors. The result provides useful summaries, for example that the maximum height in the park is 2661 \(m\) above sea level.
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.
out_image = Rasters.mask(src_nlcd; with = GO.reproject(zion; target_crs = GI.crs(src_nlcd)))
using StatsBase
counts = StatsBase.countmap(out_image)Dict{Union{Missing, UInt8}, Int64} with 8 entries:
0x05 => 203701
0x04 => 298299
0x06 => 235
0x07 => 62
0x02 => 4205
missing => 852741
0x08 => 679
0x03 => 98285According 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
fig, ax, plt = plot(src_srtm)
poly!(ax, zion.geom; color = :transparent, strokecolor = :black, strokewidth = 0.75)
display(fig)
# Categorical raster
fig, ax, plt = plot(src_nlcd; colormap = cgrad(:Set3; categorical = true), source = GI.crs(src_nlcd), axis = (; type = GeoAxis, dest = GI.crs(src_nlcd)))
poly!(ax, zion.geom; source = GI.crs(zion.geom[1]), color = :transparent, strokecolor = :black, strokewidth = 0.75)
cm = Colorbar(fig[1, 2], plt)
ax.xgridvisible = false
ax.ygridvisible = false
fig # TODO: make GeoMakie better on small extents
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.
cycle_hire_osm_projected = GO.reproject(cycle_hire_osm; target_crs = GFT.EPSG(27700))540ร6 DataFrame
534 rows omitted
| Row | geom | osm_id | name | capacity | cyclestreets_id | description |
|---|---|---|---|---|---|---|
| Tupleโฆ | String | String? | Float64? | String? | String? | |
| 1 | (5.32354e5, 1.82858e5) | 108539 | Windsor Terrace | 14.0 | missing | missing |
| 2 | (5.29848e5, 1.83337e5) | 598093293 | Pancras Road, King's Cross | missing | missing | missing |
| 3 | (5.30636e5, 182609.0) | 772536185 | Clerkenwell, Ampton Street | 11.0 | missing | missing |
| โฎ | โฎ | โฎ | โฎ | โฎ | โฎ | โฎ |
| 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.
ch_raster1 = Rasters.rasterize(
last, # reducer
cycle_hire_osm_projected; # data
fill = 1,
size = (1000, 1000) # specify size in "pixels"
)โญโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโฎ
โ 1000ร1000 Raster{Union{Missing, Int64},2} last โ
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโดโโโโโโโโโโโโโโโโโโโโโโโ dims โ
โ X Sampled{Float64} 523038.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 missingIn our second variant of point rasterization, we count the number of bike hire stations. To do that, we use the fixed value of 1 (same as in the last example), but this time combined with the reducer=sum argument. That way, multiple values burned into the same pixel are summed, rather than replaced keeping last (which is the default). The new output, ch_raster2, shows the number of cycle hire points in each grid cell.
ch_raster2 = Rasters.rasterize(
sum, # reducer
cycle_hire_osm_projected; # data
fill = 1,
size = (1000, 1000) # specify size in "pixels"
)โญโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโฎ
โ 1000ร1000 Raster{Union{Missing, Int64},2} sum โ
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโดโโโโโโโโโโโโโโโโโโโโโโโโ dims โ
โ X Sampled{Float64} 523038.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 missingThe cycle hire locations have different numbers of bicycles described by the capacity variable, raising the question, what is the capacity in each grid cell? To calculate that, in our third point rasterization variant we sum the field ('capacity') rather than the fixed values of 1.
This is extremely simple to run, but we will show how to do this two ways: first, by passing the column name in the feature collection to fill.
ch_raster3 = Rasters.rasterize(
sum, # reducer
cycle_hire_osm_projected; # data
fill = :capacity,
size = (1000, 1000) # specify size in "pixels"
)โญโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโฎ
โ 1000ร1000 Raster{Union{Missing, Float64},2} capacity โ
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโดโโโโโโโโโโโโโโโโโ dims โ
โ X Sampled{Float64} 523038.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 missingSecond, by passing the vectors of geometries and values separately.
ch_raster3 = Rasters.rasterize(
sum, # reducer
cycle_hire_osm_projected.geom; # data
fill = cycle_hire_osm_projected.capacity,
crs = GI.crs(cycle_hire_osm_projected),
size = (1000, 1000) # specify size in "pixels"
)โญโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโฎ
โ 1000ร1000 Raster{Union{Missing, Float64},2} sum โ
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโดโโโโโโโโโโโโโโโโโโโโโโ dims โ
โ X Projected{Float64} 523038.61452275474:15.684723604717:538707.653403867 ForwardOrdered Regular Intervals{Start},
โ Y Projected{Float64} 174934.65727249475:10.03675127048517:184961.37179170942 ForwardOrdered Regular Intervals{Start}
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโ metadata โค
Metadata of Dict{Any, Any}()
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโ raster โค
extent: Extent(X = (523038.61452275474, 538723.3381274717), Y = (174934.65727249475, 184971.40854297992))
missingval: missing
crs: EPSG:27700
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโ
โ โ 1.74935e5 1.74945e5 โฆ 1.84941e5 1.84951e5 1.84961e5
5.23039e5 missing missing missing missing missing
5.23054e5 missing missing missing missing missing
523070.0 missing missing missing missing missing
โฎ โฑ โฎ
5.38676e5 missing missing missing missing missing
538692.0 missing missing missing missing missing
5.38708e5 missing missing โฆ missing missing missingThe 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
nonmissing_df = dropmissing(cycle_hire_osm_projected, [:capacity, :geom])
f, a, p = scatter(nonmissing_df.geom; color = nonmissing_df.capacity)
Colorbar(f[1, 2], p)
display(f)
# Presence/Absence
plot(ch_raster1) |> display
# Point counts
plot(ch_raster2) |> display
# Summed attribute values
plot(ch_raster3)
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.
california = us_states[ us_states[!, "NAME"] .== "California", :]1ร7 DataFrame
| Row | geom | GEOID | NAME | REGION | AREA | total_pop_10 | total_pop_15 |
|---|---|---|---|---|---|---|---|
| IGeometrโฆ | String | String | String | Float64 | Float64 | Float64 | |
| 1 | Geometry: wkbMultiPolygon | 06 | California | West | 4.09747e5 | 3.66373e7 | 3.84215e7 |
Second, we obtain the borders of the polygon as a `โMultiLineStringโ
california_geom = only(california.geom)
california_borders = GI.MultiLineString(GI.LineString.(GI.getexterior.(GI.getgeom(california_geom))); crs = GI.crs(california_geom)) # TODO: make this a lot better....GeoInterface.Wrappers.MultiLineString{false, false, Vector{GeoInterface.Wrappers.LineString{false, false, ArchGDAL.IGeometry{ArchGDAL.wkbLineString}, Nothing, Nothing}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LineString{false, false, ArchGDAL.IGeometry{ArchGDAL.wkbLineString}, Nothing, Nothing}[GeoInterface.Wrappers.LineString{false, false, ArchGDAL.IGeometry{ArchGDAL.wkbLineString}, Nothing, Nothing}(Geometry: LINEARRING (-118.603375 33.478098,-118.368301 33.4 ... 8098), nothing, nothing), GeoInterface.Wrappers.LineString{false, false, ArchGDAL.IGeometry{ArchGDAL.wkbLineString}, Nothing, Nothing}(Geometry: LINEARRING (-118.605534 33.030999,-118.369984 32.8 ... 0999), nothing, nothing), GeoInterface.Wrappers.LineString{false, false, ArchGDAL.IGeometry{ArchGDAL.wkbLineString}, Nothing, Nothing}(Geometry: LINEARRING (-119.919155 34.07728,-119.755521 34.05 ... 7728), nothing, nothing), GeoInterface.Wrappers.LineString{false, false, ArchGDAL.IGeometry{ArchGDAL.wkbLineString}, Nothing, Nothing}(Geometry: LINEARRING (-120.247393 34.001911,-120.073609 34.0 ... 1911), nothing, nothing), GeoInterface.Wrappers.LineString{false, false, ArchGDAL.IGeometry{ArchGDAL.wkbLineString}, Nothing, Nothing}(Geometry: LINEARRING (-124.211605 41.99846,-123.347562 41.99 ... 9846), nothing, nothing)], nothing, WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"NAD83\",DATUM[\"North_American_Datum_1983\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"6269\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4269\"]]"))Finally, we rasterize california_borders on a grid with resolution of 0.5 degrees per pixel.
california_raster1 = Rasters.rasterize(
last,
california_borders;
fill = 1,
res = 0.5, # degrees - this is in units of GI.crs(california_borders)
boundary = :touches,
)โญโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโฎ
โ 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 missingCompare 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).
california_raster2 = Rasters.rasterize(
last,
california;
geometrycolumn = :geom,
fill = 1,
res = 0.5,
boundary = :center,
)โญโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโฎ
โ 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 missingTo 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.
dp = DimPoints(california_raster1)โญโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโฎ
โ 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.
[(x, y) for x in dims(california_raster1, X), y in dims(california_raster1, Y)]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
fig, ax, plt = plot(california_raster1; colormap = cgrad(:Set3; categorical = true))
lines!(ax, california_borders; color = :darkgrey, linewidth = 1)
scatter!(ax, vec(dp); markersize = 3, color = :black)
display(fig)
# Polygon rasterization
fig, ax, plt = plot(california_raster2; colormap = cgrad(:Set3; categorical = true))
lines!(ax, california_borders; color = :darkgrey, linewidth = 1)
scatter!(ax, vec(dp); markersize = 3, color = :black)
fig
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.
fc = GO.polygonize(src_grain)GeoInterface.Wrappers.FeatureCollection{Vector{GeoInterface.Wrappers.Feature{@NamedTuple{geometry::GeoInterface.Wrappers.MultiPolygon{false, false, Vector{GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}, value::UInt8}, WellKnownText{GeoFormatTypes.CRS}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}}}, WellKnownText{GeoFormatTypes.CRS}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}}(GeoInterface.Wrappers.Feature{@NamedTuple{geometry::GeoInterface.Wrappers.MultiPolygon{false, false, Vector{GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}, value::UInt8}, WellKnownText{GeoFormatTypes.CRS}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}}[GeoInterface.Wrappers.Feature{@NamedTuple{geometry::GeoInterface.Wrappers.MultiPolygon{false, false, Vector{GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}, value::UInt8}, WellKnownText{GeoFormatTypes.CRS}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}}((geometry = GeoInterface.Wrappers.MultiPolygon{false, false, Vector{GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}([(1.5, 1.5), (1.5, 3.5), (2.5, 3.5), (2.5, 4.5), (0.5, 4.5), (0.5, 1.5), (1.5, 1.5)], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 2.5), Y = (1.5, 4.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 2.5), Y = (1.5, 4.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]")), GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}([(1.5, 1.5), (1.5, 0.5), (2.5, 0.5), (2.5, 1.5), (1.5, 1.5)], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (1.5, 2.5), Y = (0.5, 1.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (1.5, 2.5), Y = (0.5, 1.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]")), GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}([(5.5, 5.5), (5.5, 6.5), (4.5, 6.5), (4.5, 5.5), (5.5, 5.5)], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (4.5, 5.5), Y = (5.5, 6.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (4.5, 5.5), Y = (5.5, 6.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]")), GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}([(4.5, 3.5), (3.5, 3.5), (3.5, 2.5), (2.5, 2.5), (2.5, 1.5), (4.5, 1.5), (4.5, 2.5), (5.5, 2.5), (5.5, 3.5), (4.5, 3.5)], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (2.5, 5.5), Y = (1.5, 3.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (2.5, 5.5), Y = (1.5, 3.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 5.5), Y = (0.5, 6.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]")), value = 0x00), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"), Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 5.5), Y = (0.5, 6.5)))), GeoInterface.Wrappers.Feature{@NamedTuple{geometry::GeoInterface.Wrappers.MultiPolygon{false, false, Vector{GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}, value::UInt8}, WellKnownText{GeoFormatTypes.CRS}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}}((geometry = GeoInterface.Wrappers.MultiPolygon{false, false, Vector{GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}([(0.5, 4.5), (2.5, 4.5), (2.5, 3.5), (6.5, 3.5), (6.5, 5.5), (4.5, 5.5), (4.5, 4.5), (3.5, 4.5), (3.5, 5.5), (2.5, 5.5), (2.5, 6.5), (1.5, 6.5), (1.5, 5.5), (0.5, 5.5), (0.5, 4.5)], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 6.5), Y = (3.5, 6.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 6.5), Y = (3.5, 6.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]")), GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}([(1.5, 1.5), (0.5, 1.5), (0.5, 0.5), (1.5, 0.5), (1.5, 1.5)], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 1.5), Y = (0.5, 1.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 1.5), Y = (0.5, 1.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]")), GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}([(5.5, 2.5), (5.5, 1.5), (6.5, 1.5), (6.5, 2.5), (5.5, 2.5)], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (5.5, 6.5), Y = (1.5, 2.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (5.5, 6.5), Y = (1.5, 2.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]")), GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}([(2.5, 0.5), (3.5, 0.5), (3.5, 1.5), (2.5, 1.5), (2.5, 0.5)], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (2.5, 3.5), Y = (0.5, 1.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (2.5, 3.5), Y = (0.5, 1.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 6.5), Y = (0.5, 6.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]")), value = 0x01), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"), Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 6.5), Y = (0.5, 6.5)))), GeoInterface.Wrappers.Feature{@NamedTuple{geometry::GeoInterface.Wrappers.MultiPolygon{false, false, Vector{GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}, value::UInt8}, WellKnownText{GeoFormatTypes.CRS}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}}((geometry = GeoInterface.Wrappers.MultiPolygon{false, false, Vector{GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}([(1.5, 1.5), (2.5, 1.5), (2.5, 2.5), (3.5, 2.5), (3.5, 3.5), (1.5, 3.5), (1.5, 1.5)], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (1.5, 3.5), Y = (1.5, 3.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (1.5, 3.5), Y = (1.5, 3.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]")), GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}([(5.5, 5.5), (6.5, 5.5), (6.5, 6.5), (5.5, 6.5), (5.5, 5.5)], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (5.5, 6.5), Y = (5.5, 6.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (5.5, 6.5), Y = (5.5, 6.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]")), GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}([(4.5, 6.5), (2.5, 6.5), (2.5, 5.5), (3.5, 5.5), (3.5, 4.5), (4.5, 4.5), (4.5, 6.5)], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (2.5, 4.5), Y = (4.5, 6.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (2.5, 4.5), Y = (4.5, 6.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]")), GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}([(5.5, 0.5), (6.5, 0.5), (6.5, 1.5), (5.5, 1.5), (5.5, 2.5), (4.5, 2.5), (4.5, 1.5), (3.5, 1.5), (3.5, 0.5), (5.5, 0.5)], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (3.5, 6.5), Y = (0.5, 2.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (3.5, 6.5), Y = (0.5, 2.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]")), GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}([(6.5, 3.5), (5.5, 3.5), (5.5, 2.5), (6.5, 2.5), (6.5, 3.5)], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (5.5, 6.5), Y = (2.5, 3.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (5.5, 6.5), Y = (2.5, 3.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]")), GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{GeoFormatTypes.CRS}}([(1.5, 6.5), (0.5, 6.5), (0.5, 5.5), (1.5, 5.5), (1.5, 6.5)], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 1.5), Y = (5.5, 6.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 1.5), Y = (5.5, 6.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"))], Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 6.5), Y = (0.5, 6.5))), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]")), value = 0x02), WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"), Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 6.5), Y = (0.5, 6.5))))], WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"), Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 6.5), Y = (0.5, 6.5))))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.
df = DataFrame([GI.properties(f) for f in GI.getfeature(fc)])
df.geometry = [GI.geometry(f) for f in GI.getfeature(fc)]
df3ร2 DataFrame
| Row | value | geometry |
|---|---|---|
| UInt8 | MultiPolโฆ | |
| 1 | 0 | MultiPolygon{false, false, Vector{Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}([(1.5, 1.5), (1.5, 3.5), (2.5, 3.5), (2.5, 4.5), (0.5, 4.5), (0.5, 1.5), (1.5, 1.5)], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 2.5), Y = (1.5, 4.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 2.5), Y = (1.5, 4.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]")), Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}([(1.5, 1.5), (1.5, 0.5), (2.5, 0.5), (2.5, 1.5), (1.5, 1.5)], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (1.5, 2.5), Y = (0.5, 1.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (1.5, 2.5), Y = (0.5, 1.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]")), Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}([(5.5, 5.5), (5.5, 6.5), (4.5, 6.5), (4.5, 5.5), (5.5, 5.5)], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (4.5, 5.5), Y = (5.5, 6.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (4.5, 5.5), Y = (5.5, 6.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]")), Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}([(4.5, 3.5), (3.5, 3.5), (3.5, 2.5), (2.5, 2.5), (2.5, 1.5), (4.5, 1.5), (4.5, 2.5), (5.5, 2.5), (5.5, 3.5), (4.5, 3.5)], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (2.5, 5.5), Y = (1.5, 3.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (2.5, 5.5), Y = (1.5, 3.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 5.5), Y = (0.5, 6.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]")) |
| 2 | 1 | MultiPolygon{false, false, Vector{Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}([(0.5, 4.5), (2.5, 4.5), (2.5, 3.5), (6.5, 3.5), (6.5, 5.5), (4.5, 5.5), (4.5, 4.5), (3.5, 4.5), (3.5, 5.5), (2.5, 5.5), (2.5, 6.5), (1.5, 6.5), (1.5, 5.5), (0.5, 5.5), (0.5, 4.5)], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 6.5), Y = (3.5, 6.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 6.5), Y = (3.5, 6.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]")), Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}([(1.5, 1.5), (0.5, 1.5), (0.5, 0.5), (1.5, 0.5), (1.5, 1.5)], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 1.5), Y = (0.5, 1.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 1.5), Y = (0.5, 1.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]")), Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}([(5.5, 2.5), (5.5, 1.5), (6.5, 1.5), (6.5, 2.5), (5.5, 2.5)], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (5.5, 6.5), Y = (1.5, 2.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (5.5, 6.5), Y = (1.5, 2.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]")), Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}([(2.5, 0.5), (3.5, 0.5), (3.5, 1.5), (2.5, 1.5), (2.5, 0.5)], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (2.5, 3.5), Y = (0.5, 1.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (2.5, 3.5), Y = (0.5, 1.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 6.5), Y = (0.5, 6.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]")) |
| 3 | 2 | MultiPolygon{false, false, Vector{Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}([(1.5, 1.5), (2.5, 1.5), (2.5, 2.5), (3.5, 2.5), (3.5, 3.5), (1.5, 3.5), (1.5, 1.5)], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (1.5, 3.5), Y = (1.5, 3.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (1.5, 3.5), Y = (1.5, 3.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]")), Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}([(5.5, 5.5), (6.5, 5.5), (6.5, 6.5), (5.5, 6.5), (5.5, 5.5)], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (5.5, 6.5), Y = (5.5, 6.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (5.5, 6.5), Y = (5.5, 6.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]")), Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}([(4.5, 6.5), (2.5, 6.5), (2.5, 5.5), (3.5, 5.5), (3.5, 4.5), (4.5, 4.5), (4.5, 6.5)], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (2.5, 4.5), Y = (4.5, 6.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (2.5, 4.5), Y = (4.5, 6.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]")), Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}([(5.5, 0.5), (6.5, 0.5), (6.5, 1.5), (5.5, 1.5), (5.5, 2.5), (4.5, 2.5), (4.5, 1.5), (3.5, 1.5), (3.5, 0.5), (5.5, 0.5)], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (3.5, 6.5), Y = (0.5, 2.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (3.5, 6.5), Y = (0.5, 2.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]")), Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}([(6.5, 3.5), (5.5, 3.5), (5.5, 2.5), (6.5, 2.5), (6.5, 3.5)], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (5.5, 6.5), Y = (2.5, 3.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (5.5, 6.5), Y = (2.5, 3.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]")), Polygon{false, false, Vector{LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}(LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}[LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}, WellKnownText{CRS}}([(1.5, 6.5), (0.5, 6.5), (0.5, 5.5), (1.5, 5.5), (1.5, 6.5)], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 1.5), Y = (5.5, 6.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 1.5), Y = (5.5, 6.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]"))], Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (0.5, 6.5), Y = (0.5, 6.5))), WellKnownText{CRS}(CRS(), "GEOGCS[\\"WGS 84\\",DATUM[\\"WGS_1984\\",SPHEROID[\\"WGS 84\\",6378137,298.257223563,AUTHORITY[\\"EPSG\\",\\"7030\\"]],AUTHORITY[\\"EPSG\\",\\"6326\\"]],PRIMEM[\\"Greenwich\\",0,AUTHORITY[\\"EPSG\\",\\"8901\\"]],UNIT[\\"degree\\",0.0174532925199433,AUTHORITY[\\"EPSG\\",\\"9122\\"]],AXIS[\\"Latitude\\",NORTH],AXIS[\\"Longitude\\",EAST],AUTHORITY[\\"EPSG\\",\\"4326\\"]]")) |
The polygon layer df is shown in Figure 6.7.
f, a, p = poly(df.geometry; color = df.value, strokecolor = :black, strokewidth = 0.75)
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.
dt = DimTable(Raster("output/elev.tif"))DimensionalData.DimTable with 36 rows, 3 columns, and schema:
:X Float64
:Y Float64
:layer1 UInt8Notice that this has three columns, :X, :Y, and :layer1, corresponding to the pixel centroids and elevation values. But what if we want to treat the X and Y dimensionas as point geometries?
DimTable has a mergedims keyword argument for this, which allows us to merge the X and Y dimensions into a single dimension.
dt = DimTable(Raster("output/elev.tif"), mergedims = (X, Y))DimensionalData.DimTable with 36 rows, 2 columns, and schema:
:XY Tuple{Float64, Float64}
:layer1 UInt8This 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.
df = DataFrame(dt)
DataFrames.metadata!(df, "GEOINTERFACE:geometrycolumns", (:XY,); style = :note)
df36ร2 DataFrame
30 rows omitted
| Row | XY | layer1 |
|---|---|---|
| Tupleโฆ | UInt8 | |
| 1 | (-1.5, 1.0) | 1 |
| 2 | (-1.0, 1.0) | 2 |
| 3 | (-0.5, 1.0) | 3 |
| โฎ | โฎ | โฎ |
| 34 | (-1.11022e-16, -1.5) | 34 |
| 35 | (0.5, -1.5) | 35 |
| 36 | (1.0, -1.5) | 36 |
scatter(df.XY; color = df.layer1)
We can even save this to a file trivially easily:
GeoDataFrames.write("output/elev.gpkg", df)
GeoDataFrames.read("output/elev.gpkg")36ร2 DataFrame
30 rows omitted
| Row | XY | layer1 |
|---|---|---|
| IGeometrโฆ | Int16 | |
| 1 | Geometry: wkbPoint | 1 |
| 2 | Geometry: wkbPoint | 2 |
| 3 | Geometry: wkbPoint | 3 |
| โฎ | โฎ | โฎ |
| 34 | Geometry: wkbPoint | 34 |
| 35 | Geometry: wkbPoint | 35 |
| 36 | Geometry: wkbPoint | 36 |
Figure 6.8 shows the input raster and the resulting point layer.
# Input raster
fig, ax, plt = plot(src_elev)
scatter!(ax, df.XY; color = df.layer1)
display(fig)
# Points
fig, ax, plt = plot(src_elev; alpha = 0.1)
scatter!(ax, df.XY; color = df.layer1, strokecolor = :black, strokewidth = 1)
fig
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.
f, ax, plt = contour(src_dem; levels = LinRange(0, 1200, 50), color = :black)
TODO: gdal_contour (via ArchGDAL??)
It would be good to show how to use the provided GDAL executables thoughโฆ
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
coastline_linestrings = GI.getexterior.(GI.getgeom(LibGEOS.unaryUnion(GI.GeometryCollection(nz.geom)))) .|> x -> GI.LineString(collect(GI.getpoint(x)))
coastline = GI.MultiLineString(coastline_linestrings)
coastline = GO.reproject(coastline; target_crs = GI.crs(src_nz_elev), source_crs = GI.crs(nz))GeoInterface.Wrappers.MultiLineString{false, false, Vector{GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}(GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}[GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}([(1.2299979013e6, 4.7934322369e6), (1.2216031693e6, 4.7943017685e6), (1.2177895399e6, 4.7904898277e6), (1.2311152067e6, 4.787137924e6), (1.2361839211e6, 4.7843314255e6), (1.2325099445e6, 4.7718700237e6), (1.2173644933e6, 4.7681066692e6), (1.2080683451e6, 4.7605677875e6), (1.199129799e6, 4.7634183939e6), (1.1938565725e6, 4.7574957427e6) โฆ (1.1964031733e6, 4.7848147499e6), (1.2011220591e6, 4.7884807601e6), (1.2003610396e6, 4.7969222422e6), (1.1955720567e6, 4.8084895651e6), (1.1989855975e6, 4.8153873633e6), (1.2077044806e6, 4.8171297889e6), (1.2161546733e6, 4.8132662203e6), (1.2184777781e6, 4.8064984488e6), (1.2297297347e6, 4.7987186067e6), (1.2299979013e6, 4.7934322369e6)], nothing, WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "PROJCS[\"unknown\",GEOGCS[\"unknown\",DATUM[\"Unknown based on GRS80 ellipsoid\",SPHEROID[\"GRS 1980\",6378137,298.257222101004,AUTHORITY[\"EPSG\",\"7019\"]],TOWGS84[0,0,0,0,0,0,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",173],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",1600000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]")), GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}([(1.1010682698e6, 4.9169506843e6), (1.0962776826e6, 4.9262213416e6), (1.1070468733e6, 4.9308841834e6), (1.1101943531e6, 4.9284267394e6), (1.1114045691e6, 4.9201779986e6), (1.1088006957e6, 4.914502166e6), (1.1010682698e6, 4.9169506843e6)], nothing, WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "PROJCS[\"unknown\",GEOGCS[\"unknown\",DATUM[\"Unknown based on GRS80 ellipsoid\",SPHEROID[\"GRS 1980\",6378137,298.257222101004,AUTHORITY[\"EPSG\",\"7019\"]],TOWGS84[0,0,0,0,0,0,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",173],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",1600000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]")), GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}([(1.2822947146e6, 4.8243903885e6), (1.2777592139e6, 4.8327835132e6), (1.2629460033e6, 4.8352306146e6), (1.2563895979e6, 4.832365038e6), (1.2436815529e6, 4.8358664271e6), (1.237141582e6, 4.8332854896e6), (1.2365580519e6, 4.8400088904e6), (1.2301828929e6, 4.8518834763e6), (1.2202312165e6, 4.8564065645e6), (1.2167650225e6, 4.8512224604e6) โฆ (1.3536635084e6, 4.861558587e6), (1.3540190062e6, 4.8508162544e6), (1.3384831555e6, 4.837432841e6), (1.3299957273e6, 4.8375532567e6), (1.3247154647e6, 4.8316347032e6), (1.3151459743e6, 4.8300668586e6), (1.310825215e6, 4.826096002e6), (1.3025458256e6, 4.8269904113e6), (1.2945187261e6, 4.8232725043e6), (1.2822947146e6, 4.8243903885e6)], nothing, WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "PROJCS[\"unknown\",GEOGCS[\"unknown\",DATUM[\"Unknown based on GRS80 ellipsoid\",SPHEROID[\"GRS 1980\",6378137,298.257222101004,AUTHORITY[\"EPSG\",\"7019\"]],TOWGS84[0,0,0,0,0,0,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",173],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",1600000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]")), GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}([(1.1253763259e6, 4.982869816e6), (1.1248143323e6, 4.976641044e6), (1.1282897921e6, 4.9707170423e6), (1.1254201958e6, 4.9665221375e6), (1.1188173258e6, 4.9720125867e6), (1.1253763259e6, 4.982869816e6)], nothing, WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "PROJCS[\"unknown\",GEOGCS[\"unknown\",DATUM[\"Unknown based on GRS80 ellipsoid\",SPHEROID[\"GRS 1980\",6378137,298.257222101004,AUTHORITY[\"EPSG\",\"7019\"]],TOWGS84[0,0,0,0,0,0,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",173],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",1600000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]")), GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}([(1.7015121629e6, 5.9962048743e6), (1.6940719249e6, 5.9966701668e6), (1.6984728948e6, 5.9883324659e6), (1.7037685647e6, 5.9848085328e6), (1.7060307416e6, 5.9751359042e6), (1.7018195613e6, 5.9719674757e6), (1.6933144403e6, 5.9719708739e6), (1.6883953701e6, 5.9896812008e6), (1.6779656958e6, 6.0055889233e6), (1.662101242e6, 6.0277652897e6) โฆ (1.725607508e6, 5.945892267e6), (1.7281682052e6, 5.9584185318e6), (1.7268710844e6, 5.9703482785e6), (1.7127870833e6, 5.9730682541e6), (1.7128762566e6, 5.9791634373e6), (1.7232886065e6, 5.9797407952e6), (1.7331650316e6, 5.9897140167e6), (1.7201965583e6, 5.9800776907e6), (1.709110421e6, 5.9866719829e6), (1.7015121629e6, 5.9962048743e6)], nothing, WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "PROJCS[\"unknown\",GEOGCS[\"unknown\",DATUM[\"Unknown based on GRS80 ellipsoid\",SPHEROID[\"GRS 1980\",6378137,298.257222101004,AUTHORITY[\"EPSG\",\"7019\"]],TOWGS84[0,0,0,0,0,0,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",173],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",1600000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]")), GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}([(1.6650793317e6, 5.4762796377e6), (1.6679650479e6, 5.4828286371e6), (1.6790721723e6, 5.4888496631e6), (1.6808040605e6, 5.4949012903e6), (1.6810365469e6, 5.4869445483e6), (1.6763032282e6, 5.4766970825e6), (1.6685295209e6, 5.4686427381e6), (1.6650793317e6, 5.4762796377e6)], nothing, WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "PROJCS[\"unknown\",GEOGCS[\"unknown\",DATUM[\"Unknown based on GRS80 ellipsoid\",SPHEROID[\"GRS 1980\",6378137,298.257222101004,AUTHORITY[\"EPSG\",\"7019\"]],TOWGS84[0,0,0,0,0,0,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",173],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",1600000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]")), GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}([(1.7939818643e6, 5.9319178426e6), (1.7912570345e6, 5.9199068808e6), (1.7842798937e6, 5.9215788138e6), (1.7818307229e6, 5.9258456737e6), (1.7939818643e6, 5.9319178426e6)], nothing, WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "PROJCS[\"unknown\",GEOGCS[\"unknown\",DATUM[\"Unknown based on GRS80 ellipsoid\",SPHEROID[\"GRS 1980\",6378137,298.257222101004,AUTHORITY[\"EPSG\",\"7019\"]],TOWGS84[0,0,0,0,0,0,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",173],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",1600000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]")), GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, WellKnownText{GeoFormatTypes.CRS}}([(1.8266292121e6, 5.9745333354e6), (1.8188737655e6, 5.9791838569e6), (1.8185878044e6, 5.9830890239e6), (1.8112269977e6, 5.9881418267e6), (1.8126655001e6, 5.9932466009e6), (1.8093053687e6, 5.9982083619e6), (1.8108957047e6, 6.0059053369e6), (1.8159541251e6, 6.0078779063e6), (1.8180162185e6, 5.9990038299e6), (1.823326993e6, 5.9918794741e6), (1.8227997093e6, 5.9860917836e6), (1.8283790056e6, 5.9803071686e6), (1.8266292121e6, 5.9745333354e6)], nothing, WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "PROJCS[\"unknown\",GEOGCS[\"unknown\",DATUM[\"Unknown based on GRS80 ellipsoid\",SPHEROID[\"GRS 1980\",6378137,298.257222101004,AUTHORITY[\"EPSG\",\"7019\"]],TOWGS84[0,0,0,0,0,0,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",173],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",1600000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]"))], nothing, WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "PROJCS[\"unknown\",GEOGCS[\"unknown\",DATUM[\"Unknown based on GRS80 ellipsoid\",SPHEROID[\"GRS 1980\",6378137,298.257222101004,AUTHORITY[\"EPSG\",\"7019\"]],TOWGS84[0,0,0,0,0,0,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",173],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",1600000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]"))For a โtemplateโ raster, we will aggregate the New Zealand DEM, in the nz_elev.tif file, to 5 times coarser resolution. The code section below follows the aggeregation example in ?sec-raster-agg-disagg.
factor = 2/10
r = Rasters.resample(src_nz_elev; size = round.(Int, size(src_nz_elev) .* factor), method = :average)โญโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโฎ
โ 223ร290 Raster{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 missingThe 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).
fig, ax, plt = plot(r)
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.
dp = DimPoints(r)
nonmissing_points = dp[r .=== missingval(r)]53316-element Vector{Tuple{Float64, Float64}}:
(995456.5461976258, 6.186960927303582e6)
(1.0004565461976258e6, 6.186960927303582e6)
(1.0054565461976257e6, 6.186960927303582e6)
(1.0104565461976259e6, 6.186960927303582e6)
(1.0154565461976258e6, 6.186960927303582e6)
(1.0204565461976258e6, 6.186960927303582e6)
(1.0254565461976258e6, 6.186960927303582e6)
(1.0304565461976257e6, 6.186960927303582e6)
(1.0354565461976259e6, 6.186960927303582e6)
(1.0404565461976258e6, 6.186960927303582e6)
โฎ
(2.0654565461976258e6, 4.741960927303582e6)
(2.0704565461976258e6, 4.741960927303582e6)
(2.0754565461976258e6, 4.741960927303582e6)
(2.0804565461976258e6, 4.741960927303582e6)
(2.0854565461976258e6, 4.741960927303582e6)
(2.0904565461976258e6, 4.741960927303582e6)
(2.0954565461976258e6, 4.741960927303582e6)
(2.100456546197626e6, 4.741960927303582e6)
(2.105456546197626e6, 4.741960927303582e6)The result is a vector of 2-tuples, which are recognized as GeoInterface point geometries.
We can compute the Cartesian distance from each point to the nearest line in the coastline multilinestring using the distance method from GeometryOps.
distances = GO.distance.((coastline,), nonmissing_points)53316-element Vector{Float64}:
572764.2876244619
567764.3239539921
562764.3609290748
557764.3985670706
552764.436885969
547764.4759044152
542764.5156417422
537764.5561180016
532764.5973539978
527764.639371324
โฎ
610809.9189554364
614562.8629679172
618333.460184955
622121.3896252689
625926.3366197749
629747.9927075815
633586.0555316888
637440.2287345697
641310.2218537896Finally, we rasterize (see Section 6.4.1) the distances into our raster template.
img = Rasters.rasterize(
last,
nonmissing_points;
to = r,
fill = distances,
)โญโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโฎ
โ 223ร290 Raster{Union{Missing, Float64},2} last โ
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโดโโโโโโโโโโโโโโโโโโโโโโโ dims โ
โ X Projected{Float64} LinRange{Float64}(995456.5461976258, 2.105456546197626e6, 223) ForwardOrdered Regular Intervals{Start},
โ Y Projected{Float64} LinRange{Float64}(6.186960927303582e6, 4.741960927303582e6, 290) ReverseOrdered Regular Intervals{Start}
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโ metadata โค
Metadata of Dict{Any, Any}()
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโ raster โค
extent: Extent(X = (995456.5461976258, 2.110456546197626e6), Y = (4.741960927303582e6, 6.191960927303582e6))
missingval: missing
crs: PROJCS["unknown",GEOGCS["unknown",DATUM["Unknown based on GRS80 ellipsoid",SPHEROID["GRS 1980",6378137,298.257222101004,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",173],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",1600000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโ
โ โ 6.18696e6 6.18196e6 โฆ 4.74696e6 4.74196e6
9.95457e5 5.72764e5 5.72767e5 164188.0 1.67786e5
1.00046e6 5.67764e5 5.67767e5 1.60717e5 164391.0
1.00546e6 5.62764e5 5.62767e5 1.57329e5 161080.0
โฎ โฑ โฎ
2.09546e6 3.31952e5 3.29282e5 6.30403e5 6.33586e5
2.10046e6 3.36173e5 3.33536e5 6.34276e5 6.3744e5
2.10546e6 340415.0 3.37811e5 6.38166e5 6.4131e5The final result, a raster of distances to the nearest coastline, is shown in Figure 6.10.
fig, ax, plt = plot(img)
lines!(ax, coastline; color = :red)
Colorbar(fig[1, 2], plt; label = "Distance to coastline (m)")
fig