1  Geographic data in Julia

using Pkg
Pkg.status()

Status `~/work/geocompjl/geocompjl/Project.toml`
  [c9ce4bd3] ArchGDAL v0.10.8
  [13f3f980] CairoMakie v0.13.2
  [a93c6f00] DataFrames v1.7.0
  [0703355e] DimensionalData v0.29.13
  [62cb38b5] GeoDataFrames v0.3.12
  [68eda718] GeoFormatTypes v0.4.4
  [cf35fbd7] GeoInterface v1.4.1
  [61d90e0f] GeoJSON v0.8.2
  [db073c08] GeoMakie v0.7.11
  [3251bfac] GeometryOps v0.1.15
  [a90b1aa1] LibGEOS v0.9.4
  [c94c279d] Proj v1.8.1
  [a3a2b9e3] Rasters v0.14.2
  [10745b16] Statistics v1.11.1
  [2913bbd2] StatsBase v0.34.4

"output"

1.1 Introduction

using GeoDataFrames
df = GeoDataFrames.read("data/world.gpkg")

Precompiling GeoDataFrames...
    451.6 ms  ✓ IteratorInterfaceExtensions
    561.9 ms  ✓ LaTeXStrings
    526.2 ms  ✓ GeoFormatTypes
    566.3 ms  ✓ CEnum
    606.6 ms  ✓ TensorCore
    609.7 ms  ✓ Requires
    736.7 ms  ✓ Statistics
    731.9 ms  ✓ Observables
    782.7 ms  ✓ OrderedCollections
    436.0 ms  ✓ DataValueInterfaces
    472.8 ms  ✓ Reexport
    776.8 ms  ✓ IntervalSets
    903.1 ms  ✓ InlineStrings
    546.8 ms  ✓ InvertedIndices
   1863.4 ms  ✓ OffsetArrays
    703.8 ms  ✓ Extents
    445.5 ms  ✓ DataAPI
    522.8 ms  ✓ StaticArraysCore
    524.8 ms  ✓ MappedArrays
    915.5 ms  ✓ LRUCache
    884.7 ms  ✓ MPIPreferences
    945.5 ms  ✓ OpenSSL_jll
   1074.1 ms  ✓ Qhull_jll
    950.2 ms  ✓ Bzip2_jll
   2374.1 ms  ✓ SentinelArrays
   1982.7 ms  ✓ boost_jll
   1203.8 ms  ✓ Lz4_jll
   1127.7 ms  ✓ ICU_jll
    937.6 ms  ✓ libaec_jll
    889.9 ms  ✓ Xorg_libXau_jll
    880.2 ms  ✓ libpng_jll
    849.6 ms  ✓ Hwloc_jll
    842.7 ms  ✓ Giflib_jll
    869.8 ms  ✓ LERC_jll
    864.0 ms  ✓ SQLite_jll
    912.3 ms  ✓ EarCut_jll
    868.8 ms  ✓ JpegTurbo_jll
    835.9 ms  ✓ XZ_jll
    862.6 ms  ✓ snappy_jll
    709.2 ms  ✓ Xorg_libXdmcp_jll
    911.1 ms  ✓ brotli_jll
    952.3 ms  ✓ Zstd_jll
    831.5 ms  ✓ Expat_jll
    737.8 ms  ✓ MicrosoftMPI_jll
    823.2 ms  ✓ Xorg_xtrans_jll
    894.2 ms  ✓ GEOS_jll
    824.3 ms  ✓ Kerberos_krb5_jll
   1074.5 ms  ✓ Libgpg_error_jll
    807.5 ms  ✓ Xorg_libpthread_stubs_jll
    449.2 ms  ✓ TableTraits
   1103.4 ms  ✓ XML2_jll
   2469.1 ms  ✓ RecipesBase
   1021.9 ms  ✓ IntervalSets → IntervalSetsRandomExt
    488.6 ms  ✓ IntervalSets → IntervalSetsStatisticsExt
    875.0 ms  ✓ StackViews
    852.6 ms  ✓ PaddedViews
   2865.4 ms  ✓ DataStructures
   3584.6 ms  ✓ StringManipulation
   1023.0 ms  ✓ PooledArrays
    915.7 ms  ✓ Missings
   4598.0 ms  ✓ FixedPointNumbers
   1659.1 ms  ✓ GeoInterface
   1138.8 ms  ✓ Thrift_jll
   1928.5 ms  ✓ DiskArrays
   1568.0 ms  ✓ MPItrampoline_jll
   1541.1 ms  ✓ OpenMPI_jll
   1008.7 ms  ✓ Blosc_jll
   1276.1 ms  ✓ HDF4_jll
    964.5 ms  ✓ libzip_jll
    888.4 ms  ✓ Libtiff_jll
   1976.8 ms  ✓ MPICH_jll
    855.3 ms  ✓ Libgcrypt_jll
   1226.0 ms  ✓ LibPQ_jll
    744.3 ms  ✓ IntervalSets → IntervalSetsRecipesBaseExt
    651.4 ms  ✓ MosaicViews
   1179.7 ms  ✓ Tables
    809.5 ms  ✓ SortingAlgorithms
    990.0 ms  ✓ GeoInterfaceRecipes
   1221.9 ms  ✓ Arrow_jll
   1280.1 ms  ✓ PROJ_jll
   1132.0 ms  ✓ LittleCMS_jll
    875.4 ms  ✓ XSLT_jll
   1253.1 ms  ✓ libgeotiff_jll
   4107.9 ms  ✓ ColorTypes
   2876.7 ms  ✓ HDF5_jll
   1233.5 ms  ✓ OpenJpeg_jll
   2363.0 ms  ✓ Xorg_libxcb_jll
  12653.4 ms  ✓ StaticArrays
   2092.3 ms  ✓ NetCDF_jll
   1132.1 ms  ✓ Xorg_libX11_jll
   4226.6 ms  ✓ ColorVectorSpace
   1099.7 ms  ✓ StaticArrays → StaticArraysStatisticsExt
    964.5 ms  ✓ Xorg_libXext_jll
    923.6 ms  ✓ Libglvnd_jll
   1289.1 ms  ✓ libwebp_jll
   7888.5 ms  ✓ Colors
   2759.9 ms  ✓ GDAL_jll
   4967.9 ms  ✓ GDAL
  17594.7 ms  ✓ GeometryBasics
   6684.2 ms  ✓ MakieCore
  31966.4 ms  ✓ PrettyTables
   1856.4 ms  ✓ GeoInterfaceMakie
  26724.3 ms  ✓ ImageCore
  11180.8 ms  ✓ ArchGDAL
  42008.7 ms  ✓ DataFrames
   2573.5 ms  ✓ GeoDataFrames
  106 dependencies successfully precompiled in 96 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

using CairoMakie
using GeoMakie

f, a, p = poly(df.geom)

Precompiling CairoMakie...
    464.1 ms  ✓ RangeArrays
    476.9 ms  ✓ IndirectArrays
    508.0 ms  ✓ StatsAPI
    663.9 ms  ✓ PolygonOps
    786.2 ms  ✓ AbstractFFTs
    626.1 ms  ✓ Contour
    676.2 ms  ✓ TriplotBase
    558.6 ms  ✓ InverseFunctions
    704.5 ms  ✓ EnumX
   1591.9 ms  ✓ Format
   1721.9 ms  ✓ FillArrays
   1412.2 ms  ✓ Grisu
    619.1 ms  ✓ StableRNGs
   1001.9 ms  ✓ AbstractTrees
    506.5 ms  ✓ PtrArrays
    612.9 ms  ✓ RoundingEmulator
    833.1 ms  ✓ NaNMath
    923.3 ms  ✓ TranscodingStreams
   2139.8 ms  ✓ IrrationalConstants
    591.9 ms  ✓ Ratios
    716.4 ms  ✓ LazyModules
   3939.4 ms  ✓ MacroTools
    677.4 ms  ✓ Inflate
    643.4 ms  ✓ Adapt
    797.2 ms  ✓ ConstructionBase
    924.2 ms  ✓ Statistics → SparseArraysExt
   2761.8 ms  ✓ AdaptivePredicates
    882.8 ms  ✓ SuiteSparse
    999.1 ms  ✓ WoodburyMatrices
   1161.9 ms  ✓ ProgressMeter
    881.6 ms  ✓ DocStringExtensions
    891.1 ms  ✓ Animations
    660.0 ms  ✓ SignedDistanceFields
    883.6 ms  ✓ Graphite2_jll
    745.8 ms  ✓ Libmount_jll
   2126.3 ms  ✓ ChainRulesCore
   2180.0 ms  ✓ UnicodeFun
    968.4 ms  ✓ LLVMOpenMP_jll
    946.3 ms  ✓ Rmath_jll
    874.5 ms  ✓ libfdk_aac_jll
    756.3 ms  ✓ Imath_jll
    858.4 ms  ✓ LAME_jll
    842.7 ms  ✓ CRlibm_jll
    877.1 ms  ✓ Ogg_jll
    925.2 ms  ✓ oneTBB_jll
    827.3 ms  ✓ x265_jll
    872.0 ms  ✓ x264_jll
    884.7 ms  ✓ LZO_jll
    803.9 ms  ✓ Opus_jll
   1034.9 ms  ✓ libaom_jll
    733.2 ms  ✓ Libffi_jll
    794.6 ms  ✓ FFTW_jll
    809.7 ms  ✓ OpenSpecFun_jll
    953.5 ms  ✓ isoband_jll
    810.4 ms  ✓ Libuuid_jll
    764.4 ms  ✓ FriBidi_jll
    903.9 ms  ✓ OpenBLASConsistentFPCSR_jll
   1026.1 ms  ✓ libsixel_jll
   1634.9 ms  ✓ PkgVersion
   1849.8 ms  ✓ FilePathsBase
    967.1 ms  ✓ FreeType2_jll
   2538.3 ms  ✓ JSON
   5856.5 ms  ✓ ColorSchemes
   1967.6 ms  ✓ QuadGK
   7876.0 ms  ✓ FileIO
   4214.5 ms  ✓ ImageBase
   2086.9 ms  ✓ IntelOpenMP_jll
   1930.4 ms  ✓ Packing
   2835.9 ms  ✓ ShaderAbstractions
   1375.0 ms  ✓ Gettext_jll
   1106.0 ms  ✓ Xorg_libXrender_jll
   1549.7 ms  ✓ AxisArrays
  15892.3 ms  ✓ SIMD
    866.0 ms  ✓ InverseFunctions → InverseFunctionsTestExt
   2124.7 ms  ✓ AbstractFFTs → AbstractFFTsTestExt
    609.0 ms  ✓ InverseFunctions → InverseFunctionsDatesExt
    774.9 ms  ✓ FillArrays → FillArraysStatisticsExt
    915.1 ms  ✓ FillArrays → FillArraysSparseArraysExt
    738.4 ms  ✓ Showoff
    926.3 ms  ✓ AliasTables
    952.6 ms  ✓ Graphics
   1020.2 ms  ✓ Ratios → RatiosFixedPointNumbersExt
   1841.9 ms  ✓ SimpleTraits
   1341.1 ms  ✓ Adapt → AdaptStaticArraysExt
   8099.8 ms  ✓ PNGFiles
    883.8 ms  ✓ OffsetArrays → OffsetArraysAdaptExt
   1066.0 ms  ✓ Adapt → AdaptSparseArraysExt
    682.8 ms  ✓ ConstructionBase → ConstructionBaseLinearAlgebraExt
    736.2 ms  ✓ ConstructionBase → ConstructionBaseIntervalSetsExt
   1252.8 ms  ✓ ConstructionBase → ConstructionBaseStaticArraysExt
  14025.5 ms  ✓ GridLayoutBase
   1296.0 ms  ✓ PDMats
   1411.1 ms  ✓ AxisAlgorithms
   1398.5 ms  ✓ LogExpFunctions
    755.9 ms  ✓ AbstractFFTs → AbstractFFTsChainRulesCoreExt
   1170.6 ms  ✓ StaticArrays → StaticArraysChainRulesCoreExt
   1316.1 ms  ✓ ChainRulesCore → ChainRulesCoreSparseArraysExt
    976.7 ms  ✓ Pixman_jll
   1113.3 ms  ✓ libvorbis_jll
   1495.7 ms  ✓ Rmath
    712.0 ms  ✓ Isoband
   1376.9 ms  ✓ OpenEXR_jll
    977.1 ms  ✓ FilePathsBase → FilePathsBaseMmapExt
   1547.1 ms  ✓ FilePaths
   2396.7 ms  ✓ FilePathsBase → FilePathsBaseTestExt
   1419.3 ms  ✓ Fontconfig_jll
   2189.7 ms  ✓ FreeType
   1280.8 ms  ✓ ColorBrewer
   6682.4 ms  ✓ IntervalArithmetic
  41765.3 ms  ✓ Unitful
   5170.1 ms  ✓ Sixel
   2584.1 ms  ✓ QOI
   6272.1 ms  ✓ JpegTurbo
   1325.2 ms  ✓ Glib_jll
   3551.4 ms  ✓ WebP
   2913.1 ms  ✓ ImageAxes
   6603.2 ms  ✓ MKL_jll
   1285.2 ms  ✓ StructArrays
   1091.6 ms  ✓ FillArrays → FillArraysPDMatsExt
    743.6 ms  ✓ LogExpFunctions → LogExpFunctionsInverseFunctionsExt
   2580.6 ms  ✓ LogExpFunctions → LogExpFunctionsChainRulesCoreExt
   5302.4 ms  ✓ SpecialFunctions
   4535.9 ms  ✓ StatsBase
  22634.3 ms  ✓ PlotUtils
   2777.2 ms  ✓ OpenEXR
   4371.9 ms  ✓ Interpolations
   1214.1 ms  ✓ IntervalArithmetic → IntervalArithmeticIntervalSetsExt
   1239.9 ms  ✓ Unitful → ConstructionBaseUnitfulExt
   1078.0 ms  ✓ Unitful → InverseFunctionsUnitfulExt
   4676.4 ms  ✓ FreeTypeAbstraction
   1656.1 ms  ✓ Cairo_jll
   2437.3 ms  ✓ ImageMetadata
    781.7 ms  ✓ StructArrays → StructArraysAdaptExt
  24891.8 ms  ✓ Automa
   1397.0 ms  ✓ StructArrays → StructArraysSparseArraysExt
   1338.6 ms  ✓ StructArrays → StructArraysStaticArraysExt
    986.3 ms  ✓ StructArrays → StructArraysLinearAlgebraExt
   9693.1 ms  ✓ ExactPredicates
   1025.9 ms  ✓ ColorVectorSpace → SpecialFunctionsExt
   1818.1 ms  ✓ HypergeometricFunctions
   3060.8 ms  ✓ SpecialFunctions → SpecialFunctionsChainRulesCoreExt
    979.3 ms  ✓ HarfBuzz_jll
   1739.3 ms  ✓ Interpolations → InterpolationsUnitfulExt
   3316.5 ms  ✓ Netpbm
   3463.7 ms  ✓ StatsFuns
   1385.9 ms  ✓ libass_jll
   1314.4 ms  ✓ Pango_jll
   1063.6 ms  ✓ StatsFuns → StatsFunsInverseFunctionsExt
  10255.5 ms  ✓ DelaunayTriangulation
  12737.1 ms  ✓ MathTeXEngine
   2791.4 ms  ✓ StatsFuns → StatsFunsChainRulesCoreExt
  20894.2 ms  ✓ FFTW
   1389.7 ms  ✓ FFMPEG_jll
   1906.7 ms  ✓ Cairo
   6602.4 ms  ✓ Distributions
   1351.1 ms  ✓ Distributions → DistributionsChainRulesCoreExt
   1409.1 ms  ✓ Distributions → DistributionsTestExt
   1322.5 ms  ✓ KernelDensity
  82983.4 ms  ✓ TiffImages
    963.7 ms  ✓ ImageIO
 130031.5 ms  ✓ Makie
  40801.5 ms  ✓ CairoMakie
  162 dependencies successfully precompiled in 305 seconds. 110 already precompiled.
  1 dependency had output during precompilation:
┌ MKL_jll
│   Downloading artifact: IntelOpenMP
└  
Precompiling ParsersExt...
    413.1 ms  ✓ InlineStrings → ParsersExt
  1 dependency successfully precompiled in 0 seconds. 9 already precompiled.
Precompiling SerializationExt...
    335.3 ms  ✓ LRUCache → SerializationExt
  1 dependency successfully precompiled in 0 seconds. 2 already precompiled.
Precompiling IntervalArithmeticRecipesBaseExt...
    662.7 ms  ✓ IntervalArithmetic → IntervalArithmeticRecipesBaseExt
  1 dependency successfully precompiled in 1 seconds. 33 already precompiled.
Precompiling ArchGDALMakieExt...
   2101.3 ms  ✓ ArchGDAL → ArchGDALMakieExt
  1 dependency successfully precompiled in 3 seconds. 306 already precompiled.
Precompiling GeoMakie...
    528.5 ms  ✓ SortTileRecursiveTree
    939.5 ms  ✓ GeometryOpsCore
    991.9 ms  ✓ CoordinateTransformations
   2357.9 ms  ✓ Geodesy
   2585.6 ms  ✓ Proj
   4088.6 ms  ✓ GeometryOps
   1186.5 ms  ✓ GeometryOps → GeometryOpsProjExt
  24586.5 ms  ✓ GeoJSON
   1074.2 ms  ✓ NaturalEarth
   1321.4 ms  ✓ GeoJSON → GeoJSONMakieExt
   8679.3 ms  ✓ GeoMakie
  11 dependencies successfully precompiled in 35 seconds. 276 already precompiled.

1.1.1 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.

using Rasters
import GeoFormatTypes as GFT

Precompiling Rasters...
    668.8 ms  ✓ Interfaces
    673.1 ms  ✓ FieldMetadata
    830.2 ms  ✓ ArrayInterface
   1153.4 ms  ✓ CFTime
    803.5 ms  ✓ Flatten
    638.4 ms  ✓ ArrayInterface → ArrayInterfaceStaticArraysCoreExt
   1562.1 ms  ✓ Setfield
    848.8 ms  ✓ ArrayInterface → ArrayInterfaceSparseArraysExt
   1397.5 ms  ✓ CommonDataModel
  14260.9 ms  ✓ DimensionalData
    754.7 ms  ✓ DimensionalData → DimensionalDataDiskArraysExt
   5360.0 ms  ✓ Rasters
  12 dependencies successfully precompiled in 22 seconds. 58 already precompiled.
Precompiling ArrayInterfaceChainRulesCoreExt...
    364.4 ms  ✓ ArrayInterface → ArrayInterfaceChainRulesCoreExt
  1 dependency successfully precompiled in 0 seconds. 11 already precompiled.
Precompiling DimensionalDataMakie...
    806.6 ms  ✓ DimensionalData → DimensionalDataStatsBase
   6423.3 ms  ✓ DimensionalData → DimensionalDataMakie
  2 dependencies successfully precompiled in 7 seconds. 277 already precompiled.
Precompiling RastersMakieExt...
   1292.0 ms  ✓ Rasters → RastersStatsBaseExt
   6214.7 ms  ✓ Rasters → RastersMakieExt
  2 dependencies successfully precompiled in 7 seconds. 291 already precompiled.
Precompiling RastersProjExt...
   1225.8 ms  ✓ Rasters → RastersCoordinateTransformationsExt
   1340.7 ms  ✓ Rasters → RastersProjExt
  2 dependencies successfully precompiled in 1 seconds. 92 already precompiled.
Precompiling RastersArchGDALExt...
   2226.8 ms  ✓ Rasters → RastersArchGDALExt
  1 dependency successfully precompiled in 2 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.

elev = reshape(UInt8(1):UInt8(36), (6, 6))

6×6 reshape(::UnitRange{UInt8}, 6, 6) with eltype UInt8:
 0x01  0x07  0x0d  0x13  0x19  0x1f
 0x02  0x08  0x0e  0x14  0x1a  0x20
 0x03  0x09  0x0f  0x15  0x1b  0x21
 0x04  0x0a  0x10  0x16  0x1c  0x22
 0x05  0x0b  0x11  0x17  0x1d  0x23
 0x06  0x0c  0x12  0x18  0x1e  0x24

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.

v = UInt8[
  1, 0, 1, 2, 2, 2, 
  0, 2, 0, 0, 2, 1, 
  0, 2, 2, 0, 0, 2, 
  0, 0, 1, 1, 1, 1, 
  1, 1, 1, 2, 1, 1, 
  2, 1, 2, 2, 0, 2
]
grain = reshape(v, (6, 6))

6×6 Matrix{UInt8}:
 0x01  0x00  0x00  0x00  0x01  0x02
 0x00  0x02  0x02  0x00  0x01  0x01
 0x01  0x00  0x02  0x01  0x01  0x02
 0x02  0x00  0x00  0x01  0x02  0x02
 0x02  0x02  0x00  0x01  0x01  0x00
 0x02  0x01  0x02  0x01  0x01  0x02

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.

new_x = X(range(-1.5, step=0.5, length=6))
new_y = Y(range(1.0, step=-0.5, length=6))

Y 1.0:-0.5:-1.5

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.

elev_raster = Raster(elev, (new_x, new_y); crs = GFT.EPSG(4326))

┌ 6×6 Raster{UInt8, 2} ┐
├──────────────────────┴───────────────────────────────────────────────── dims ┐
  ↓ X Projected{Float64} -1.5:0.5:1.0 ForwardOrdered Regular Points,
  → Y Projected{Float64} 1.0:-0.5:-1.5 ReverseOrdered Regular Points
├────────────────────────────────────────────────────────────────────── raster ┤
  extent: Extent(X = (-1.5, 1.0), Y = (-1.5, 1.0))
  crs: EPSG:4326
└──────────────────────────────────────────────────────────────────────────────┘
  ↓ →     1.0     0.5     0.0    -0.5    -1.0    -1.5
 -1.5  0x01    0x07    0x0d    0x13    0x19    0x1f
 -1.0  0x02    0x08    0x0e    0x14    0x1a    0x20
 -0.5  0x03    0x09    0x0f    0x15    0x1b    0x21
  0.0  0x04    0x0a    0x10    0x16    0x1c    0x22
  0.5  0x05    0x0b    0x11    0x17    0x1d    0x23
  1.0  0x06    0x0c    0x12    0x18    0x1e    0x24

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.1).

plot(elev_raster)

Figure 1.1: 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.2).

plot(Raster(grain, (new_x, new_y); crs = GFT.EPSG(4326)))

Figure 1.2: 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:

1. Create a raster file (where we set the lookups and the CRS, among other settings)
2. Write the array with raster values into the connection
3. 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.