This chapter outlines two fundamental geographic data models โ vector and raster โ and introduces the main Python packages for working with them. Before demonstrating their implementation in Python, we will introduce the theory behind each data model and the disciplines in which they predominate.
The vector data model (Section 1.2) represents the world using points, lines, and polygons. These have discrete, well-defined borders, meaning that vector datasets usually have a high level of precision (but not necessarily accuracy). The raster data model (?sec-raster-data), on the other hand, divides the surface up into cells of constant size. Raster datasets are the basis of background images used in web-mapping and have been a vital source of geographic data since the origins of aerial photography and satellite-based remote sensing devices. Rasters aggregate spatially specific features to a given resolution, meaning that they are consistent over space and scalable, with many worldwide raster datasets available.
Which to use? The answer likely depends on your domain of application, and the datasets you have access to:
Vector datasets and methods dominate the social sciences because human settlements and and processes (e.g., transport infrastructure) tend to have discrete borders.
Raster datasets and methods dominate many environmental sciences because of the reliance on remote sensing data.
Julia has strong support for both data models. We will focus on GeoDataFrames.jl and the GeoInterface.jl ecosystem for working with vector data, including the packages GeometryOps.jl and LibGEOS.jl. We will focus on the Rasters.jl package for working with rasters.
TODO: alternatives, geostats, etc.
There is much overlap in some fields and raster and vector datasets can be used together: ecologists and demographers, for example, commonly use both vector and raster data. Furthermore, it is possible to convert between the two forms (see Chapter 5). Whether your work involves more use of vector or raster datasets, it is worth understanding the underlying data models before using them, as discussed in subsequent chapters.
1.2 Vector data
The geographic vector data model is based on points located within a coordinate reference system (CRS). Points can represent self-standing features (e.g., the location of a bus stop), or they can be linked together to form more complex geometries such as lines and polygons. Most point geometries contain only two dimensions (3-dimensional CRSs may contain an additional \(z\) value, typically representing height above sea level).
In this system, London, for example, can be represented by the coordinates (-0.1,51.5). This means that its location is -0.1 degrees east and 51.5 degrees north of the origin. The origin, in this case, is at 0 degrees longitude (a prime meridian located at Greenwich) and 0 degrees latitude (the Equator) in a geographic (โlon/latโ) CRS (Figure 1.1, left panel). The same point could also be approximated in a projected CRS with โEasting/Northingโ values of (530000,180000) in the British National Grid, meaning that London is located 530 \(km\) East and 180 \(km\) North of the origin of the CRS (Figure 1.1, right panel). The location of National Gridโs origin, in the sea beyond South West Peninsular, ensures that most locations in the UK have positive Easting and Northing values.
Figure 1.1: Illustration of vector (point) data in which location of London (the red X) is represented with reference to an origin (the blue circle). The left plot represents a geographic CRS with an origin at 0ยฐ longitude and latitude. The right plot represents a projected CRS with an origin located in the sea west of the South West Peninsula.
There is more to CRSs, as described in ?sec-coordinate-reference-systems-intro and Chapter 6 but, for the purposes of this section, it is sufficient to know that coordinates consist of two numbers representing the distance from an origin, usually in \(x\) and \(y\) dimensions.
TODO: explain the JuliaGeo ecosystem like they explain geopandas E.g GeoInterface defines how to access any geometry, then LibGEOS (wrapping GEOS), GeometryOps, Proj, etc consume such geometries.
1.2.1 Vector data classes
Juliaโs geographic vector data model is based on the Simple Features standard, which is an ISO standard for representing vector data. Simple Features defines types of geometries (points, lines, polygons, multipolygons, etc.), as well as โfeature collectionsโ that are basically tables of geometries associated with some data.
Starting with the highest level class, feature collections come in two flavours: - Loaded from file (Shapefile.Table, GeoJSON.FeatureCollection, โฆ) - Tables with geometry columns (e.g. DataFrame, but can be any Julia table)
One can easily convert from feature collections to
Precompiling GeoDataFrames...
515.5 ms โ IteratorInterfaceExtensions
566.7 ms โ LaTeXStrings
530.8 ms โ CEnum
573.4 ms โ TensorCore
687.2 ms โ GeoFormatTypes
660.4 ms โ Observables
719.9 ms โ Statistics
721.8 ms โ Requires
1108.2 ms โ OrderedCollections
482.3 ms โ DataValueInterfaces
574.3 ms โ Reexport
1635.2 ms โ OffsetArrays
570.4 ms โ InvertedIndices
682.5 ms โ Extents
998.6 ms โ InlineStrings
826.8 ms โ IntervalSets
642.4 ms โ LRUCache
517.1 ms โ StaticArraysCore
622.9 ms โ DataAPI
650.8 ms โ MappedArrays
777.7 ms โ MPIPreferences
814.3 ms โ OpenSSL_jll
739.9 ms โ Bzip2_jll
1383.4 ms โ Qhull_jll
1972.4 ms โ boost_jll
1303.9 ms โ Lz4_jll
2609.6 ms โ SentinelArrays
1025.4 ms โ ICU_jll
968.8 ms โ libaec_jll
785.1 ms โ libpng_jll
906.2 ms โ Xorg_libXau_jll
948.6 ms โ Hwloc_jll
805.3 ms โ Giflib_jll
868.6 ms โ LERC_jll
933.5 ms โ SQLite_jll
935.7 ms โ EarCut_jll
930.5 ms โ JpegTurbo_jll
783.0 ms โ snappy_jll
927.3 ms โ XZ_jll
832.8 ms โ Xorg_libXdmcp_jll
894.9 ms โ brotli_jll
797.8 ms โ Expat_jll
860.4 ms โ MicrosoftMPI_jll
976.7 ms โ Zstd_jll
746.8 ms โ Xorg_xtrans_jll
870.0 ms โ GEOS_jll
790.4 ms โ Libgpg_error_jll
788.8 ms โ Kerberos_krb5_jll
1026.1 ms โ Xorg_libpthread_stubs_jll
486.6 ms โ TableTraits
904.5 ms โ XML2_jll
539.8 ms โ StackViews
702.7 ms โ PaddedViews
2865.3 ms โ RecipesBase
1089.6 ms โ IntervalSets โ IntervalSetsRandomExt
806.1 ms โ IntervalSets โ IntervalSetsStatisticsExt
3665.4 ms โ StringManipulation
3162.3 ms โ DataStructures
1064.6 ms โ PooledArrays
1942.2 ms โ DiskArrays
4900.5 ms โ FixedPointNumbers
689.9 ms โ Missings
1979.7 ms โ GeoInterface
954.8 ms โ Thrift_jll
1770.0 ms โ MPItrampoline_jll
1715.7 ms โ OpenMPI_jll
1168.7 ms โ HDF4_jll
1849.6 ms โ MPICH_jll
863.6 ms โ Blosc_jll
1115.1 ms โ libzip_jll
982.2 ms โ Libtiff_jll
854.0 ms โ LibPQ_jll
1154.2 ms โ Libgcrypt_jll
604.4 ms โ MosaicViews
888.5 ms โ IntervalSets โ IntervalSetsRecipesBaseExt
875.9 ms โ SortingAlgorithms
1432.7 ms โ Tables
858.8 ms โ GeoInterfaceRecipes
998.8 ms โ Arrow_jll
1074.9 ms โ LittleCMS_jll
1315.4 ms โ PROJ_jll
2945.9 ms โ HDF5_jll
3660.6 ms โ ColorTypes
1104.7 ms โ XSLT_jll
918.5 ms โ OpenJpeg_jll
1269.0 ms โ libgeotiff_jll
2312.5 ms โ NetCDF_jll
13683.7 ms โ StaticArrays
2219.9 ms โ Xorg_libxcb_jll
4351.9 ms โ ColorVectorSpace
1163.3 ms โ StaticArrays โ StaticArraysStatisticsExt
981.1 ms โ Xorg_libX11_jll
770.2 ms โ Xorg_libXext_jll
1041.2 ms โ Libglvnd_jll
1362.8 ms โ libwebp_jll
8277.4 ms โ Colors
2753.0 ms โ GDAL_jll
4962.6 ms โ GDAL
19233.2 ms โ GeometryBasics
31036.0 ms โ PrettyTables
6832.2 ms โ MakieCore
1699.5 ms โ GeoInterfaceMakie
27289.7 ms โ ImageCore
11887.2 ms โ ArchGDAL
44723.6 ms โ DataFrames
2628.0 ms โ GeoDataFrames
106 dependencies successfully precompiled in 99 seconds. 37 already precompiled.
177ร11 DataFrame
152 rows omitted
Row
geom
iso_a2
name_long
continent
region_un
subregion
type
area_km2
pop
lifeExp
gdpPercap
IGeometrโฆ
String?
String
String
String
String
String
Float64
Float64?
Float64?
Float64?
1
Geometry: wkbMultiPolygon
FJ
Fiji
Oceania
Oceania
Melanesia
Sovereign country
19290.0
885806.0
69.96
8222.25
2
Geometry: wkbMultiPolygon
TZ
Tanzania
Africa
Africa
Eastern Africa
Sovereign country
9.32746e5
5.22349e7
64.163
2402.1
3
Geometry: wkbMultiPolygon
EH
Western Sahara
Africa
Africa
Northern Africa
Indeterminate
96270.6
missing
missing
missing
4
Geometry: wkbMultiPolygon
CA
Canada
North America
Americas
Northern America
Sovereign country
1.0036e7
3.55353e7
81.953
43079.1
5
Geometry: wkbMultiPolygon
US
United States
North America
Americas
Northern America
Country
9.51074e6
3.18623e8
78.8415
51922.0
6
Geometry: wkbMultiPolygon
KZ
Kazakhstan
Asia
Asia
Central Asia
Sovereign country
2.72981e6
1.72883e7
71.62
23587.3
7
Geometry: wkbMultiPolygon
UZ
Uzbekistan
Asia
Asia
Central Asia
Sovereign country
4.6141e5
3.07577e7
71.039
5370.87
8
Geometry: wkbMultiPolygon
PG
Papua New Guinea
Oceania
Oceania
Melanesia
Sovereign country
4.6452e5
7.75578e6
65.23
3709.08
9
Geometry: wkbMultiPolygon
ID
Indonesia
Asia
Asia
South-Eastern Asia
Sovereign country
1.81925e6
2.55131e8
68.856
10003.1
10
Geometry: wkbMultiPolygon
AR
Argentina
South America
Americas
South America
Sovereign country
2.78447e6
4.29815e7
76.252
18797.5
11
Geometry: wkbMultiPolygon
CL
Chile
South America
Americas
South America
Sovereign country
8.14844e5
1.76138e7
79.117
22195.3
12
Geometry: wkbMultiPolygon
CD
Democratic Republic of the Congo
Africa
Africa
Middle Africa
Sovereign country
2.32349e6
7.37229e7
58.782
785.347
13
Geometry: wkbMultiPolygon
SO
Somalia
Africa
Africa
Eastern Africa
Sovereign country
4.84333e5
1.35131e7
55.467
missing
โฎ
โฎ
โฎ
โฎ
โฎ
โฎ
โฎ
โฎ
โฎ
โฎ
โฎ
โฎ
166
Geometry: wkbMultiPolygon
ET
Ethiopia
Africa
Africa
Eastern Africa
Sovereign country
1.13239e6
9.73668e7
64.535
1424.53
167
Geometry: wkbMultiPolygon
DJ
Djibouti
Africa
Africa
Eastern Africa
Sovereign country
21880.3
912164.0
62.006
missing
168
Geometry: wkbMultiPolygon
missing
Somaliland
Africa
Africa
Eastern Africa
Indeterminate
1.6735e5
missing
missing
missing
169
Geometry: wkbMultiPolygon
UG
Uganda
Africa
Africa
Eastern Africa
Sovereign country
2.45768e5
3.88333e7
59.224
1637.28
170
Geometry: wkbMultiPolygon
RW
Rwanda
Africa
Africa
Eastern Africa
Sovereign country
23365.4
1.13454e7
66.188
1629.87
171
Geometry: wkbMultiPolygon
BA
Bosnia and Herzegovina
Europe
Europe
Southern Europe
Sovereign country
50605.1
3.566e6
76.561
10516.8
172
Geometry: wkbMultiPolygon
MK
Macedonia
Europe
Europe
Southern Europe
Sovereign country
25062.3
2.0775e6
75.384
12298.5
173
Geometry: wkbMultiPolygon
RS
Serbia
Europe
Europe
Southern Europe
Sovereign country
76388.6
7.13058e6
75.3366
13112.9
174
Geometry: wkbMultiPolygon
ME
Montenegro
Europe
Europe
Southern Europe
Sovereign country
13443.7
621810.0
76.712
14796.6
175
Geometry: wkbMultiPolygon
XK
Kosovo
Europe
Europe
Southern Europe
Sovereign country
11230.3
1.8218e6
71.0976
8698.29
176
Geometry: wkbMultiPolygon
TT
Trinidad and Tobago
North America
Americas
Caribbean
Sovereign country
7737.81
1.35449e6
70.426
31181.8
177
Geometry: wkbMultiPolygon
SS
South Sudan
Africa
Africa
Eastern Africa
Sovereign country
6.24909e5
1.1531e7
55.817
1935.88
usingCairoMakieusingGeoMakief, a, p =poly(df.geom)
Precompiling CairoMakie...
466.2 ms โ RangeArrays
528.3 ms โ IndirectArrays
563.2 ms โ PolygonOps
606.5 ms โ StatsAPI
765.3 ms โ AbstractFFTs
589.8 ms โ Contour
708.2 ms โ TriplotBase
527.6 ms โ InverseFunctions
593.3 ms โ EnumX
1755.9 ms โ FillArrays
1902.9 ms โ Format
750.8 ms โ StableRNGs
1414.5 ms โ Grisu
782.4 ms โ AbstractTrees
614.8 ms โ RoundingEmulator
628.6 ms โ PtrArrays
751.5 ms โ NaNMath
971.9 ms โ TranscodingStreams
2318.3 ms โ IrrationalConstants
3926.9 ms โ MacroTools
923.6 ms โ LazyModules
863.7 ms โ Ratios
716.9 ms โ Inflate
635.9 ms โ ConstructionBase
662.2 ms โ Adapt
2587.1 ms โ AdaptivePredicates
951.3 ms โ Statistics โ SparseArraysExt
784.9 ms โ SuiteSparse
970.3 ms โ WoodburyMatrices
950.0 ms โ DocStringExtensions
1412.9 ms โ ProgressMeter
917.3 ms โ Animations
596.6 ms โ SignedDistanceFields
886.8 ms โ Graphite2_jll
763.2 ms โ Libmount_jll
2170.6 ms โ ChainRulesCore
879.8 ms โ LLVMOpenMP_jll
2335.9 ms โ UnicodeFun
1167.6 ms โ Rmath_jll
724.2 ms โ Imath_jll
1008.1 ms โ libfdk_aac_jll
863.8 ms โ LAME_jll
825.8 ms โ CRlibm_jll
905.5 ms โ Ogg_jll
1135.8 ms โ oneTBB_jll
951.7 ms โ x265_jll
1022.9 ms โ x264_jll
837.1 ms โ libaom_jll
819.1 ms โ LZO_jll
936.6 ms โ Opus_jll
1064.6 ms โ Libffi_jll
1007.7 ms โ isoband_jll
839.2 ms โ FFTW_jll
987.6 ms โ OpenSpecFun_jll
823.7 ms โ Libuuid_jll
950.2 ms โ FriBidi_jll
997.4 ms โ libsixel_jll
771.7 ms โ OpenBLASConsistentFPCSR_jll
1711.7 ms โ PkgVersion
1948.4 ms โ FilePathsBase
1326.2 ms โ FreeType2_jll
2525.8 ms โ JSON
6101.1 ms โ ColorSchemes
1795.8 ms โ QuadGK
7848.5 ms โ FileIO
2462.0 ms โ IntelOpenMP_jll
4520.0 ms โ ImageBase
2325.1 ms โ Packing
2283.9 ms โ ShaderAbstractions
1086.0 ms โ Xorg_libXrender_jll
1219.3 ms โ Gettext_jll
1646.9 ms โ AxisArrays
16391.8 ms โ SIMD
798.1 ms โ InverseFunctions โ InverseFunctionsTestExt
594.8 ms โ InverseFunctions โ InverseFunctionsDatesExt
2423.8 ms โ AbstractFFTs โ AbstractFFTsTestExt
833.1 ms โ FillArrays โ FillArraysStatisticsExt
856.1 ms โ Showoff
1672.4 ms โ FillArrays โ FillArraysSparseArraysExt
1263.4 ms โ AliasTables
1500.0 ms โ Graphics
1049.7 ms โ Ratios โ RatiosFixedPointNumbersExt
2036.0 ms โ SimpleTraits
1088.8 ms โ ConstructionBase โ ConstructionBaseIntervalSetsExt
739.9 ms โ ConstructionBase โ ConstructionBaseLinearAlgebraExt
9543.9 ms โ PNGFiles
14592.1 ms โ GridLayoutBase
1357.8 ms โ ConstructionBase โ ConstructionBaseStaticArraysExt
1335.8 ms โ Adapt โ AdaptSparseArraysExt
1534.5 ms โ Adapt โ AdaptStaticArraysExt
815.2 ms โ OffsetArrays โ OffsetArraysAdaptExt
1168.3 ms โ AxisAlgorithms
1185.9 ms โ LogExpFunctions
1574.5 ms โ PDMats
1232.9 ms โ ChainRulesCore โ ChainRulesCoreSparseArraysExt
878.5 ms โ AbstractFFTs โ AbstractFFTsChainRulesCoreExt
1055.4 ms โ Pixman_jll
1435.9 ms โ StaticArrays โ StaticArraysChainRulesCoreExt
1524.0 ms โ Rmath
884.5 ms โ Isoband
1399.3 ms โ OpenEXR_jll
1194.7 ms โ libvorbis_jll
783.2 ms โ FilePathsBase โ FilePathsBaseMmapExt
1432.1 ms โ FilePaths
2179.1 ms โ FilePathsBase โ FilePathsBaseTestExt
1917.9 ms โ FreeType
1524.9 ms โ Fontconfig_jll
1166.0 ms โ ColorBrewer
40287.9 ms โ Unitful
6387.1 ms โ IntervalArithmetic
2444.7 ms โ QOI
4392.7 ms โ Sixel
5098.7 ms โ JpegTurbo
1236.9 ms โ Glib_jll
4164.5 ms โ WebP
6315.4 ms โ MKL_jll
2882.4 ms โ ImageAxes
789.1 ms โ LogExpFunctions โ LogExpFunctionsInverseFunctionsExt
1157.0 ms โ StructArrays
2698.4 ms โ LogExpFunctions โ LogExpFunctionsChainRulesCoreExt
5009.9 ms โ SpecialFunctions
1206.6 ms โ FillArrays โ FillArraysPDMatsExt
4316.5 ms โ StatsBase
19859.1 ms โ PlotUtils
2630.5 ms โ OpenEXR
3812.6 ms โ Interpolations
1204.8 ms โ Unitful โ ConstructionBaseUnitfulExt
1013.7 ms โ Unitful โ InverseFunctionsUnitfulExt
1047.1 ms โ IntervalArithmetic โ IntervalArithmeticIntervalSetsExt
4620.6 ms โ FreeTypeAbstraction
1431.7 ms โ Cairo_jll
2229.5 ms โ ImageMetadata
817.9 ms โ StructArrays โ StructArraysAdaptExt
1205.5 ms โ StructArrays โ StructArraysSparseArraysExt
1523.1 ms โ StructArrays โ StructArraysStaticArraysExt
1101.0 ms โ StructArrays โ StructArraysLinearAlgebraExt
26784.1 ms โ Automa
10395.2 ms โ ExactPredicates
981.8 ms โ ColorVectorSpace โ SpecialFunctionsExt
1824.5 ms โ HypergeometricFunctions
1317.5 ms โ HarfBuzz_jll
3543.0 ms โ SpecialFunctions โ SpecialFunctionsChainRulesCoreExt
2460.6 ms โ Interpolations โ InterpolationsUnitfulExt
3997.3 ms โ Netpbm
3716.2 ms โ StatsFuns
1378.4 ms โ libass_jll
1229.7 ms โ Pango_jll
21400.8 ms โ FFTW
10329.3 ms โ DelaunayTriangulation
961.1 ms โ StatsFuns โ StatsFunsInverseFunctionsExt
12504.4 ms โ MathTeXEngine
2342.5 ms โ StatsFuns โ StatsFunsChainRulesCoreExt
1622.0 ms โ FFMPEG_jll
1642.9 ms โ Cairo
6148.3 ms โ Distributions
1386.5 ms โ Distributions โ DistributionsChainRulesCoreExt
1458.5 ms โ Distributions โ DistributionsTestExt
1332.9 ms โ KernelDensity
84424.7 ms โ TiffImages
1013.2 ms โ ImageIO
131731.8 ms โ Makie
41970.0 ms โ CairoMakie
162 dependencies successfully precompiled in 310 seconds. 110 already precompiled.
1 dependency had output during precompilation:
โ MKL_jll
โ Downloading artifact: IntelOpenMP
โ
Precompiling ParsersExt...
410.6 ms โ InlineStrings โ ParsersExt
1 dependency successfully precompiled in 0 seconds. 9 already precompiled.
Precompiling SerializationExt...
365.6 ms โ LRUCache โ SerializationExt
1 dependency successfully precompiled in 0 seconds. 2 already precompiled.
Precompiling IntervalArithmeticRecipesBaseExt...
671.5 ms โ IntervalArithmetic โ IntervalArithmeticRecipesBaseExt
1 dependency successfully precompiled in 1 seconds. 33 already precompiled.
Precompiling ArchGDALMakieExt...
2153.4 ms โ ArchGDAL โ ArchGDALMakieExt
1 dependency successfully precompiled in 3 seconds. 306 already precompiled.
Precompiling GeoMakie...
568.6 ms โ SortTileRecursiveTree
879.2 ms โ GeometryOpsCore
1108.9 ms โ CoordinateTransformations
2499.4 ms โ Geodesy
2530.7 ms โ Proj
4161.5 ms โ GeometryOps
1196.4 ms โ GeometryOps โ GeometryOpsProjExt
25375.2 ms โ GeoJSON
1122.9 ms โ NaturalEarth
1343.3 ms โ GeoJSON โ GeoJSONMakieExt
8964.5 ms โ GeoMakie
11 dependencies successfully precompiled in 36 seconds. 276 already precompiled.
1.2.2 Raster from scratch
In this section, we are going to demonstrate the creation of rasters from scratch. We will construct two small rasters, elev and grain, which we will use in examples later in the book. Unlike creating a vector layer (see ?sec-vector-layer-from-scratch), creating a raster from scratch is rarely needed in practice because aligning a raster with the proper spatial extent is challenging to do programmatically (โgeoreferencingโ tools in GIS software are a better fit for the job). Nevertheless, the examples will be helpful to become more familiar with the Rasters.jl data structures.
usingRastersimportGeoFormatTypes as GFT
Precompiling Rasters...
627.7 ms โ FieldMetadata
632.6 ms โ Interfaces
811.6 ms โ ArrayInterface
621.4 ms โ Flatten
1239.0 ms โ CFTime
511.1 ms โ ArrayInterface โ ArrayInterfaceStaticArraysCoreExt
816.5 ms โ ArrayInterface โ ArrayInterfaceSparseArraysExt
1665.3 ms โ Setfield
1311.6 ms โ CommonDataModel
14662.1 ms โ DimensionalData
772.7 ms โ DimensionalData โ DimensionalDataDiskArraysExt
5468.8 ms โ Rasters
12 dependencies successfully precompiled in 23 seconds. 58 already precompiled.
Precompiling ArrayInterfaceChainRulesCoreExt...
370.8 ms โ ArrayInterface โ ArrayInterfaceChainRulesCoreExt
1 dependency successfully precompiled in 0 seconds. 11 already precompiled.
Precompiling DimensionalDataMakie...
813.2 ms โ DimensionalData โ DimensionalDataStatsBase
6630.9 ms โ DimensionalData โ DimensionalDataMakie
2 dependencies successfully precompiled in 7 seconds. 277 already precompiled.
Precompiling RastersMakieExt...
1301.8 ms โ Rasters โ RastersStatsBaseExt
6446.6 ms โ Rasters โ RastersMakieExt
2 dependencies successfully precompiled in 7 seconds. 291 already precompiled.
Precompiling RastersProjExt...
1279.5 ms โ Rasters โ RastersCoordinateTransformationsExt
1377.6 ms โ Rasters โ RastersProjExt
2 dependencies successfully precompiled in 1 seconds. 92 already precompiled.
Precompiling RastersArchGDALExt...
2315.4 ms โ Rasters โ RastersArchGDALExt
1 dependency successfully precompiled in 3 seconds. 165 already precompiled.
Conceptually, a raster is an array combined with georeferencing information, whereas the latter comprises:
Lookup vectors for the axes, encoding the spatial coordinates for each grid cell. These take the form of the X and Y dimensions in the raster that youโll see below.
A coordinate reference system (CRS) definition, specifying the association of the rasterโs coordinates with the surface of the earth.
Therefore, to create a raster, we first need to have an array with the values, and then supplement it with the georeferencing information. Letโs create the arrays elev and grain. The elev array is a \(6 \times 6\) array with sequential values from 1 to 36. It can be created as follows using base Julia functions.
The grain array represents a categorical raster with values 0, 1, 2, corresponding to categories โclayโ, โsiltโ, โsandโ, respectively. We will create it from a specific arrangement of pixel values, using reshape.
Note that in both cases, we are using the uint8 (unsigned integer in 8 bits, i.e., 0-255) data type, which is sufficient to represent all possible values of the given rasters (see ?tbl-numpy-data-types). This is the recommended approach for a minimal memory footprint.
What is missing now is the georeferencing information (see ?sec-using-rasters-jl). In this case, since the rasters are arbitrary, we also set up arbitrary dimension lookups for the x and y axes, where:
The origin (\(x_{min}\), \(y_{max}\)) is at -1.5,1.5
The raster resolution (\(delta_{x}\), \(delta_{y}\)) is 0.5,-0.5
We can add this information using rasterio.transform.from_origin, and specifying west, north, xsize, and ysize parameters. The resulting transformation matrix object is hereby named new_transform.
We can now construct a Raster object, from the elev array and the dimensions new_x and new_y.
We assign to it a CRS of EPSG:4326 (which encodes that the coordinate system is longitude/latitude on the โstandardโ WGS84 definition of the Earthโs curvature).
Here, we use the GFT.EPSG(code) constructor to create an object that encodes a reference code under the European Petroleum Survey Group (EPSG) authorityโs database of coordinate reference systems.
The raster can now be plotted in its coordinate system, passing the array elev along with the transformation matrix new_transform to rasterio.plot.show (Figure 1.2).
plot(elev_raster)
Figure 1.2: Plot of the elev raster, a minimal example of a continuous raster, created from scratch
The grain raster can be plotted the same way, as we are going to use the same transformation matrix for it as well (Figure 1.3).
Figure 1.3: Plot of the grain raster, a minimal example of a categorical raster, created from scratch
At this point, we have two rasters, each composed of an array and related dimension lookups. We can work with the raster using Rasters.jl by:
Keeping in mind that any other layer we use in the analysis is in the same CRS
Finally, to export the raster for permanent storage, along with the spatial metadata, we need to go through the following steps:
Create a raster file (where we set the lookups and the CRS, among other settings)
Write the array with raster values into the connection
Close the connection
Donโt worry if the code below is unclear; the concepts related to writing raster data to file will be explained in ?sec-data-output-raster. For now, for completeness, and also to use these rasters in subsequent chapters without having to re-create them from scratch, we just provide the code for exporting the elev and grain rasters into the output directory. In the case of elev, we do it as follows with the Rasters.write functions and methods of the Rasters.jl package.
write("output/elev.tif", elev_raster; force =true)
"output/elev.tif"
Note that the CRS we (arbitrarily) set for the elev raster is WGS84, defined using crs=4326 according to the EPSG code.
Exporting the grain raster is done in the same way, with the only differences being the file name and the array we write into the connection.
write("output/grain.tif", Raster(grain, (new_x, new_y); crs = GFT.EPSG(4326)); force =true)
"output/grain.tif"
As a result, the files elev.tif and grain.tif are written into the output directory. We are going to use these small raster files later on in the examples (for example, ?sec-raster-subsetting).
Note that the transform matrices and dimensions of elev and grain are identical. This means that the rasters are overlapping, and can be combined into one two-band raster, processed in raster algebra operations (?sec-map-algebra), etc.
