Examples and Plotting

Plots, simple

Plots.jl and Makie.jl are fully supported by Rasters.jl, with recipes for plotting Raster and RasterStack provided. plot will plot a heatmap with axes matching dimension values. If mappedcrs is used, converted values will be shown on axes instead of the underlying crs values. contourf will similarly plot a filled contour plot.

Pixel resolution is limited to allow loading very large files quickly. max_res specifies the maximum pixel resolution to show on the longest axis of the array. It can be set manually to change the resolution (e.g. for large or high-quality plots):

using Rasters, RasterDataSources, ArchGDAL, Plots
A = Raster(WorldClim{BioClim}, 5)
plot(A; max_res=3000)

For Makie, plot functions in a similar way. plot will only accept two-dimensional rasters. You can invoke contour, contourf, heatmap, surface or any Makie plotting function which supports surface-like data on a 2D raster.

To obtain tiled plots for 3D rasters and RasterStacks, use the function Rasters.rplot([gridposition], raster; kw_args...). This is an unexported function, since we're not sure how the API will change going forward.

Makie, simple

using CairoMakie
CairoMakie.activate!(px_per_unit = 2)
using Rasters, CairoMakie, RasterDataSources, ArchGDAL
A = Raster(WorldClim{BioClim}, 5)
Makie.plot(A)

Loading data

Our first example simply loads a file from disk and plots it.

This netcdf file only has one layer, if it has more we could use RasterStack instead.

using Rasters, NCDatasets, Plots
url = "https://www.unidata.ucar.edu/software/netcdf/examples/tos_O1_2001-2002.nc";
filename = download(url, "tos_O1_2001-2002.nc");
A = Raster(filename)

╭──────────────────────────────────────────────────╮
│ 180×170×24 Raster{Union{Missing, Float32},3} tos │
├──────────────────────────────────────────────────┴───────────────────── dims ┐
  ↓ X  Mapped{Float64} [1.0, 3.0, …, 357.0, 359.0] ForwardOrdered Regular Intervals{Center},
  → Y  Mapped{Float64} [-79.5, -78.5, …, 88.5, 89.5] ForwardOrdered Regular Intervals{Center},
  ↗ Ti Sampled{DateTime360Day} [DateTime360Day(2001-01-16T00:00:00), …, DateTime360Day(2002-12-16T00:00:00)] ForwardOrdered Explicit Intervals{Center}
├──────────────────────────────────────────────────────────────────── metadata ┤
  Metadata{Rasters.NCDsource} of Dict{String, Any} with 9 entries:
  "units"          => "K"
  "missing_value"  => 1.0f20
  "original_units" => "degC"
  "cell_methods"   => "time: mean (interval: 30 minutes)"
  "history"        => " At   16:37:23 on 01/11/2005: CMOR altered the data in t…
  "long_name"      => "Sea Surface Temperature"
  "standard_name"  => "sea_surface_temperature"
  "_FillValue"     => 1.0f20
  "original_name"  => "sosstsst"
├────────────────────────────────────────────────────────────────────── raster ┤
  extent: Extent(X = (0.0, 360.0), Y = (-80.0, 90.0), Ti = (DateTime360Day(2001-01-01T00:00:00), DateTime360Day(2003-01-01T00:00:00)))
  missingval: missing
  crs: EPSG:4326
  mappedcrs: EPSG:4326
└──────────────────────────────────────────────────────────────────────────────┘
[:, :, 1]
 ⋮      ⋱  

Objects with Dimensions other than X and Y will produce multi-pane plots. Here we plot every third month in the first year in one plot:

A[Ti=1:3:12] |> plot

Now plot the ocean temperatures around the Americas in the first month of 2001. Notice we are using lat/lon coordinates and date/time instead of regular indices. The time dimension uses DateTime360Day, so we need to load CFTime.jl to index it with Near.

using CFTime
A[Ti(Near(DateTime360Day(2001, 01, 17))), Y(-60.0 .. 90.0), X(45.0 .. 190.0)] |> plot

Now get the mean over the timespan, then save it to disk, and plot it as a filled contour.

Other plot functions and sliced objects that have only one X/Y/Z dimension fall back to generic DimensionalData.jl plotting, which will still correctly label plot axes.

using Statistics
# Take the mean
mean_tos = mean(A; dims=Ti)

╭─────────────────────────────────────────────────╮
│ 180×170×1 Raster{Union{Missing, Float32},3} tos │
├─────────────────────────────────────────────────┴────────────────────── dims ┐
  ↓ X  Mapped{Float64} [1.0, 3.0, …, 357.0, 359.0] ForwardOrdered Regular Intervals{Center},
  → Y  Mapped{Float64} [-79.5, -78.5, …, 88.5, 89.5] ForwardOrdered Regular Intervals{Center},
  ↗ Ti Sampled{DateTime360Day} DateTime360Day(2002-01-16T00:00:00):Millisecond(2592000000):DateTime360Day(2002-01-16T00:00:00) ForwardOrdered Explicit Intervals{Center}
├──────────────────────────────────────────────────────────────────── metadata ┤
  Metadata{Rasters.NCDsource} of Dict{String, Any} with 9 entries:
  "units"          => "K"
  "missing_value"  => 1.0f20
  "original_units" => "degC"
  "cell_methods"   => "time: mean (interval: 30 minutes)"
  "history"        => " At   16:37:23 on 01/11/2005: CMOR altered the data in t…
  "long_name"      => "Sea Surface Temperature"
  "standard_name"  => "sea_surface_temperature"
  "_FillValue"     => 1.0f20
  "original_name"  => "sosstsst"
├────────────────────────────────────────────────────────────────────── raster ┤
  extent: Extent(X = (0.0, 360.0), Y = (-80.0, 90.0), Ti = (DateTime360Day(2001-01-01T00:00:00), DateTime360Day(2003-01-01T00:00:00)))
  missingval: missing
  crs: EPSG:4326
  mappedcrs: EPSG:4326
└──────────────────────────────────────────────────────────────────────────────┘
[:, :, 1]
 ⋮      ⋱  

Plot a contour plot

using Plots
Plots.contourf(mean_tos; dpi=300, size=(800, 400))

write to disk

Write the mean values to disk

write("mean_tos.nc", mean_tos)

"mean_tos.nc"

Plotting recipes in DimensionalData.jl are the fallback for Rasters.jl when the object doesn't have 2 X/Y/Z dimensions, or a non-spatial plot command is used. So (as a random example) we could plot a transect of ocean surface temperature at 20 degree latitude :

A[Y(Near(20.0)), Ti(1)] |> plot

Polygon masking, mosaic and plot

In this example we will mask the Scandinavian countries with border polygons, then mosaic together to make a single plot.

First, get the country boundary shape files using GADM.jl.

using Rasters, RasterDataSources, ArchGDAL, Shapefile, Plots, Dates, Downloads, NCDatasets

WARNING: using Downloads.download in module Main conflicts with an existing identifier.

Download the shapefile

shapefile_url = "https://github.com/nvkelso/natural-earth-vector/raw/master/10m_cultural/ne_10m_admin_0_countries.shp"
shapefile_name = "boundary_lines.shp"
Downloads.download(shapefile_url, shapefile_name)

"boundary_lines.shp"

Load using Shapefile.jl

shapes = Shapefile.Handle(shapefile_name)
denmark_border = shapes.shapes[71]
norway_border = shapes.shapes[53]
sweden_border = shapes.shapes[54]

Polygon(4665 Points)

Then load raster data. We load some worldclim layers using RasterDataSources via Rasters.jl:

climate = RasterStack(WorldClim{Climate}, (:tmin, :tmax, :prec, :wind); month=July)

╭───────────────────────╮
│ 2160×1080 RasterStack │
├───────────────────────┴──────────────────────────────────────────────── dims ┐
  ↓ X Projected{Float64} LinRange{Float64}(-180.0, 179.83333333333331, 2160) ForwardOrdered Regular Intervals{Start},
  → Y Projected{Float64} LinRange{Float64}(89.83333333333333, -90.0, 1080) ReverseOrdered Regular Intervals{Start}
├────────────────────────────────────────────────────────────────────── layers ┤
  :tmin eltype: Float32 dims: X, Y size: 2160×1080
  :tmax eltype: Float32 dims: X, Y size: 2160×1080
  :prec eltype: Int16 dims: X, Y size: 2160×1080
  :wind eltype: Float32 dims: X, Y size: 2160×1080
├────────────────────────────────────────────────────────────────────── raster ┤
  extent: Extent(X = (-180.0, 179.99999999999997), Y = (-90.0, 90.0))
  missingval: (tmin = -3.4f38, tmax = -3.4f38, prec = -32768, wind = -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"]]
└──────────────────────────────────────────────────────────────────────────────┘

mask Denmark, Norway and Sweden from the global dataset using their border polygon, then trim the missing values. We pad trim with a 10 pixel margin.

mask_trim(climate, poly) = trim(mask(climate; with=poly); pad=10)

denmark = mask_trim(climate, denmark_border)
norway = mask_trim(climate, norway_border)
sweden = mask_trim(climate, sweden_border)

╭────────────────────╮
│ 98×102 RasterStack │
├────────────────────┴─────────────────────────────────────────────────── dims ┐
  ↓ X Projected{Float64} LinRange{Float64}(9.499999999999972, 25.666666666666643, 98) ForwardOrdered Regular Intervals{Start},
  → Y Projected{Float64} LinRange{Float64}(70.5, 53.66666666666667, 102) ReverseOrdered Regular Intervals{Start}
├────────────────────────────────────────────────────────────────────── layers ┤
  :tmin eltype: Float32 dims: X, Y size: 98×102
  :tmax eltype: Float32 dims: X, Y size: 98×102
  :prec eltype: Int16 dims: X, Y size: 98×102
  :wind eltype: Float32 dims: X, Y size: 98×102
├────────────────────────────────────────────────────────────────────── raster ┤
  extent: Extent(X = (9.499999999999972, 25.83333333333331), Y = (53.66666666666667, 70.66666666666667))
  missingval: (tmin = -3.4f38, tmax = -3.4f38, prec = -32768, wind = -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"]]
└──────────────────────────────────────────────────────────────────────────────┘

Then load raster data. We load some worldclim layers using RasterDataSources via Rasters.jl:

climate = RasterStack(WorldClim{Climate}, (:tmin, :tmax, :prec, :wind); month=July)

╭───────────────────────╮
│ 2160×1080 RasterStack │
├───────────────────────┴──────────────────────────────────────────────── dims ┐
  ↓ X Projected{Float64} LinRange{Float64}(-180.0, 179.83333333333331, 2160) ForwardOrdered Regular Intervals{Start},
  → Y Projected{Float64} LinRange{Float64}(89.83333333333333, -90.0, 1080) ReverseOrdered Regular Intervals{Start}
├────────────────────────────────────────────────────────────────────── layers ┤
  :tmin eltype: Float32 dims: X, Y size: 2160×1080
  :tmax eltype: Float32 dims: X, Y size: 2160×1080
  :prec eltype: Int16 dims: X, Y size: 2160×1080
  :wind eltype: Float32 dims: X, Y size: 2160×1080
├────────────────────────────────────────────────────────────────────── raster ┤
  extent: Extent(X = (-180.0, 179.99999999999997), Y = (-90.0, 90.0))
  missingval: (tmin = -3.4f38, tmax = -3.4f38, prec = -32768, wind = -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"]]
└──────────────────────────────────────────────────────────────────────────────┘

mask Denmark, Norway and Sweden from the global dataset using their border polygon, then trim the missing values. We pad trim with a 10 pixel margin.

mask_trim(climate, poly) = trim(mask(climate; with=poly); pad=10)

denmark = mask_trim(climate, denmark_border)
norway = mask_trim(climate, norway_border)
sweden = mask_trim(climate, sweden_border)

╭────────────────────╮
│ 98×102 RasterStack │
├────────────────────┴─────────────────────────────────────────────────── dims ┐
  ↓ X Projected{Float64} LinRange{Float64}(9.499999999999972, 25.666666666666643, 98) ForwardOrdered Regular Intervals{Start},
  → Y Projected{Float64} LinRange{Float64}(70.5, 53.66666666666667, 102) ReverseOrdered Regular Intervals{Start}
├────────────────────────────────────────────────────────────────────── layers ┤
  :tmin eltype: Float32 dims: X, Y size: 98×102
  :tmax eltype: Float32 dims: X, Y size: 98×102
  :prec eltype: Int16 dims: X, Y size: 98×102
  :wind eltype: Float32 dims: X, Y size: 98×102
├────────────────────────────────────────────────────────────────────── raster ┤
  extent: Extent(X = (9.499999999999972, 25.83333333333331), Y = (53.66666666666667, 70.66666666666667))
  missingval: (tmin = -3.4f38, tmax = -3.4f38, prec = -32768, wind = -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"]]
└──────────────────────────────────────────────────────────────────────────────┘

Plotting with Plots

First define a function to add borders to all subplots.

function borders!(p, poly)
    for i in 1:length(p)
        Plots.plot!(p, poly; subplot=i, fillalpha=0, linewidth=0.6)
    end
    return p
end

borders! (generic function with 1 method)

Now we can plot the individual countries.

dp = plot(denmark)
borders!(dp, denmark_border)

and sweden

sp = plot(sweden)
borders!(sp, sweden_border)

and norway

np = plot(norway)
borders!(np, norway_border)

The Norway shape includes a lot of islands. Lets crop them out using .. intervals:

norway_region = climate[X(0..40), Y(55..73)]
plot(norway_region)

And mask it with the border again:

norway = mask_trim(norway_region, norway_border)
np = plot(norway)
borders!(np, norway_border)

Now we can combine the countries into a single raster using mosaic. first will take the first value if/when there is an overlap.

scandinavia = mosaic(first, denmark, norway, sweden)

╭─────────────────────╮
│ 177×119 RasterStack │
├─────────────────────┴────────────────────────────────────────────────── dims ┐
  ↓ X Projected{Float64} 3.1666666666666443:0.16666666666666666:32.49999999999998 ForwardOrdered Regular Intervals{Start},
  → Y Projected{Float64} 72.66666666666666:-0.16666666666666666:52.99999999999999 ReverseOrdered Regular Intervals{Start}
├────────────────────────────────────────────────────────────────────── layers ┤
  :tmin eltype: Float32 dims: X, Y size: 177×119
  :tmax eltype: Float32 dims: X, Y size: 177×119
  :prec eltype: Int16 dims: X, Y size: 177×119
  :wind eltype: Float32 dims: X, Y size: 177×119
├────────────────────────────────────────────────────────────────────── raster ┤
  extent: Extent(X = (3.1666666666666443, 32.66666666666664), Y = (52.99999999999999, 72.83333333333333))
  missingval: (tmin = -3.4f38, tmax = -3.4f38, prec = -32768, wind = -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"]]
└──────────────────────────────────────────────────────────────────────────────┘

And plot scandinavia, with all borders included:

p = plot(scandinavia)
borders!(p, denmark_border)
borders!(p, norway_border)
borders!(p, sweden_border)
p

And save to netcdf - a single multi-layered file, and tif, which will write a file for each stack layer.

write("scandinavia.nc", scandinavia)
write("scandinavia.tif", scandinavia)

(tmin = "scandinavia_tmin.tif", tmax = "scandinavia_tmax.tif", prec = "scandinavia_prec.tif", wind = "scandinavia_wind.tif")

Rasters.jl provides a range of other methods that are being added to over time. Where applicable these methods read and write lazily to and from disk-based arrays of common raster file types. These methods also work for entire RasterStacks and RasterSeries using the same syntax.

---

This page was generated using Literate.jl.