Quick start
Quick start
Install the package by typing:
]
add Rasters
then do
using Rasters
Using Rasters to read GeoTiff or NetCDF files will output something similar to the following toy examples. This is possible because Rasters.jl extends DimensionalData.jl so that spatial data can be indexed using named dimensions like X, Y and Ti (time) and e.g. spatial coordinates.
using Rasters, Dates
lon, lat = X(25:1:30), Y(25:1:30)
ti = Ti(DateTime(2001):Month(1):DateTime(2002))
ras = Raster(rand(lon, lat, ti)) # this generates random numbers with the dimensions given
╭──────────────────────────╮
│ 6×6×13 Raster{Float64,3} │
├──────────────────────────┴───────────────────────────────────────────── dims ┐
↓ X Sampled{Int64} 25:1:30 ForwardOrdered Regular Points,
→ Y Sampled{Int64} 25:1:30 ForwardOrdered Regular Points,
↗ Ti Sampled{DateTime} DateTime("2001-01-01T00:00:00"):Month(1):DateTime("2002-01-01T00:00:00") ForwardOrdered Regular Points
├────────────────────────────────────────────────────────────────────── raster ┤
extent: Extent(X = (25, 30), Y = (25, 30), Ti = (DateTime("2001-01-01T00:00:00"), DateTime("2002-01-01T00:00:00")))
└──────────────────────────────────────────────────────────────────────────────┘
[:, :, 1]
↓ → 25 26 27 28 29 30
25 0.755422 0.799349 0.186365 0.851987 0.445698 0.582066
26 0.911574 0.591608 0.185338 0.132995 0.239419 0.752367
27 0.0831978 0.310477 0.38834 0.772493 0.147646 0.867995
28 0.0432791 0.349809 0.480449 0.907303 0.12202 0.702061
29 0.677934 0.781828 0.584605 0.13499 0.0296433 0.951586
30 0.823539 0.18972 0.771593 0.22449 0.257423 0.375458
Getting the lookup array from dimensions¤
lon = lookup(ras, X) # if X is longitude
lat = lookup(ras, Y) # if Y is latitude
Sampled{Int64} ForwardOrdered Regular Points
wrapping: 25:1:30
Select by index¤
Selecting a time slice by index is done via
ras[Ti(1)]
╭───────────────────────╮
│ 6×6 Raster{Float64,2} │
├───────────────────────┴───────────────────────────── dims ┐
↓ X Sampled{Int64} 25:1:30 ForwardOrdered Regular Points,
→ Y Sampled{Int64} 25:1:30 ForwardOrdered Regular Points
├─────────────────────────────────────────────────── raster ┤
extent: Extent(X = (25, 30), Y = (25, 30))
└───────────────────────────────────────────────────────────┘
↓ → 25 26 27 28 29 30
25 0.755422 0.799349 0.186365 0.851987 0.445698 0.582066
26 0.911574 0.591608 0.185338 0.132995 0.239419 0.752367
27 0.0831978 0.310477 0.38834 0.772493 0.147646 0.867995
28 0.0432791 0.349809 0.480449 0.907303 0.12202 0.702061
29 0.677934 0.781828 0.584605 0.13499 0.0296433 0.951586
30 0.823539 0.18972 0.771593 0.22449 0.257423 0.375458
also
ras[Ti=1]
╭───────────────────────╮
│ 6×6 Raster{Float64,2} │
├───────────────────────┴───────────────────────────── dims ┐
↓ X Sampled{Int64} 25:1:30 ForwardOrdered Regular Points,
→ Y Sampled{Int64} 25:1:30 ForwardOrdered Regular Points
├─────────────────────────────────────────────────── raster ┤
extent: Extent(X = (25, 30), Y = (25, 30))
└───────────────────────────────────────────────────────────┘
↓ → 25 26 27 28 29 30
25 0.755422 0.799349 0.186365 0.851987 0.445698 0.582066
26 0.911574 0.591608 0.185338 0.132995 0.239419 0.752367
27 0.0831978 0.310477 0.38834 0.772493 0.147646 0.867995
28 0.0432791 0.349809 0.480449 0.907303 0.12202 0.702061
29 0.677934 0.781828 0.584605 0.13499 0.0296433 0.951586
30 0.823539 0.18972 0.771593 0.22449 0.257423 0.375458
or and interval of indices using the syntax =a:b or (a:b)
ras[Ti(1:10)]
╭──────────────────────────╮
│ 6×6×10 Raster{Float64,3} │
├──────────────────────────┴───────────────────────────────────────────── dims ┐
↓ X Sampled{Int64} 25:1:30 ForwardOrdered Regular Points,
→ Y Sampled{Int64} 25:1:30 ForwardOrdered Regular Points,
↗ Ti Sampled{DateTime} DateTime("2001-01-01T00:00:00"):Month(1):DateTime("2001-10-01T00:00:00") ForwardOrdered Regular Points
├────────────────────────────────────────────────────────────────────── raster ┤
extent: Extent(X = (25, 30), Y = (25, 30), Ti = (DateTime("2001-01-01T00:00:00"), DateTime("2001-10-01T00:00:00")))
└──────────────────────────────────────────────────────────────────────────────┘
[:, :, 1]
↓ → 25 26 27 28 29 30
25 0.755422 0.799349 0.186365 0.851987 0.445698 0.582066
26 0.911574 0.591608 0.185338 0.132995 0.239419 0.752367
27 0.0831978 0.310477 0.38834 0.772493 0.147646 0.867995
28 0.0432791 0.349809 0.480449 0.907303 0.12202 0.702061
29 0.677934 0.781828 0.584605 0.13499 0.0296433 0.951586
30 0.823539 0.18972 0.771593 0.22449 0.257423 0.375458
Select by value¤
ras[Ti=At(DateTime(2001))]
╭───────────────────────╮
│ 6×6 Raster{Float64,2} │
├───────────────────────┴───────────────────────────── dims ┐
↓ X Sampled{Int64} 25:1:30 ForwardOrdered Regular Points,
→ Y Sampled{Int64} 25:1:30 ForwardOrdered Regular Points
├─────────────────────────────────────────────────── raster ┤
extent: Extent(X = (25, 30), Y = (25, 30))
└───────────────────────────────────────────────────────────┘
↓ → 25 26 27 28 29 30
25 0.755422 0.799349 0.186365 0.851987 0.445698 0.582066
26 0.911574 0.591608 0.185338 0.132995 0.239419 0.752367
27 0.0831978 0.310477 0.38834 0.772493 0.147646 0.867995
28 0.0432791 0.349809 0.480449 0.907303 0.12202 0.702061
29 0.677934 0.781828 0.584605 0.13499 0.0296433 0.951586
30 0.823539 0.18972 0.771593 0.22449 0.257423 0.375458
More options are available, like Near, Contains and Where.
Dimensions
Rasters uses X, Y, and Z dimensions from DimensionalData to represent spatial directions like longitude, latitude and the vertical dimension, and subset data with them. Ti is used for time, and Band represent bands. Other dimensions can have arbitrary names, but will be treated generically. See DimensionalData for more details on how they work.
Lookup Arrays
These specify properties of the index associated with e.g. the X and Y dimension. Rasters.jl defines additional lookup arrays: Projected to handle dimensions with projections, and Mapped where the projection is mapped to another projection like EPSG(4326). Mapped is largely designed to handle NetCDF dimensions, especially with Explicit spans.
Subsetting an object¤
Regular getindex (e.g. A[1:100, :]) and view work on all objects just as with an Array. view is always lazy, and reads from disk are deferred until getindex is used. DimensionalData.jl Dimensions and Selectors are the other way to subset an object, making use of the objects index to find values at e.g. certain X/Y coordinates. The available selectors are listed here:
Selectors |
Description |
|---|---|
At(x) |
get the index exactly matching the passed in value(s). |
Near(x) |
get the closest index to the passed in value(s). |
Where(f::Function) |
filter the array axis by a function of the dimension index values. |
a..b/Between(a, b) |
get all indices between two values, excluding the high value. |
Contains(x) |
get indices where the value x falls within an interval. |
Info
- Use the
..selector to take a view of madagascar:
using Rasters, RasterDataSources
const RS = Rasters
using CairoMakie
CairoMakie.activate!()
A = Raster(WorldClim{BioClim}, 5)
madagascar = view(A, X(43.25 .. 50.48), Y(-25.61 .. -12.04))
╭──────────────────────────────╮
│ 42×80 Raster{Float32,2} bio5 │
├──────────────────────────────┴───────────────────────────────────────── dims ┐
↓ X Projected{Float64} LinRange{Float64}(43.33333333333334, 50.16666666666667, 42) ForwardOrdered Regular Intervals{Start},
→ Y Projected{Float64} LinRange{Float64}(-12.333333333333343, -25.5, 80) ReverseOrdered Regular Intervals{Start}
├──────────────────────────────────────────────────────────────────── metadata ┤
Metadata{Rasters.GDALsource} of Dict{String, Any} with 4 entries:
"units" => ""
"offset" => 0.0
"filepath" => "/tmp/WorldClim/BioClim/wc2.1_10m_bio_5.tif"
"scale" => 1.0
├────────────────────────────────────────────────────────────────────── raster ┤
extent: Extent(X = (43.33333333333334, 50.333333333333336), Y = (-25.5, -12.166666666666677))
missingval: -3.4f38
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"]]
└──────────────────────────────────────────────────────────────────────────────┘
↓ → -12.3333 -12.5 -12.6667 … -25.1667 -25.3333 -25.5
43.3333 -3.4f38 -3.4f38 -3.4f38 -3.4f38 -3.4f38 -3.4f38
⋮ ⋱ ⋮
50.1667 -3.4f38 -3.4f38 -3.4f38 -3.4f38 -3.4f38 -3.4f38
Note the space between .. -12
Makie.plot(madagascar)
This page was generated using Literate.jl.