API Reference

Index

• Rasters.Rasters
• Rasters.AbstractProjected
• Rasters.AbstractRaster
• Rasters.AbstractRasterSeries
• Rasters.AbstractRasterStack
• Rasters.Band
• Rasters.FileArray
• Rasters.FileStack
• Rasters.Mapped
• Rasters.OpenStack
• Rasters.Projected
• Rasters.Raster
• Rasters.RasterDiskArray
• Rasters.RasterSeries
• Rasters.RasterStack
• Rasters.aggregate
• Rasters.aggregate!
• Rasters.boolmask
• Rasters.cellsize
• Rasters.classify
• Rasters.classify!
• Rasters.combine
• Rasters.convertlookup
• Rasters.coverage
• Rasters.coverage!
• Rasters.crop
• Rasters.disaggregate
• Rasters.disaggregate!
• Rasters.extend
• Rasters.extract
• Rasters.mappedbounds
• Rasters.mappedcrs
• Rasters.mappedindex
• Rasters.mask
• Rasters.mask!
• Rasters.missingmask
• Rasters.missingval
• Rasters.mosaic
• Rasters.mosaic!
• Rasters.points
• Rasters.rasterize
• Rasters.rasterize!
• Rasters.replace_missing
• Rasters.reproject
• Rasters.reproject
• Rasters.resample
• Rasters.rplot
• Rasters.rplot
• Rasters.setcrs
• Rasters.setmappedcrs
• Rasters.slice
• Rasters.trim
• Rasters.warp
• Rasters.zonal

Reference - Exported functions

# Rasters.Rasters — Module.

Rasters

Rasters.jl defines common types and methods for reading, writing and manipulating rasterized spatial data.

These currently include raster arrays like GeoTIFF and NetCDF, R grd files, multi-layered stacks, and multi-file series of arrays and stacks.

:warning: Packages extensions and Rasters 0.8 and onwards

On Julia 1.9 we can put additional packages in extensions, so the code only loads when you load a specific package. Rasters.jl was always intended to work like this, and its finally possible. This reduced package using time from many seconds to well under a second.

But, it means you have to manually load packages you need for each backend or additional functionality.

For example, to use the GDAL backend, and download files, you now need to do:

using Rasters, ArchGDAL, RasterDataSources

where previously it was just using Rasters.

Sources and packages needed:

• :gdal: using ArchGDAL
• :netcdf: using NCDatasets
• :grd: built-in.
• :smap: using HDF5
• :grib: using GRIBDatasets.

Other functionality in extensions:

• Raster data downloads, like Worldclim{Climate}: using RasterDataSources
• Makie plots: using Makie
• Coordinate transformations for gdal rasters: using CoordinateTransformations

Quick start

Install the package by typing:

]
add Rasters

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} with dimensions: 
  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
extent: Extent(X = (25, 30), Y = (25, 30), Ti = (DateTime("2001-01-01T00:00:00"), DateTime("2002-01-01T00:00:00")))
missingval: missing
values: [:, :, 1]
     25         26          27          28         29          30
 25   0.9063     0.427328    0.0320967   0.297023   0.0571002   0.891377
 26   0.443494   0.867547    0.350546    0.150155   0.24565     0.711039
 27   0.745673   0.0991336   0.930332    0.893537   0.805931    0.360583
 28   0.512083   0.125287    0.959434    0.354868   0.337824    0.259563
 29   0.253849   0.692209    0.774092    0.131798   0.823656    0.390013
 30   0.334152   0.136551    0.183555    0.941133   0.450484    0.461862
[and 12 more slices...]

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} with dimensions: 
  X Sampled{Int64} 25:1:30 ForwardOrdered Regular Points,
  Y Sampled{Int64} 25:1:30 ForwardOrdered Regular Points
and reference dimensions: 
  Ti Sampled{DateTime} DateTime("2001-01-01T00:00:00"):Month(1):DateTime("2001-01-01T00:00:00") ForwardOrdered Regular Points
extent: Extent(X = (25, 30), Y = (25, 30))
missingval: missing
values:      25         26          27          28         29          30
 25   0.9063     0.427328    0.0320967   0.297023   0.0571002   0.891377
 26   0.443494   0.867547    0.350546    0.150155   0.24565     0.711039
 27   0.745673   0.0991336   0.930332    0.893537   0.805931    0.360583
 28   0.512083   0.125287    0.959434    0.354868   0.337824    0.259563
 29   0.253849   0.692209    0.774092    0.131798   0.823656    0.390013
 30   0.334152   0.136551    0.183555    0.941133   0.450484    0.461862

ras[Ti=1]

6×6 Raster{Float64,2} with dimensions: 
  X Sampled{Int64} 25:1:30 ForwardOrdered Regular Points,
  Y Sampled{Int64} 25:1:30 ForwardOrdered Regular Points
and reference dimensions: 
  Ti Sampled{DateTime} DateTime("2001-01-01T00:00:00"):Month(1):DateTime("2001-01-01T00:00:00") ForwardOrdered Regular Points
extent: Extent(X = (25, 30), Y = (25, 30))
missingval: missing
values:      25         26          27          28         29          30
 25   0.9063     0.427328    0.0320967   0.297023   0.0571002   0.891377
 26   0.443494   0.867547    0.350546    0.150155   0.24565     0.711039
 27   0.745673   0.0991336   0.930332    0.893537   0.805931    0.360583
 28   0.512083   0.125287    0.959434    0.354868   0.337824    0.259563
 29   0.253849   0.692209    0.774092    0.131798   0.823656    0.390013
 30   0.334152   0.136551    0.183555    0.941133   0.450484    0.461862

or and interval of indices using the syntax =a:b or (a:b)

ras[Ti(1:10)]

6×6×10 Raster{Float64,3} with dimensions: 
  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
extent: Extent(X = (25, 30), Y = (25, 30), Ti = (DateTime("2001-01-01T00:00:00"), DateTime("2001-10-01T00:00:00")))
missingval: missing
values: [:, :, 1]
     25         26          27          28         29          30
 25   0.9063     0.427328    0.0320967   0.297023   0.0571002   0.891377
 26   0.443494   0.867547    0.350546    0.150155   0.24565     0.711039
 27   0.745673   0.0991336   0.930332    0.893537   0.805931    0.360583
 28   0.512083   0.125287    0.959434    0.354868   0.337824    0.259563
 29   0.253849   0.692209    0.774092    0.131798   0.823656    0.390013
 30   0.334152   0.136551    0.183555    0.941133   0.450484    0.461862
[and 9 more slices...]

Select by value

ras[Ti=At(DateTime(2001))]

6×6 Raster{Float64,2} with dimensions: 
  X Sampled{Int64} 25:1:30 ForwardOrdered Regular Points,
  Y Sampled{Int64} 25:1:30 ForwardOrdered Regular Points
and reference dimensions: 
  Ti Sampled{DateTime} DateTime("2001-01-01T00:00:00"):Month(1):DateTime("2001-01-01T00:00:00") ForwardOrdered Regular Points
extent: Extent(X = (25, 30), Y = (25, 30))
missingval: missing
values:      25         26          27          28         29          30
 25   0.9063     0.427328    0.0320967   0.297023   0.0571002   0.891377
 26   0.443494   0.867547    0.350546    0.150155   0.24565     0.711039
 27   0.745673   0.0991336   0.930332    0.893537   0.805931    0.360583
 28   0.512083   0.125287    0.959434    0.354868   0.337824    0.259563
 29   0.253849   0.692209    0.774092    0.131798   0.823656    0.390013
 30   0.334152   0.136551    0.183555    0.941133   0.450484    0.461862

More options are available, like Near, Contains and Where. For more details go here.

Dimensions can also be used in most Base and Statistics methods like mean and reduce where dims arguments are required. Much of the behaviour is covered in the DimensionalData docs.

See the docs for more details and examples for Rasters.jl.

Data-source abstraction

Rasters provides a standardised interface that allows many source data types to be used with identical syntax.

• Scripts and packages building on Rasters.jl can treat Raster, RasterStack, and RasterSeries as black boxes.

  • The data could hold GeoTiff or NetCDF files, Arrays in memory or CuArrays on the GPU - they will all behave in the same way.
  • RasterStack can be backed by a Netcdf or HDF5 file, or a NamedTuple of Raster holding .tif files, or all Raster in memory.
  • Users do not have to deal with the specifics of spatial file types.
  • Projected lookups with Cylindrical projections can by indexed using other Cylindrical projections by setting the mappedcrs keyword on construction. You don't need to know the underlying projection, the conversion is handled automatically. This means lat/lon EPSG(4326) can be used seamlessly if you need that.

Bugs, errors and making issues for Rasters.jl

Raster data is complicated and there are many places for subtle or not-so-subtle bugs to creep in.

We need bug reports to reduce how often they occur over time. But also, we need issues that are easy to reproduce or it isn't practically possible to fix them.

Because there are so many raster file types and variations of them, most of the time we need the exact file that caused your problem to know how to fix it, and be sure that we have actually fixed it when we are done. So fixing a Rasters.jl bug nearly always involves downloading some file and running some code that breaks with it (if you can trigger the bug without a file, thats great! but its not always possible).

To make an issue we can fix quickly (or at all) there are three key steps:

1. Include the file in an accessible place on web without autentication or any other work on our part, so we can just get it and find your bug. You can put it on a file hosting platform (e.g. google drive, drop box, whatever you use) and share the url.
2. Add a minimum working example to the issue template that first downloads the file, then runs the function that triggers the bug.
3. Paste the complete stack trace of the error it produces, right to the bottom, into the issue template. Then we can be sure we reproduced the same problem.

Good issues are really appreciated, but they do take just a little extra effort with Rasters.jl because of this need for files.

source

# Rasters.AbstractRaster — Type.

AbstractRaster <: DimensionalData.AbstractDimArray

Abstract supertype for objects that wrap an array (or location of an array) and metadata about its contents. It may be memory or hold a FileArray, which holds the filename, and is only opened when required.

AbstractRasters inherit from AbstractDimArray from DimensionalData.jl. They can be indexed as regular Julia arrays or with DimensionalData.jl Dimensions. They will plot as a heatmap in Plots.jl with correct coordinates and labels, even after slicing with getindex or view. getindex on a AbstractRaster will always return a memory-backed Raster.

source

# Rasters.AbstractRasterSeries — Type.

AbstractRasterSeries <: DimensionalData.AbstractDimensionalArray

Abstract supertype for high-level DimensionalArray that hold RasterStacks, Rasters, or the paths they can be loaded from. RasterSeries are indexed with dimensions as with a AbstractRaster. This is useful when you have multiple files containing rasters or stacks of rasters spread over dimensions like time and elevation.

As much as possible, implementations should facilitate loading entire directories and detecting the dimensions from metadata.

This allows syntax like below for a series of stacks of arrays:

RasterSeries[Time(Near(DateTime(2001, 1))][:temp][Y(Between(70, 150)), X(Between(-20,20))] |> plot`

RasterSeries is the concrete implementation.

source

# Rasters.AbstractRasterStack — Type.

AbstractRasterStack

Abstract supertype for objects that hold multiple AbstractRasters that share spatial dimensions.

They are NamedTuple-like structures that may either contain NamedTuple of AbstractRasters, string paths that will load AbstractRasters, or a single path that points to a file containing multiple layers, like NetCDF or HDF5. Use and syntax is similar or identical for all cases.

AbstractRasterStack can hold layers that share some or all of their dimensions. They cannot have the same dimension with different length or spatial extent as another layer.

getindex on an AbstractRasterStack generally returns a memory backed standard Raster. raster[:somelayer] |> plot plots the layers array, while raster[:somelayer, X(1:100), Band(2)] |> plot will plot the subset without loading the whole array.

getindex on an AbstractRasterStack with a key returns another stack with getindex applied to all the arrays in the stack.

source

# Rasters.Band — Type.

Band <: Dimension

Band(val=:)

Band Dimension for multi-band rasters.

Example:

banddim = Band(10:10:100)
# Or
val = A[Band(1)]
# Or
mean(A; dims=Band)

source

# Rasters.Mapped — Type.

Mapped <: AbstractProjected

Mapped(order, span, sampling, crs, mappedcrs)
Mapped(; order=AutoOrder(), span=AutoSpan(), sampling=AutoSampling(), crs=nothing, mappedcrs)

An AbstractSampled Lookup, where the dimension index has been mapped to another projection, usually lat/lon or EPSG(4326). Mapped matches the dimension format commonly used in netcdf files.

Fields and behaviours are identical to Sampled with the addition of crs and mappedcrs fields.

The mapped dimension index will be used as for Sampled, but to save in another format the underlying crs may be used to convert it.

source

# Rasters.Projected — Type.

Projected <: AbstractProjected

Projected(order, span, sampling, crs, mappedcrs)
Projected(; order=AutoOrder(), span=AutoSpan(), sampling=AutoSampling(), crs, mappedcrs=nothing)

An AbstractSampled Lookup with projections attached.

Fields and behaviours are identical to Sampled with the addition of crs and mappedcrs fields.

If both crs and mappedcrs fields contain CRS data (in a GeoFormat wrapper from GeoFormatTypes.jl) the selector inputs and plot axes will be converted from and to the specified mappedcrs projection automatically. A common use case would be to pass mappedcrs=EPSG(4326) to the constructor when loading eg. a GDALarray:

GDALarray(filename; mappedcrs=EPSG(4326))

The underlying crs will be detected by GDAL.

If mappedcrs is not supplied (ie. mappedcrs=nothing), the base index will be shown on plots, and selectors will need to use whatever format it is in.

source

# Rasters.Raster — Type.

Raster <: AbsractRaster

Raster(filepath::String; kw...)
Raster(A::AbstractDimArray; kw...)
Raster(A::AbstractArray, dims; kw...)

A generic AbstractRaster for spatial/raster array data. It may hold memory-backed arrays or FileArray, that simply holds the String path to an unopened file. This will only be opened lazily when it is indexed with getindex or when read(A) is called. Broadcasting, taking a view, reversing and most other methods do not load data from disk: they are applied later, lazily.

An AbatractArray for spatial/raster data.

It may hold memory-backed arrays or, when lazy=true a FileArray that simply holds the String path to an unopened file.

WIth lazy=true the file will be opened lazily when it is indexed with getindex or when read(A) is called. Broadcasting, taking a view, reversing and most other methods will not load data from disk: they are applied later, lazily.

Arguments

• dims: Tuple of Dimensions needed when an AbstractArray is used.

Keywords

• name: a Symbol name for the array, which will also retreive named layers if Raster is used on a multi-layered file like a NetCDF. name becomes the layer name if the Raster is combined into a RasterStack.

• group: the group in the dataset where name can be found. Only needed for nested datasets. A String or Symbol will select a single group. Pairs can also used to access groups at any nested depth, i.e group=:group1 => :group2 => :group3.

• missingval: value reprsenting missing data, normally detected from the file. Set manually when you know the value is not specified or is incorrect. This will not change any values in the raster, it simply assigns which value is treated as missing. To replace all of the missing values in the raster, use replace_missing.

• metadata: Dict or Metadata object for the array, or NoMetadata().

• crs: the coordinate reference system of the objects XDim/YDim dimensions. Only set this if you know the detected crs is incrorrect, or it is not present in the file. The crs is expected to be a GeoFormatTypes.jl CRS or Mixed mode GeoFormat object, like EPSG(4326).

• mappedcrs: the mapped coordinate reference system of the objects XDim/YDim dimensions. for Mapped lookups these are the actual values of the index. For Projected lookups this can be used to index in eg. EPSG(4326) lat/lon values, having it converted automatically. Only set this if the detected mappedcrs in incorrect, or the file does not have a mappedcrs, e.g. a tiff. The mappedcrs is expected to be a GeoFormatTypes.jl CRS or Mixed mode GeoFormat type.

• refdims: Tuple of position Dimensions the array was sliced from, defaulting to (). Usually not needed.

When a filepath String is used:

• dropband: drop single band dimensions when creating stacks from filenames. true by default.
• lazy: A Bool specifying if to load data lazily from disk. false by default.
• replace_missing: replace missingval with missing. This is done lazily if lazy=true. Note that currently for NetCDF and GRIB files replace_missing is always true. In future replace_missing=false will also work for these data sources.
• source: Usually automatically detected from filepath extension. To manually force, a Symbol can be passed :gdal, :netcdf, :grd, :grib. The internal Rasters.Source objects, such as Rasters.GDALsource(), Rasters.GRIBsource() or Rasters.NCDsource() can also be used.
• write: defines the default write keyword value when calling open on the Raster. false by default. Only makes sense to use when lazy=true.

When A is an AbstractDimArray:

• data: can replace the data in an existing AbstractRaster

source

# Rasters.RasterSeries — Type.

RasterSeries <: AbstractRasterSeries

RasterSeries(rasters::AbstractArray{<:AbstractRaster}, dims; [refdims])
RasterSeries(stacks::AbstractArray{<:AbstractRasterStack}, dims; [refdims]) 

RasterSeries(paths::AbstractArray{<:AbstractString}, dims; child, duplicate_first, kw...)
RasterSeries(path:::AbstractString, dims; ext, separator, child, duplicate_first, kw...)

Concrete implementation of AbstractRasterSeries.

A RasterSeries is an array of Rasters or RasterStacks, along some dimension(s).

Existing Raster RasterStack can be wrapped in a RasterSeries, or new files can be loaded from an array of String or from a single String.

A single String can refer to a whole directory, or the name of a series of files in a directory, sharing a common stem. The differnce between the filenames can be used as the lookup for the series.

For example, with some tifs at these paths :

"series_dir/myseries_2001-01-01T00:00:00.tif"
"series_dir/myseries_2002-01-01T00:00:00.tif"

We can load a RasterSeries with a DateTime lookup:

julia> ser = RasterSeries("series_dir/myseries.tif", Ti(DateTime))
2-element RasterSeries{Raster,1} with dimensions: 
  Ti Sampled{DateTime} DateTime[DateTime("2001-01-01T00:00:00"), DateTime("2002-01-01T00:00:00")] ForwardOrdered Irregular Points

The DateTime suffix is parsed from the filenames. Using Ti(Int) would try to parse integers intead.

Just using the directory will also work, unless there are other files mixed in it:

julia> ser = RasterSeries("series_dir", Ti(DateTime))
2-element RasterSeries{Raster,1} with dimensions: 
  Ti Sampled{DateTime} DateTime[DateTime("2001-01-01T00:00:00"), DateTime("2002-01-01T00:00:00")] ForwardOrdered Irregular Points

Arguments

• dims: series dimension/s.

Keywords

When loading a series from a Vector of String paths or a single String path:

• child: constructor of child objects for use when filenames are passed in, can be Raster or RasterStack. Defaults to Raster.
• duplicate_first::Bool: wether to duplicate the dimensions and metadata of the first file with all other files. This can save load time with a large series where dimensions are identical. false by default.
• lazy: A Bool specifying if to load data lazily from disk. false by default.
• kw: keywords passed to the child constructor Raster or RasterStack.

When loading a series from a single String path:

• separator: separator used to split lookup elements from the rest of a filename. '_' by default.

Others:

• refdims: existing reference dimension/s, normally not required.

source

# Rasters.RasterStack — Type.

RasterStack <: AbstrackRasterStack

RasterStack(data...; name, kw...)
RasterStack(data::Union{Vector,Tuple}; name, kw...)
RasterStack(data::NamedTuple; kw...))
RasterStack(data::RasterStack; kw...)
RasterStack(data::Raster; layersfrom=Band, kw...)
RasterStack(filepath::AbstractString; kw...)

Load a file path or a NamedTuple of paths as a RasterStack, or convert arguments, a Vector or NamedTuple of Rasters to RasterStack.

Arguments

• data: A NamedTuple of Rasters or String, or a Vector, Tuple or splatted arguments of Raster. The latter options must pass a name keyword argument.
• filepath: A file (such as netcdf or tif) to be loaded as a stack, or a directory path containing multiple files.

Keywords

• name: Used as stack layer names when a Tuple, Vector or splat of Raster is passed in. Has no effect when NameTuple is used - the NamedTuple keys are the layer names.

• group: the group in the dataset where name can be found. Only needed for nested datasets. A String or Symbol will select a single group. Pairs can also used to access groups at any nested depth, i.e group=:group1 => :group2 => :group3.

• metadata: A Dict or DimensionalData.Metadata object.

• missingval: a single value for all layers or a NamedTuple of missingval for each layer. nothing specifies no missing value.

• crs: the coordinate reference system of the objects XDim/YDim dimensions. Only set this if you know the detected crs is incrorrect, or it is not present in the file. The crs is expected to be a GeoFormatTypes.jl CRS or Mixed mode GeoFormat object, like EPSG(4326).

• mappedcrs: the mapped coordinate reference system of the objects XDim/YDim dimensions. for Mapped lookups these are the actual values of the index. For Projected lookups this can be used to index in eg. EPSG(4326) lat/lon values, having it converted automatically. Only set this if the detected mappedcrs in incorrect, or the file does not have a mappedcrs, e.g. a tiff. The mappedcrs is expected to be a GeoFormatTypes.jl CRS or Mixed mode GeoFormat type.

• refdims: Tuple of Dimension that the stack was sliced from.

For when one or multiple filepaths are used:

• dropband: drop single band dimensions when creating stacks from filenames. true by default.
• lazy: A Bool specifying if to load data lazily from disk. false by default.
• replace_missing: replace missingval with missing. This is done lazily if lazy=true. Note that currently for NetCDF and GRIB files replace_missing is always true. In future replace_missing=false will also work for these data sources.
• source: Usually automatically detected from filepath extension. To manually force, a Symbol can be passed :gdal, :netcdf, :grd, :grib. The internal Rasters.Source objects, such as Rasters.GDALsource(), Rasters.GRIBsource() or Rasters.NCDsource() can also be used.

For when a single Raster is used:

• layersfrom: Dimension to source stack layers from if the file is not already multi-layered. nothing is default, so that a single RasterStack(raster) is a single layered stack. RasterStack(raster; layersfrom=Band) will use the bands as layers.

files = (temp="temp.tif", pressure="pressure.tif", relhum="relhum.tif")
stack = RasterStack(files; mappedcrs=EPSG(4326))
stack[:relhum][Lat(Contains(-37), Lon(Contains(144))

source

# DimensionalData.modify — Method.

modify(f, series::AbstractRasterSeries)

Apply function f to the data of the child object. If the child is an AbstractRasterStack the function will be passed on to its child AbstractRasters.

f must return an idenically sized array.

This method triggers a complete rebuild of all objects, and disk based objects will be transferred to memory.

An example of the usefulnesss of this is for swapping out array backend for an entire series to CuArray from CUDA.jl to copy data to a GPU.

source

# GeoInterface.crs — Method.

crs(x::Raster)

Get the projected coordinate reference system of a Y or X Dimension, or of the Y/X dims of an AbstractRaster.

For Mapped lookup this may be nothing as there may be no projected coordinate reference system at all. See setcrs to set it manually.

source

# Rasters.aggregate — Function.

aggregate(method, object, scale; filename, progress, skipmissing)

Aggregate a Raster, or all arrays in a RasterStack or RasterSeries, by scale using method.

Arguments

• method: a function such as mean or sum that can combine the value of multiple cells to generate the aggregated cell, or a Locus like Start() or Center() that specifies where to sample from in the interval.
• object: Object to aggregate, like AbstractRasterSeries, AbstractStack, AbstractRaster or Dimension.
• scale: the aggregation factor, which can be an integer, a tuple of integers for each dimension, or any Dimension, Selector or Int combination you can usually use in getindex. Using a Selector will determine the scale by the distance from the start of the index.

When the aggregation scale of is larger than the array axis, the length of the axis is used.

Keywords

• skipmissingval: if true, any missingval will be skipped during aggregation, so that only areas of all missing values will be aggregated to missingval. If false, any aggegrated area containing a missingval will be assigned missingval.
• filename: a filename to write to directly, useful for large files.
• suffix: a string or value to append to the filename. A tuple of suffix will be applied to stack layers. keys(stack) are the default.
• progress: show a progress bar, true by default, false to hide.

Example

using Rasters, RasterDataSources, Statistics, Plots
using Rasters: Center
st = read(RasterStack(WorldClim{Climate}; month=1))
ag = aggregate(Center(), st, (Y(20), X(20)); skipmissingval=true, progress=false)
plot(ag)
savefig("docs/build/aggregate_example.png"); nothing
# output

Note: currently it is faster to aggregate over memory-backed arrays. Use read on src before use where required.

source

# Rasters.aggregate! — Method.

aggregate!(method, dst::AbstractRaster, src::AbstractRaster, scale; skipmissingval=false)

Aggregate array src to array dst by scale, using method.

Arguments

• method: a function such as mean or sum that can combine the value of multiple cells to generate the aggregated cell, or a Locus like Start() or Center() that species where to sample from in the interval.
• scale: the aggregation factor, which can be an integer, a tuple of integers for each dimension, or any Dimension, Selector or Int combination you can usually use in getindex. Using a Selector will determine the scale by the distance from the start of the index in the src array.

When the aggregation scale of is larger than the array axis, the length of the axis is used.

Keywords

• progress: show a progress bar.
• skipmissingval: if true, any missingval will be skipped during aggregation, so that only areas of all missing values will be aggregated to missingval. If false, any aggegrated area containing a missingval will be assigned missingval.

Note: currently it is much faster to aggregate over memory-backed arrays. Use read on src before use where required.

source

# Rasters.boolmask — Function.

boolmask(obj::Raster; [missingval])
boolmask(obj; [to, res, size])
boolmask(obj::RasterStack; alllayers = true, kw...)

Create a mask array of Bool values, from another Raster. AbstractRasterStack or AbstractRasterSeries are also accepted, if alllayers is true (the default), a mask is taken for all layers, otherwise only the first layer is used.

The array returned from calling boolmask on a AbstractRaster is a Raster with the same dimensions as the original array and a missingval of false.

Arguments

• obj: a Raster, a GeoInterface.jl geometry, or a vector or table of geometries.

Raster / RasterStack Keywords

• missingval: The missing value of the source array, with default missingval(raster).

Keywords

• alllayers
• to: a Raster, RasterStack, Tuple of Dimension or Extents.Extent. If no to object is provided the extent will be calculated from the geometries, Additionally, when no to object or an Extent is passed for to, the size or res keyword must also be used.
• res: the resolution of the dimensions (often in meters or degrees), a Real or Tuple{<:Real,<:Real}. Only required when to is not used or is an Extents.Extent, and size is not used.
• size: the size of the output array, as a Tuple{Int,Int} or single Int for a square. Only required when to is not used or is an Extents.Extent, and res is not used.
• crs: a crs which will be attached to the resulting raster when to not passed or is an Extent. Otherwise the crs from to is used.
• shape: Force data to be treated as :polygon, :line or :point geometries. using points or lines as polygons may have unexpected results.

• boundary: for polygons, include pixels where the :center is inside the polygon, where the polygon :touches the pixel, or that are completely :inside the polygon. The default is :center.

• threaded: run operations in parallel, false by default. In some circumstances threaded can give large speedups over single-threaded operation. This can be true for complicated geometries written into low-resolution rasters, but may not be for simple geometries with high-resolution rasters. With very large rasters threading may be counter productive due to excessing memory use. Caution should also be used: threaded should not be used in in-place functions wrinting to BitArray or other arrays where race conditions can occur.

• progress: show a progress bar, true by default, false to hide.

And specifically for shape=:polygon:

• boundary: include pixels where the :center is inside the polygon, where the line :touches the pixel, or that are completely :inside inside the polygon. The default is :center.

For tabular data, feature collections and other iterables

• collapse: if true, collapse all geometry masks into a single mask. Otherwise return a Raster with an additional geometry dimension, so that each slice along this axis is the mask of the geometry opbject of each row of the table, feature in the feature collection, or just each geometry in the iterable.

Example

using Rasters, RasterDataSources, ArchGDAL, Plots, Dates
wc = Raster(WorldClim{Climate}, :prec; month=1)
boolmask(wc) |> plot

savefig("docs/build/boolmask_example.png"); nothing

# output

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source

# Rasters.cellsize — Method.

cellsize(x)

Gives the approximate size of each cell in square km. This function works for any projection, using an algorithm for polygons on a sphere. It approximates the true size to about 0.1%, depending on latitude.

Run using ArchGDAL to make this method available.

Arguments

• x: A Raster or a Tuple of X and Y dimensions.

Example

using Rasters, ArchGDAL, Rasters.Lookups
dimz = X(Projected(90.0:10.0:120; sampling=Intervals(Start()), order=ForwardOrdered(), span=Regular(10.0), crs=EPSG(4326))),
       Y(Projected(0.0:10.0:50; sampling=Intervals(Start()), order=ForwardOrdered(), span=Regular(10.0), crs=EPSG(4326)))

cs = cellsize(dimz)

# output
4×6 Raster{Float64,2} with dimensions:
  X Projected{Float64} 90.0:10.0:120.0 ForwardOrdered Regular Intervals{Start} crs: EPSG,
  Y Projected{Float64} 0.0:10.0:50.0 ForwardOrdered Regular Intervals{Start} crs: EPSG
extent: Extent(X = (90.0, 130.0), Y = (0.0, 60.0))
missingval: missing
crs: EPSG:4326
parent:
        0.0       10.0       20.0        30.0        40.0            50.0
  90.0  1.2332e6   1.1952e6   1.12048e6   1.01158e6   8.72085e5  706488.0
 100.0  1.2332e6   1.1952e6   1.12048e6   1.01158e6   8.72085e5  706488.0
 110.0  1.2332e6   1.1952e6   1.12048e6   1.01158e6   8.72085e5  706488.0
 120.0  1.2332e6   1.1952e6   1.12048e6   1.01158e6   8.72085e5  706488.0

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source

# Rasters.classify — Function.

classify(x, pairs; lower=(>=), upper=(<), others=nothing)
classify(x, pairs...; lower, upper, others)

Create a new array with values in x classified by the values in pairs.

pairs can hold tuples fo values (2, 3), a Fix2 function e.g. <=(1), a Tuple of Fix2 e.g. (>=(4), <(7)), or an IntervalSets.jl interval, e.g. 3..9 or OpenInterval(10, 12). pairs can also be a n * 3 matrix where each row is lower bounds, upper bounds, replacement.

If tuples or a Matrix are used, the lower and upper keywords define how the lower and upper boundaries are chosen.

If others is set other values not covered in pairs will be set to that values.

Arguments

• x: a Raster or RasterStack
• pairs: each pair contains a value and a replacement, a tuple of lower and upper range and a replacement, or a Tuple of Fix2 like (>(x), <(y).

Keywords

• lower: Which comparison (< or <=) to use for lower values, if Fix2 are not used.
• upper: Which comparison (> or >=) to use for upper values, if Fix2 are not used.
• others: A value to assign to all values not included in pairs. Passing nothing (the default) will leave them unchanged.

Example

using Rasters, RasterDataSources, ArchGDAL, Plots
A = Raster(WorldClim{Climate}, :tavg; month=1)
classes = <=(15) => 10,
          15..25 => 20,
          25..35 => 30,
          >(35) => 40
classified = classify(A, classes; others=0, missingval=0)
plot(classified; c=:magma)

savefig("docs/build/classify_example.png"); nothing

# output

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source

# Rasters.classify! — Method.

classify!(x, pairs...; lower, upper, others)
classify!(x, pairs; lower, upper, others)

Classify the values of x in-place, by the values in pairs.

If Fix2 is not used, the lower and upper keywords

If others is set other values not covered in pairs will be set to that values.

Arguments

• x: a Raster or RasterStack
• pairs: each pair contains a value and a replacement, a tuple of lower and upper range and a replacement, or a Tuple of Fix2 like (>(x), <(y).

Keywords

• lower: Which comparison (< or <=) to use for lower values, if Fix2 are not used.
• upper: Which comparison (> or >=) to use for upper values, if Fix2 are not used.
• others: A value to assign to all values not included in pairs. Passing nothing (the default) will leave them unchanged.

Example

classify! to disk, with key steps:

• copying a tempory file so we don't write over the RasterDataSources.jl version.
• use open with write=true to open the file with disk-write permissions.
• use Float32 like 10.0f0 for all our replacement values and other, because the file is stored as Float32. Attempting to write some other type will fail.

using Rasters, RasterDataSources, ArchGDAL, Plots
# Download and copy the file
filename = getraster(WorldClim{Climate}, :tavg; month=6)
tempfile = tempname() * ".tif"
cp(filename, tempfile)
# Define classes
classes = (5, 15) => 10,
          (15, 25) => 20,
          (25, 35) => 30,
          >=(35) => 40
# Open the file with write permission
open(Raster(tempfile); write=true) do A
    classify!(A, classes; others=0)
end
# Open it again to plot the changes
plot(Raster(tempfile); c=:magma)

savefig("docs/build/classify_bang_example.png"); nothing

# output

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source

# Rasters.combine — Method.

combine(A::Union{AbstractRaster,AbstractRasterStack,AbstracRasterSeries}, [dims]) => Raster

Combine a RasterSeries along some dimension/s, creating a new Raster or RasterStack, depending on the contents of the series.

If dims are passed, only the specified dimensions will be combined with a RasterSeries returned, unless dims is all the dims in the series.

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source

<a id='Rasters.convertlookup-Tuple{Type{<:Lookup}, AbstractDimArray}' href='#Rasters.convertlookup-Tuple{Type{<:Lookup}, AbstractDimArray}'># Rasters.convertlookup — Method.

convertlookup(dstlookup::Type{<:Lookup}, x)

Convert the dimension lookup between Projected and Mapped. Other dimension lookups pass through unchanged.

This is used to e.g. save a netcdf file to GeoTiff.

source

# Rasters.coverage! — Method.

coverage!(A, geom; [mode, scale])

Calculate the area of a raster covered by GeoInterface.jl compatible geomtry geom, as a fraction.

Each pixel is assigned a grid of points (by default 10 x 10) that are each checked to be inside the geometry. The sum divided by the number of points to give coverage.

In pracice, most pixel coverage is not calculated this way - shortcuts that produce the same result are taken wherever possible.

If geom is an AbstractVector or table, the mode keyword will determine how coverage is combined.

Keywords

• mode: method for combining multiple geometries - union or sum.

  • union (the default) gives the areas covered by all geometries. Usefull in spatial coverage where overlapping regions should not be counted twice. The returned raster will contain Float64 values between 0.0 and 1.0.
  • sum gives the summed total of the areas covered by all geometries, as in taking the sum of running coverage separately on all geometries. The returned values are positive Float64.

  For a single geometry, the mode keyword has no effect - the result is the same. * scale: Integer scale of pixel subdivision. The default of 10 means each pixel has 10 x 10 or 100 points that contribute to coverage. Using 100 means 10,000 points contribute. Performance will decline as scale increases. Memory use will grow by scale^2 when mode=:union. * threaded: run operations in parallel, false by default. In some circumstances threaded can give large speedups over single-threaded operation. This can be true for complicated geometries written into low-resolution rasters, but may not be for simple geometries with high-resolution rasters. With very large rasters threading may be counter productive due to excessing memory use. Caution should also be used: threaded should not be used in in-place functions wrinting to BitArray or other arrays where race conditions can occur. * progress: show a progress bar, true by default, false to hide. * vebose: whether to print messages about potential problems. true by default.

source

# Rasters.coverage — Method.

coverage(mode, geom; [to, res, size, scale, verbose, progress])
coverage(geom; [to, mode, res, size, scale, verbose, progress])

Calculate the area of a raster covered by GeoInterface.jl compatible geomtry geom, as a fraction.

Each pixel is assigned a grid of points (by default 10 x 10) that are each checked to be inside the geometry. The sum divided by the number of points to give coverage.

In pracice, most pixel coverage is not calculated this way - shortcuts that produce the same result are taken wherever possible.

If geom is an AbstractVector or table, the mode keyword will determine how coverage is combined.

Keywords

• mode: method for combining multiple geometries - union or sum.

  • union (the default) gives the areas covered by all geometries. Usefull in spatial coverage where overlapping regions should not be counted twice. The returned raster will contain Float64 values between 0.0 and 1.0.
  • sum gives the summed total of the areas covered by all geometries, as in taking the sum of running coverage separately on all geometries. The returned values are positive Float64.

  For a single geometry, the mode keyword has no effect - the result is the same. * scale: Integer scale of pixel subdivision. The default of 10 means each pixel has 10 x 10 or 100 points that contribute to coverage. Using 100 means 10,000 points contribute. Performance will decline as scale increases. Memory use will grow by scale^2 when mode=:union. * threaded: run operations in parallel, false by default. In some circumstances threaded can give large speedups over single-threaded operation. This can be true for complicated geometries written into low-resolution rasters, but may not be for simple geometries with high-resolution rasters. With very large rasters threading may be counter productive due to excessing memory use. Caution should also be used: threaded should not be used in in-place functions wrinting to BitArray or other arrays where race conditions can occur. * progress: show a progress bar, true by default, false to hide. * vebose: whether to print messages about potential problems. true by default.

• to: a Raster, RasterStack, Tuple of Dimension or Extents.Extent. If no to object is provided the extent will be calculated from the geometries, Additionally, when no to object or an Extent is passed for to, the size or res keyword must also be used.

• size: the size of the output array, as a Tuple{Int,Int} or single Int for a square. Only required when to is not used or is an Extents.Extent, and res is not used.
• res: the resolution of the dimensions (often in meters or degrees), a Real or Tuple{<:Real,<:Real}. Only required when to is not used or is an Extents.Extent, and size is not used.

source

# Rasters.crop — Function.

crop(x; to)
crop(xs...; to)

Crop one or multiple AbstractRaster or AbstractRasterStack x to match the size of the object to, or smallest of any dimensions that are shared.

crop is lazy, using a view into the object rather than alocating new memory.

Keywords

• to: the object to crop to. If no to keyword is passed, the smallest shared area of all xs is used.
• touches: true or false. Whether to use Touches wraper on the object extent. When lines need to be included in e.g. zonal statistics, true should be used.

As crop is lazy, filename and suffix keywords are not used.

Example

Crop to another raster:

using Rasters, RasterDataSources, Plots
evenness = Raster(EarthEnv{HabitatHeterogeneity}, :evenness)
rnge = Raster(EarthEnv{HabitatHeterogeneity}, :range)

# Roughly cut out New Zealand from the evenness raster
nz_bounds = X(165 .. 180), Y(-50 .. -32)
nz_evenness = evenness[nz_bounds...]

# Crop range to match evenness
nz_range = crop(rnge; to=nz_evenness)
plot(nz_range)

savefig("docs/build/nz_crop_example.png")
nothing

# output

Crop to a polygon:

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

# Download a borders shapefile
shapefile_url = "https://github.com/nvkelso/natural-earth-vector/raw/master/10m_cultural/ne_10m_admin_0_countries.shp"
shapefile_name = "boundary.shp"
isfile(shapefile_name) || Downloads.download(shapefile_url, shapefile_name)
shp = Shapefile.Handle(shapefile_name).shapes[6]

evenness = Raster(EarthEnv{HabitatHeterogeneity}, :evenness)
argentina_evenness = crop(evenness; to=shp)
plot(argentina_evenness)

savefig("docs/build/argentina_crop_example.png"); nothing

# output

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source

# Rasters.disaggregate — Function.

disaggregate(method, object, scale; filename, progress, keys)

Disaggregate array, or all arrays in a stack or series, by some scale.

Arguments

• method: a function such as mean or sum that can combine the value of multiple cells to generate the aggregated cell, or a Locus like Start() or Center() that species where to sample from in the interval.
• object: Object to aggregate, like AbstractRasterSeries, AbstractStack, AbstractRaster or a Dimension.
• scale: the aggregation factor, which can be an integer, a tuple of integers for each dimension, or any Dimension, Selector or Int combination you can usually use in getindex. Using a Selector will determine the scale by the distance from the start of the index.

Keywords

• progress: show a progress bar.

Note: currently it is faster to aggregate over memory-backed arrays. Use read on src before use where required.

source

# Rasters.disaggregate! — Method.

disaggregate!(method, dst::AbstractRaster, src::AbstractRaster, filename, scale)

Disaggregate array src to array dst by some scale, using method.

• method: a function such as mean or sum that can combine the value of multiple cells to generate the aggregated cell, or a Locus like Start() or Center() that species where to sample from in the interval.
• scale: the aggregation factor, which can be an integer, a tuple of integers for each dimension, or any Dimension, Selector or Int combination you can usually use in getindex. Using a Selector will determine the scale by the distance from the start of the index in the src array.

Note: currently it is faster to aggregate over memory-backed arrays. Use read on src before use where required.

source

# Rasters.extend — Function.

extend(xs...; [to])
extend(xs; [to])
extend(x::Union{AbstractRaster,AbstractRasterStack}; to, kw...)

Extend one or multiple AbstractRaster to match the area covered by all xs, or by the keyword argument to.

Keywords

• to: the Raster or dims to extend to. If no to keyword is passed, the largest shared area of all xs is used.
• touches: true or false. Whether to use Touches wraper on the object extent. When lines need to be included in e.g. zonal statistics, true shoudle be used.
• filename: a filename to write to directly, useful for large files.
• suffix: a string or value to append to the filename. A tuple of suffix will be applied to stack layers. keys(stack) are the default.

using Rasters, RasterDataSources, Plots
evenness = Raster(EarthEnv{HabitatHeterogeneity}, :evenness)
rnge = Raster(EarthEnv{HabitatHeterogeneity}, :range)

# Roughly cut out South America
sa_bounds = X(-88 .. -32), Y(-57 .. 13)
sa_evenness = evenness[sa_bounds...]

# Extend range to match the whole-world raster
sa_range = extend(sa_evenness; to=rnge)
plot(sa_range)

savefig("docs/build/extend_example.png")
nothing
# output

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source

# Rasters.extract — Function.

extract(x, geoms; atol)

Extracts the value of Raster or RasterStack at given points, returning an iterable of NamedTuple with properties for :geometry and raster or stack layer values.

Note that if objects have more dimensions than the length of the point tuples, sliced arrays or stacks will be returned instead of single values.

Arguments

• x: a Raster or RasterStack to extract values from.
• geoms: GeoInterface.jl compatible geometries, or tables or iterables of geometries.

Keywords

• geometry: include :geometry in retured NamedTuple, true by default.
• index: include :index of the CartesianIndex in retured NamedTuple, false by default.
• names: Tuple of Symbol corresponding to layers of a RasterStack to extract. All layers by default.
• skipmissing: skip missing points automatically.
• atol: a tolorerance for floating point lookup values for when the Lookup contains Points. atol is ignored for Intervals.

Example

Here we extact points matching the occurrence of the Mountain Pygmy Possum, Burramis parvus. This could be used to fit a species distribution model.

using Rasters, RasterDataSources, ArchGDAL, GBIF2, CSV

# Get a stack of BioClim layers, and replace missing values with `missing`
st = RasterStack(WorldClim{BioClim}, (1, 3, 5, 7, 12)) |> replace_missing

# Download some occurrence data
obs = GBIF2.occurrence_search("Burramys parvus"; limit=5, year="2009")

# use `extract` to get values for all layers at each observation point.
# We `collect` to get a `Vector` from the lazy iterator.
extract(st, obs; skipmissing=true)

# output
5-element Vector{NamedTuple{(:geometry, :bio1, :bio3, :bio5, :bio7, :bio12)}}:
 (geometry = (0.21, 40.07), bio1 = 17.077084f0, bio3 = 41.20417f0, bio5 = 30.1f0, bio7 = 24.775f0, bio12 = 446.0f0)
 (geometry = (0.03, 39.97), bio1 = 17.076923f0, bio3 = 39.7983f0, bio5 = 29.638462f0, bio7 = 24.153847f0, bio12 = 441.0f0)
 (geometry = (0.03, 39.97), bio1 = 17.076923f0, bio3 = 39.7983f0, bio5 = 29.638462f0, bio7 = 24.153847f0, bio12 = 441.0f0)
 (geometry = (0.52, 40.37), bio1 = missing, bio3 = missing, bio5 = missing, bio7 = missing, bio12 = missing)
 (geometry = (0.32, 40.24), bio1 = 16.321388f0, bio3 = 41.659454f0, bio5 = 30.029825f0, bio7 = 25.544561f0, bio12 = 480.0f0)

Note: passing in arrays, geometry collections or feature collections containing a mix of points and other geometries has undefined results.

source

# Rasters.mappedbounds — Function.

mappedbounds(x)

Get the bounds converted to the mappedcrs value.

Whithout ArchGDAL loaded, this is just the regular bounds.

source

# Rasters.mappedcrs — Function.

mappedcrs(x)

Get the mapped coordinate reference system for the Y/X dims of an array.

In Projected lookup this is used to convert Selector values form the mappedcrs defined projection to the underlying projection, and to show plot axes in the mapped projection.

In Mapped lookup this is the coordinate reference system of the index values. See setmappedcrs to set it manually.

source

# Rasters.mappedindex — Function.

mappedindex(x)

Get the index value of a dimension converted to the mappedcrs value.

Whithout ArchGDAL loaded, this is just the regular dim value.

source

# Rasters.mask! — Function.

mask!(x; with, missingval=missingval(A))

Mask A by the missing values of with, or by all values outside with if it is a polygon.

If with is a polygon, creates a new array where points falling outside the polygon have been replaced by missingval(A).

Return a new array with values of A masked by the missing values of with, or by a polygon.

Arguments

• x: a Raster or RasterStack.

Keywords

• with: another AbstractRaster, a AbstractVector of Tuple points, or any GeoInterface.jl AbstractGeometry. The coordinate reference system of the point must match crs(A).
• missingval: the missing value to write to A in masked areas, by default missingval(A).

Example

Mask an unmasked AWAP layer with a masked WorldClim layer, by first resampling the mask to match the size and projection.

using Rasters, RasterDataSources, ArchGDAL, Plots, Dates

# Load and plot the file
awap = read(RasterStack(AWAP, (:tmin, :tmax); date=DateTime(2001, 1, 1)))
a = plot(awap; clims=(10, 45), c=:imola)

# Create a mask my resampling a worldclim file
wc = Raster(WorldClim{Climate}, :prec; month=1)
wc_mask = resample(wc; to=awap)

# Mask
mask!(awap; with=wc_mask)
b = plot(awap; clims=(10, 45))

savefig(a, "docs/build/mask_bang_example_before.png");
savefig(b, "docs/build/mask_bang_example_after.png"); nothing

# output

Before mask!:

After mask!:

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source

# Rasters.mask — Method.

mask(A:AbstractRaster; with, missingval=missingval(A))
mask(x; with)

Return a new array with values of A masked by the missing values of with, or by the shape of with, if with is a geometric object.

Arguments

• x: a Raster or RasterStack

Keywords

• with: an AbstractRaster, or any GeoInterface.jl compatible objects or table. The coordinate reference system of the point must match crs(A).
• missingval: the missing value to use in the returned file.
• filename: a filename to write to directly, useful for large files.
• suffix: a string or value to append to the filename. A tuple of suffix will be applied to stack layers. keys(stack) are the default.

Geometry keywords

These can be used when with is a GeoInterface.jl compatible object:

• shape: Force data to be treated as :polygon, :line or :point geometries. using points or lines as polygons may have unexpected results.
• boundary: for polygons, include pixels where the :center is inside the polygon, where the polygon :touches the pixel, or that are completely :inside the polygon. The default is :center.

Example

Mask an unmasked AWAP layer with a masked WorldClim layer, by first resampling the mask.

using Rasters, RasterDataSources, ArchGDAL, Plots, Dates

# Load and plot the file
awap = read(Raster(AWAP, :tmax; date=DateTime(2001, 1, 1)))
a = plot(awap; clims=(10, 45))

# Create a mask my resampling a worldclim file
wc = Raster(WorldClim{Climate}, :prec; month=1)
wc_mask = resample(wc; to=awap)

# Mask
awap_masked = mask(awap; with=wc_mask)
b = plot(awap_masked; clims=(10, 45))

savefig(a, "docs/build/mask_example_before.png");
savefig(b, "docs/build/mask_example_after.png"); nothing
# output

Before mask:

After mask:

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source

# Rasters.missingmask — Method.

missingmask(obj::Raster; kw...)
missingmask(obj; [to, res, size, collapse])
missingmask(obj::RasterStack; alllayers = true, kw...)

Create a mask array of missing and true values, from another Raster. AbstractRasterStack or AbstractRasterSeries are also accepted, if alllayers is true (the default), a mask is taken for all layers, otherwise only the first layer is used.

For AbstractRaster the default missingval is missingval(A), but others can be chosen manually.

The array returned from calling missingmask on a AbstractRaster is a Raster with the same size and fields as the original array.

Keywords

• to: a Raster, RasterStack, Tuple of Dimension or Extents.Extent. If no to object is provided the extent will be calculated from the geometries, Additionally, when no to object or an Extent is passed for to, the size or res keyword must also be used.
• res: the resolution of the dimensions (often in meters or degrees), a Real or Tuple{<:Real,<:Real}. Only required when to is not used or is an Extents.Extent, and size is not used.
• size: the size of the output array, as a Tuple{Int,Int} or single Int for a square. Only required when to is not used or is an Extents.Extent, and res is not used.
• crs: a crs which will be attached to the resulting raster when to not passed or is an Extent. Otherwise the crs from to is used.
• shape: Force data to be treated as :polygon, :line or :point geometries. using points or lines as polygons may have unexpected results.
• boundary: for polygons, include pixels where the :center is inside the polygon, where the polygon :touches the pixel, or that are completely :inside the polygon. The default is :center.

Example

using Rasters, RasterDataSources, ArchGDAL, Plots, Dates
wc = Raster(WorldClim{Climate}, :prec; month=1)
missingmask(wc) |> plot

savefig("docs/build/missingmask_example.png"); nothing

# output

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source

# Rasters.missingval — Function.

missingval(x)

Returns the value representing missing data in the dataset

source

# Rasters.mosaic! — Method.

mosaic!(f, x, regions...; missingval, atol)
mosaic!(f, x, regions::Tuple; missingval, atol)

Combine regions in x using the function f.

Arguments

• f a function (e.g. mean, sum, first or last) that is applied to values where regions overlap.
• x: A Raster or RasterStack. May be a an opened disk-based Raster, the result will be written to disk. With the current algorithm, the read speed is slow.
• regions: source objects to be joined. These should be memory-backed (use read first), or may experience poor performance. If all objects have the same extent, mosaic is simply a merge.

Keywords

• missingval: Fills empty areas, and defualts to the `missingval/ of the first layer.
• atol: Absolute tolerance for comparison between index values. This is often required due to minor differences in range values due to floating point error. It is not applied to non-float dimensions. A tuple of tolerances may be passed, matching the dimension order.

Example

Cut out Australia and Africa stacks, then combined them into a single stack.

using Rasters, RasterDataSources, ArchGDAL, Statistics, Plots
st = read(RasterStack(WorldClim{Climate}; month=1))
aus = st[X=100.0 .. 160.0, Y=-50.0 .. -10.0]
africa = st[X=-20.0 .. 60.0, Y=-40.0 .. 35.0]
mosaic!(first, st, aus, africa)
plot(st)
savefig("docs/build/mosaic_bang_example.png")
nothing
# output

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source

# Rasters.mosaic — Method.

mosaic(f, regions...; missingval, atol)
mosaic(f, regions; missingval, atol)

Combine regions into a single raster.

Arguments

• f: A reducing function (mean, sum, first, last etc.) for values where regions overlap.
• regions: Iterable or splatted Raster or RasterStack.

Keywords

• missingval: Fills empty areas, and defualts to the missingval of the first region.
• atol: Absolute tolerance for comparison between index values. This is often required due to minor differences in range values due to floating point error. It is not applied to non-float dimensions. A tuple of tolerances may be passed, matching the dimension order.
• filename: a filename to write to directly, useful for large files.
• suffix: a string or value to append to the filename. A tuple of suffix will be applied to stack layers. keys(stack) are the default.

If your mosaic has has apparent line errors, increase the atol value.

Example

Here we cut out Australia and Africa from a stack, and join them with mosaic.

using Rasters, RasterDataSources, ArchGDAL, Plots
st = RasterStack(WorldClim{Climate}; month=1);

africa = st[X(-20.0 .. 60.0), Y(-40.0 .. 35.0)]
a = plot(africa)

aus = st[X(100.0 .. 160.0), Y(-50.0 .. -10.0)]
b = plot(aus)

# Combine with mosaic
mos = mosaic(first, aus, africa)
c = plot(mos)

savefig(a, "docs/build/mosaic_example_africa.png")
savefig(b, "docs/build/mosaic_example_aus.png")
savefig(c, "docs/build/mosaic_example_combined.png")
nothing
# output

Individual continents

Mosaic of continents

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source

# Rasters.points — Method.

points(A::AbstractRaster; dims=(YDim, XDim), ignore_missing) => Array{Tuple}

Returns a generator of the points in A for dimensions in dims, where points are a tuple of the values in each specified dimension index.

Keywords

• dims the dimensions to return points from. The first slice of other layers will be used.
• ignore_missing: wether to ignore missing values in the array when considering points. If true, all points in the dimensions will be returned, if false only the points that are not === missingval(A) will be returned.

The order of dims determines the order of the points.

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source

# Rasters.rasterize — Function.

rasterize([reducer], data; kw...)

Rasterize a GeoInterface.jl compatable geometry or feature, or a Tables.jl table with a :geometry column of GeoInterface.jl objects, or X, Y points columns.

Arguments

• reducer: a reducing function to reduce the fill value for all geometries that cover or touch a pixel down to a single value. The default is last. Any that takes an iterable and returns a single value will work, including custom functions. However, there are optimisations for built-in methods including sum, first, last, minimum, maximum, extrema and Statistics.mean. These may be an order of magnitude or more faster than count is a special-cased as it does not need a fill value.
• data: a GeoInterface.jl AbstractGeometry, or a nested Vector of AbstractGeometry, or a Tables.jl compatible object containing a :geometry column or points and values columns.

Keywords

These are detected automatically from data where possible.

• to: a Raster, RasterStack, Tuple of Dimension or Extents.Extent. If no to object is provided the extent will be calculated from the geometries, Additionally, when no to object or an Extent is passed for to, the size or res keyword must also be used.
• res: the resolution of the dimensions (often in meters or degrees), a Real or Tuple{<:Real,<:Real}. Only required when to is not used or is an Extents.Extent, and size is not used.
• size: the size of the output array, as a Tuple{Int,Int} or single Int for a square. Only required when to is not used or is an Extents.Extent, and res is not used.
• crs: a crs which will be attached to the resulting raster when to not passed or is an Extent. Otherwise the crs from to is used.
• shape: Force data to be treated as :polygon, :line or :point geometries. using points or lines as polygons may have unexpected results.

• boundary: for polygons, include pixels where the :center is inside the polygon, where the polygon :touches the pixel, or that are completely :inside the polygon. The default is :center.

• fill: the value or values to fill a polygon with. A Symbol or tuple of Symbol will be used to retrieve properties from features or column values from table rows. An array or other iterable will be used for each geometry, in order. fill can also be a function of the current value, e.g. x -> x + 1.

• op: A reducing function that accepts two values and returns one, like min to minimum. For common methods this will be assigned for you, or is not required. But you can use it instead of a reducer as it will usually be faster.
• shape: force data to be treated as :polygon, :line or :point, where possible Points can't be treated as lines or polygons, and lines may not work as polygons, but an attempt will be made.
• geometrycolumn: Symbol to manually select the column the geometries are in when data is a Tables.jl compatible table, or a tuple of Symbol for columns of point coordinates.
• progress: show a progress bar, true by default, false to hide..
• verbose: print information and warnings whne there are problems with the rasterisation. true by default.
• threaded: run operations in parallel, false by default. In some circumstances threaded can give large speedups over single-threaded operation. This can be true for complicated geometries written into low-resolution rasters, but may not be for simple geometries with high-resolution rasters. With very large rasters threading may be counter productive due to excessing memory use. Caution should also be used: threaded should not be used in in-place functions wrinting to BitArray or other arrays where race conditions can occur.
• threadsafe: specify that custom reducer and/or op functions are thread-safe, in that the order of operation or blocking does not matter. For example, sum and maximum are thread-safe, because the answer is approximately (besides floating point error) the same after running on nested blocks, or on all the data. In contrast, median or last are not, because the blocking (median) or order (last) matters.
• filename: a filename to write to directly, useful for large files.
• suffix: a string or value to append to the filename. A tuple of suffix will be applied to stack layers. keys(stack) are the default.

Note on threading. Performance may be much better with threaded=false if reducer/op are not threadsafe. sum, prod, maximum, minimum count and mean (by combining sum and count) are threadsafe. If you know your algorithm is threadsafe, use threadsafe=true to allow all optimisations. Functions passed to fill are always threadsafe, and ignore the threadsafe argument.

Example

Rasterize a shapefile for China and plot, with a border.

using Rasters, RasterDataSources, ArchGDAL, Plots, Dates, Shapefile, Downloads
using Rasters.Lookups

# Download a borders shapefile
shapefile_url = "https://github.com/nvkelso/natural-earth-vector/raw/master/10m_cultural/ne_10m_admin_0_countries.shp"
shapefile_name = "country_borders.shp"
isfile(shapefile_name) || Downloads.download(shapefile_url, shapefile_name)

# Load the shapes for china
china_border = Shapefile.Handle(shapefile_name).shapes[10]

# Rasterize the border polygon
china = rasterize(last, china_border; res=0.1, missingval=0, fill=1, boundary=:touches, progress=false)

# And plot
p = plot(china; color=:spring, legend=false)
plot!(p, china_border; fillalpha=0, linewidth=0.6)

savefig("docs/build/china_rasterized.png"); nothing

# output

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source

# Rasters.rasterize! — Function.

rasterize!([reducer], dest, data; kw...)

Rasterize the geometries in data into the Raster or RasterStack dest, using the values specified by fill.

Arguments

• dest: a Raster or RasterStack to rasterize into.
• reducer: a reducing function to reduce the fill value for all geometries that cover or touch a pixel down to a single value. The default is last. Any that takes an iterable and returns a single value will work, including custom functions. However, there are optimisations for built-in methods including sum, first, last, minimum, maximum, extrema and Statistics.mean. These may be an order of magnitude or more faster than count is a special-cased as it does not need a fill value.
• data: a GeoInterface.jl AbstractGeometry, or a nested Vector of AbstractGeometry, or a Tables.jl compatible object containing a :geometry column or points and values columns.

Keywords

These are detected automatically from A and data where possible.

• fill: the value or values to fill a polygon with. A Symbol or tuple of Symbol will be used to retrieve properties from features or column values from table rows. An array or other iterable will be used for each geometry, in order. fill can also be a function of the current value, e.g. x -> x + 1.
• op: A reducing function that accepts two values and returns one, like min to minimum. For common methods this will be assigned for you, or is not required. But you can use it instead of a reducer as it will usually be faster.
• shape: force data to be treated as :polygon, :line or :point, where possible Points can't be treated as lines or polygons, and lines may not work as polygons, but an attempt will be made.
• geometrycolumn: Symbol to manually select the column the geometries are in when data is a Tables.jl compatible table, or a tuple of Symbol for columns of point coordinates.
• progress: show a progress bar, true by default, false to hide..
• verbose: print information and warnings whne there are problems with the rasterisation. true by default.
• threaded: run operations in parallel, false by default. In some circumstances threaded can give large speedups over single-threaded operation. This can be true for complicated geometries written into low-resolution rasters, but may not be for simple geometries with high-resolution rasters. With very large rasters threading may be counter productive due to excessing memory use. Caution should also be used: threaded should not be used in in-place functions wrinting to BitArray or other arrays where race conditions can occur.
• threadsafe: specify that custom reducer and/or op functions are thread-safe, in that the order of operation or blocking does not matter. For example, sum and maximum are thread-safe, because the answer is approximately (besides floating point error) the same after running on nested blocks, or on all the data. In contrast, median or last are not, because the blocking (median) or order (last) matters.
• to: a Raster, RasterStack, Tuple of Dimension or Extents.Extent. If no to object is provided the extent will be calculated from the geometries, Additionally, when no to object or an Extent is passed for to, the size or res keyword must also be used.
• res: the resolution of the dimensions (often in meters or degrees), a Real or Tuple{<:Real,<:Real}. Only required when to is not used or is an Extents.Extent, and size is not used.
• size: the size of the output array, as a Tuple{Int,Int} or single Int for a square. Only required when to is not used or is an Extents.Extent, and res is not used.
• crs: a crs which will be attached to the resulting raster when to not passed or is an Extent. Otherwise the crs from to is used.
• shape: Force data to be treated as :polygon, :line or :point geometries. using points or lines as polygons may have unexpected results.
• boundary: for polygons, include pixels where the :center is inside the polygon, where the polygon :touches the pixel, or that are completely :inside the polygon. The default is :center.

Example

using Rasters, RasterDataSources, ArchGDAL, Plots, Dates, Shapefile, GeoInterface, Downloads
using Rasters.Lookups

# Download a borders shapefile
shapefile_url = "https://github.com/nvkelso/natural-earth-vector/raw/master/10m_cultural/ne_10m_admin_0_countries.shp"
shapefile_name = "country_borders.shp"
isfile(shapefile_name) || Downloads.download(shapefile_url, shapefile_name)

# Load the shapes for indonesia
indonesia_border = Shapefile.Handle(shapefile_name).shapes[1]

# Make an empty EPSG 4326 projected Raster of the area of Indonesia
dimz = X(Projected(90.0:0.1:145; sampling=Intervals(), crs=EPSG(4326))),
       Y(Projected(-15.0:0.1:10.9; sampling=Intervals(), crs=EPSG(4326)))

A = zeros(UInt32, dimz; missingval=UInt32(0))

# Rasterize each indonesian island with a different number. The islands are
# rings of a multi-polygon, so we use `GI.getring` to get them all separately.
islands = collect(GeoInterface.getring(indonesia_border))
rasterize!(last, A, islands; fill=1:length(islands), progress=false)

# And plot
p = plot(Rasters.trim(A); color=:spring)
plot!(p, indonesia_border; fillalpha=0, linewidth=0.7)

savefig("docs/build/indonesia_rasterized.png"); nothing

# output

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source

# Rasters.replace_missing — Method.

replace_missing(a::AbstractRaster, newmissingval)
replace_missing(a::AbstractRasterStack, newmissingval)

Replace missing values in the array or stack with a new missing value, also updating the missingval field/s.

Keywords

• filename: a filename to write to directly, useful for large files.
• suffix: a string or value to append to the filename. A tuple of suffix will be applied to stack layers. keys(stack) are the default.

Example

using Rasters, RasterDataSources, ArchGDAL
A = Raster(WorldClim{Climate}, :prec; month=1) |> replace_missing
missingval(A)
# output
missing

source

# Rasters.reproject — Method.

reproject(source::GeoFormat, target::GeoFormat, dim::Dimension, val)

reproject uses ArchGDAL.reproject, but implemented for a reprojecting a value array of values, a single dimension at a time.

source

# Rasters.reproject — Method.

reproject(obj; crs)

Reproject the lookups of obj to a different crs.

This is a lossless operation for the raster data, as only the lookup values change. This is only possible when the axes of source and destination projections are alligned: the change is usually from a Regular and an Irregular lookup spans.

For converting between projections that are rotated, skewed or warped in any way, use resample.

Dimensions without an AbstractProjected lookup (such as a Ti dimension) are silently returned without modification.

Arguments

• obj: a Lookup, Dimension, Tuple of Dimension, Raster or RasterStack.
• crs: a crs which will be attached to the resulting raster when to not passed or is an Extent. Otherwise the crs from to is used.

source

# Rasters.resample — Method.

resample(x; kw...)
resample(xs...; to=first(xs), kw...)

resample uses warp (which uses GDALs gdalwarp) to resample a Raster or RasterStack to a new resolution and optionally new crs, or to snap to the bounds, resolution and crs of the object to.

Dimensions without an AbstractProjected lookup (such as a Ti dimension) are iteratively resampled with GDAL and joined back into a single array.

If projections can be converted for each axis independently, it may be faster and more accurate to use reproject.

Run using ArchGDAL to make this method available.

Arguments

• x: the object/s to resample.

Keywords

• to: a Raster, RasterStack, Tuple of Dimension or Extents.Extent. If no to object is provided the extent will be calculated from x,
• res: the resolution of the dimensions (often in meters or degrees), a Real or Tuple{<:Real,<:Real}. Only required when to is not used or is an Extents.Extent, and size is not used.
• size: the size of the output array, as a Tuple{Int,Int} or single Int for a square. Only required when to is not used or is an Extents.Extent, and res is not used.
• crs: a crs which will be attached to the resulting raster when to not passed or is an Extent. Otherwise the crs from to is used.

• method: A Symbol or String specifying the method to use for resampling. From the docs for gdalwarp:

  • :near: nearest neighbour resampling (default, fastest algorithm, worst interpolation quality).
  • :bilinear: bilinear resampling.
  • :cubic: cubic resampling.
  • :cubicspline: cubic spline resampling.
  • :lanczos: Lanczos windowed sinc resampling.
  • :average: average resampling, computes the weighted average of all non-NODATA contributing pixels. rms root mean square / quadratic mean of all non-NODATA contributing pixels (GDAL >= 3.3)
  • :mode: mode resampling, selects the value which appears most often of all the sampled points.
  • :max: maximum resampling, selects the maximum value from all non-NODATA contributing pixels.
  • :min: minimum resampling, selects the minimum value from all non-NODATA contributing pixels.
  • :med: median resampling, selects the median value of all non-NODATA contributing pixels.
  • :q1: first quartile resampling, selects the first quartile value of all non-NODATA contributing pixels.
  • :q3: third quartile resampling, selects the third quartile value of all non-NODATA contributing pixels.
  • :sum: compute the weighted sum of all non-NODATA contributing pixels (since GDAL 3.1)

  Where NODATA values are set to missingval. * filename: a filename to write to directly, useful for large files. * suffix: a string or value to append to the filename. A tuple of suffix will be applied to stack layers. keys(stack) are the default.

Note:

• GDAL may cause some unexpected changes in the raster, such as changing the crs type from EPSG to WellKnownText (it will represent the same CRS).

Example

Resample a WorldClim layer to match an EarthEnv layer:

using Rasters, RasterDataSources, ArchGDAL, Plots
A = Raster(WorldClim{Climate}, :prec; month=1)
B = Raster(EarthEnv{HabitatHeterogeneity}, :evenness)

a = plot(A)
b = plot(resample(A; to=B))

savefig(a, "docs/build/resample_example_before.png");
savefig(b, "docs/build/resample_example_after.png"); nothing

# output

Before resample:

After resample:

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source

# Rasters.setcrs — Method.

setcrs(x, crs)

Set the crs of a Raster, RasterStack, Tuple of Dimension, or a Dimension. The crs is expected to be a GeoFormatTypes.jl CRS or Mixed GeoFormat type

source

# Rasters.setmappedcrs — Method.

setmappedcrs(x, crs)

Set the mapped crs of a Raster, a RasterStack, a Tuple of Dimension, or a Dimension. The crs is expected to be a GeoFormatTypes.jl CRS or Mixed GeoFormat type

source

# Rasters.slice — Method.

slice(A::Union{AbstractRaster,AbstractRasterStack,AbstracRasterSeries}, dims) => RasterSeries

Slice views along some dimension/s to obtain a RasterSeries of the slices.

For a Raster or RasterStack this will return a RasterSeries of Raster or RasterStack that are slices along the specified dimensions.

For a RasterSeries, the output is another series where the child objects are sliced and the series dimensions index is now of the child dimensions combined. slice on a RasterSeries with no dimensions will slice along the dimensions shared by both the series and child object.

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source

# Rasters.trim — Method.

trim(x; dims::Tuple, pad::Int)

Trim missingval(x) from x for axes in dims, returning a view of x.

Arguments

• x: A Raster or RasterStack. For stacks, all layers must having missing values for a pixel for it to be trimmed.

Keywords

• dims: By default dims=(XDim, YDim), so that trimming keeps the area of X and Y that contains non-missing values along all other dimensions.
• pad: The trimmed size will be padded by pad on all sides, although padding will not be added beyond the original extent of the array.

As trim is lazy, filename and suffix keywords are not used.

Example

Create trimmed layers of Australian habitat heterogeneity.

using Rasters, RasterDataSources, Plots
layers = (:evenness, :range, :contrast, :correlation)
st = RasterStack(EarthEnv{HabitatHeterogeneity}, layers)

# Roughly cut out australia
ausbounds = X(100 .. 160), Y(-50 .. -10)
aus = st[ausbounds...]
a = plot(aus)

# Trim missing values and plot
b = plot(trim(aus))

savefig(a, "docs/build/trim_example_before.png");
savefig(b, "docs/build/trim_example_after.png"); nothing

# output

Before trim:

After trim:

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source

# Rasters.warp — Method.

warp(A::AbstractRaster, flags::Dict; kw...)

Gives access to the GDALs gdalwarp method given a Dict of flag => value arguments that can be converted to strings, or vectors where multiple space-separated arguments are required.

Arrays with additional dimensions not handled by GDAL (other than X, Y, Band) are sliced, warped, and then combined to match the original array dimensions. These slices will not be written to disk and loaded lazily at this stage - you will need to do that manually if required.

See the gdalwarp docs for a list of arguments.

Run using ArchGDAL to make this method available.

Keywords

• filename: a filename to write to directly, useful for large files.
• suffix: a string or value to append to the filename. A tuple of suffix will be applied to stack layers. keys(stack) are the default.

Any additional keywords are passed to ArchGDAL.Dataset.

Example

This simply resamples the array with the :tr (output file resolution) and :r flags, giving us a pixelated version:

using Rasters, RasterDataSources, Plots
A = Raster(WorldClim{Climate}, :prec; month=1)
a = plot(A)

flags = Dict(
    :tr => [2.0, 2.0],
    :r => :near,
)
b = plot(warp(A, flags))

savefig(a, "docs/build/warp_example_before.png");
savefig(b, "docs/build/warp_example_after.png"); nothing

# output

Before warp:

After warp:

In practise, prefer resample for this. But warp may be more flexible.

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source

# Rasters.zonal — Method.

zonal(f, x::Union{Raster,RasterStack}; of, kw...)

Calculate zonal statistics for the the zone of a Raster or RasterStack covered by the of object/s.

Arguments

• f: any function that reduces an iterable to a single value, such as sum or Statistics.mean
• x: A Raster or RasterStack
• of: A Raster, RasterStack, dim tuple, extent, GeoInterface.jl compatible geometry, Tables.jl compatible table of a :geometry column, or an AbstractVector of any of these objects..

Keywords

These can be used when of is a GeoInterface.jl compatible object:

• shape: Force data to be treated as :polygon, :line or :point, where possible.
• boundary: for polygons, include pixels where the :center is inside the polygon, where the line :touches the pixel, or that are completely :inside inside the polygon. The default is :center.
• progress: show a progress bar, true by default, false to hide..

Example

``jldoctest using Rasters, RasterDataSources, ArchGDAL, Shapefile, DataFrames, Downloads, Statistics, Dates

Download a borders shapefile

neurl = "https://github.com/nvkelso/natural-earth-vector/raw/master/10mcultural/ne10madmin0countries" shpurl, dbfurl = neurl * ".shp", neurl * ".dbf" shpname, dbfname = "countryborders.shp", "countryborders.dbf" isfile(shpname) || Downloads.download(shpurl, shpname) isfile(dbfurl) || Downloads.download(dbfurl, dbfname)

Download and read a raster stack from WorldClim

st = RasterStack(WorldClim{Climate}; month=Jan, lazy=false)

Load the shapes for world countries

countries = Shapefile.Table(shp_name) |> DataFrame

Calculate the january mean of all climate variables for all countries

january_stats = zonal(mean, st; of=countries, boundary=:touches, progress=false) |> DataFrame

Add the country name column (natural earth has some string errors it seems)

insertcols!(january_stats, 1, :country => first.(split.(countries.ADMIN, r"[^A-Za-z ]")))

output

258×8 DataFrame Row │ country tmin tmax tavg prec ⋯ │ SubStrin… Float32 Float32 Float32 Float64 ⋯ ─────┼────────────────────────────────────────────────────────────────────────── 1 │ Indonesia 21.5447 29.1864 25.3656 271.063 ⋯ 2 │ Malaysia 21.3087 28.4291 24.8688 273.381 3 │ Chile 7.24534 17.9263 12.5858 78.1287 4 │ Bolivia 17.2065 27.7454 22.4759 192.542 5 │ Peru 15.0273 25.5504 20.2888 180.007 ⋯ 6 │ Argentina 13.6751 27.6715 20.6732 67.1837 7 │ Dhekelia Sovereign Base Area 5.87126 15.8991 10.8868 76.25 8 │ Cyprus 5.65921 14.6665 10.1622 97.4474 ⋮ │ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱ 252 │ Spratly Islands 25.0 29.2 27.05 70.5 ⋯ 253 │ Clipperton Island 21.5 33.2727 27.4 6.0 254 │ Macao S 11.6694 17.7288 14.6988 28.0 255 │ Ashmore and Cartier Islands NaN NaN NaN NaN 256 │ Bajo Nuevo Bank NaN NaN NaN NaN ⋯ 257 │ Serranilla Bank NaN NaN NaN NaN 258 │ Scarborough Reef NaN NaN NaN NaN 3 columns and 243 rows omitted ```

source

Reference - Internal functions

# Rasters.AbstractProjected — Type.

AbstractProjected <: AbstractSampled

Abstract supertype for projected index lookups.

source

# Rasters.FileArray — Type.

FileArray{S} <: DiskArrays.AbstractDiskArray

Filearray is a DiskArrays.jl AbstractDiskArray. Instead of holding an open object, it just holds a filename string that is opened lazily when it needs to be read.

source

# Rasters.FileStack — Type.

FileStack{S,Na}

FileStack{S,Na}(filename, types, sizes, eachchunk, haschunks, write)

A wrapper object that holds file pointer and size/chunking metadata for a multi-layered stack stored in a single file, typically netcdf or hdf5.

S is a backend type like NCDsource, and Na is a tuple of Symbol keys.

source

# Rasters.OpenStack — Type.

OpenStack{X,K}

OpenStack{X,K}(dataset)

A wrapper for any stack-like opened dataset that can be indexed with Symbol keys to retrieve AbstractArray layers.

OpenStack is usually hidden from users, wrapped in a regular RasterStack passed as the function argument in open(stack) when the stack is contained in a single file.

X is a backend type like NCDsource, and K is a tuple of Symbol keys.

source

# Rasters.RasterDiskArray — Type.

RasterDiskArray <: DiskArrays.AbstractDiskArray

A basic DiskArrays.jl wrapper for objects that don't have one defined yet. When we open a FileArray it is replaced with a RasterDiskArray.

source

# Base.open — Method.

open(f, A::AbstractRaster; write=false)

open is used to open any lazy=true AbstractRaster and do multiple operations on it in a safe way. The write keyword opens the file in write lookup so that it can be altered on disk using e.g. a broadcast.

f is a method that accepts a single argument - an Raster object which is just an AbstractRaster that holds an open disk-based object. Often it will be a do block:

lazy=false (in-memory) rasters will ignore open and pass themselves to f.

# A is an `Raster` wrapping the opened disk-based object.
open(Raster(filepath); write=true) do A
    mask!(A; with=maskfile)
    A[I...] .*= 2
    # ...  other things you need to do with the open file
end

By using a do block to open files we ensure they are always closed again after we finish working with them.

source

# Base.read! — Method.

read!(src::Union{AbstractString,AbstractRaster}, dst::AbstractRaster)
read!(src::Union{AbstractString,AbstractRasterStack}, dst::AbstractRasterStack)
read!(scr::AbstractRasterSeries, dst::AbstractRasterSeries)

read! will copy the data from src to the object dst.

src can be an object or a file-path String.

source

# Base.read — Method.

read(A::AbstractRaster)
read(A::AbstractRasterStack)
read(A::AbstractRasterSeries)

read will move a Rasters.jl object completely to memory.

source

# Base.skipmissing — Method.

skipmissing(itr::Raster)

Returns an iterable over the elements in a Raster object, skipping any values equal to either the missingval or missing.

source

# Base.write — Method.

Base.write(filepath::AbstractString, s::AbstractRasterSeries; kw...)

Write any AbstractRasterSeries to multiple files, guessing the backend from the file extension.

The lookup values of the series will be appended to the filepath (before the extension), separated by underscores.

All keywords are passed through to these Raster and RasterStack methods.

Keywords

• chunks: a NTuple{N,Int} specifying the chunk size for each dimension. To specify only specific dimensions, a Tuple of Dimension wrapping Int or a NamedTuple of Int can be used. Other dimensions will have a chunk size of 1. true can be used to mean: use the original chunk size of the lazy Raster being written or X and Y of 256 by 256. false means dont use chunks at all.
• ext: filename extension such as ".tiff" or ".nc". Used to specify specific files if only a directory path is used.
• force: false by default. If true it force writing to a file destructively, even if it already exists.

• missingval: set the missing value (i.e. FillValue / nodataval) of the written raster, as Julias missing cannot be stored. If not passed in, missingval will be detected from metadata or a default will be chosen.

  For series with RasterStack child objects, this may be a NamedTuple, one for each layer. * source: Usually automatically detected from filepath extension. To manually force, a Symbol can be passed :gdal, :netcdf, :grd, :grib. The internal Rasters.Source objects, such as Rasters.GDALsource(), Rasters.GRIBsource() or Rasters.NCDsource() can also be used. * vebose: whether to print messages about potential problems. true by default.

source

# Base.write — Method.

Base.write(filename::AbstractString, s::AbstractRasterStack; kw...)

Write any AbstractRasterStack to one or multiple files, depending on the backend. Backend is guessed from the filename extension or forced with the source keyword.

If the source can't be saved as a stack-like object, individual array layers will be saved.

Keywords

• chunks: a NTuple{N,Int} specifying the chunk size for each dimension. To specify only specific dimensions, a Tuple of Dimension wrapping Int or a NamedTuple of Int can be used. Other dimensions will have a chunk size of 1. true can be used to mean: use the original chunk size of the lazy Raster being written or X and Y of 256 by 256. false means dont use chunks at all.
• ext: filename extension such as ".tiff" or ".nc". Used to specify specific files if only a directory path is used.
• force: false by default. If true it force writing to a file destructively, even if it already exists.

• missingval: set the missing value (i.e. FillValue / nodataval) of the written raster, as Julias missing cannot be stored. If not passed in, missingval will be detected from metadata or a default will be chosen.

  For RasterStack this may be a NamedTuple, one for each layer. * source: Usually automatically detected from filepath extension. To manually force, a Symbol can be passed :gdal, :netcdf, :grd, :grib. The internal Rasters.Source objects, such as Rasters.GDALsource(), Rasters.GRIBsource() or Rasters.NCDsource() can also be used. * suffix: a string or value to append to the filename. A tuple of suffix will be applied to stack layers. keys(stack) are the default. * vebose: whether to print messages about potential problems. true by default.

Other keyword arguments are passed to the write method for the backend.

NetCDF keywords

• append: If true, the variable of the current Raster will be appended to filename, if it actually exists.
• deflatelevel: Compression level: 0 (default) means no compression and 9 means maximum compression. Each chunk will be compressed individually.
• shuffle: If true, the shuffle filter is activated which can improve the compression ratio.
• checksum: The checksum method can be :fletcher32 or :nochecksum, the default.
• typename: The name of the NetCDF type required for vlen arrays (https://web.archive.org/save/https://www.unidata.ucar.edu/software/netcdf/netcdf-4/newdocs/netcdf-c/nc005fdef005fvlen.html)

GDAL Keywords

• force: false by default. If true it force writing to a file destructively, even if it already exists.
• driver: A GDAL driver name String or a GDAL driver retrieved via ArchGDAL.getdriver(drivername). By default driver is guessed from the filename extension.
• options::Dict{String,String}: A dictionary containing the dataset creation options passed to the driver. For example: Dict("COMPRESS" => "DEFLATE").

Valid options for each specific driver can be found at: https://gdal.org/drivers/raster/index.html

Source comments

R grd/grid files

Write a Raster to a .grd file with a .gri header file. Returns the base of filename with a .grd extension.

GDAL (tiff, and everything else)

Used if you write a Raster with a filename extension that no other backend can write. GDAL is the fallback, and writes a lot of file types, but is not guaranteed to work.

source

# Base.write — Method.

Base.write(filename::AbstractString, A::AbstractRaster; [source], kw...)

Write an AbstractRaster to file, guessing the backend from the file extension or using the source keyword.

Keywords

• chunks: a NTuple{N,Int} specifying the chunk size for each dimension. To specify only specific dimensions, a Tuple of Dimension wrapping Int or a NamedTuple of Int can be used. Other dimensions will have a chunk size of 1. true can be used to mean: use the original chunk size of the lazy Raster being written or X and Y of 256 by 256. false means dont use chunks at all.
• force: false by default. If true it force writing to a file destructively, even if it already exists.
• missingval: set the missing value (i.e. FillValue / nodataval) of the written raster, as Julias missing cannot be stored. If not passed in, missingval will be detected from metadata or a default will be chosen.
• source: Usually automatically detected from filepath extension. To manually force, a Symbol can be passed :gdal, :netcdf, :grd, :grib. The internal Rasters.Source objects, such as Rasters.GDALsource(), Rasters.GRIBsource() or Rasters.NCDsource() can also be used.

Other keyword arguments are passed to the write method for the backend.

NetCDF keywords

• append: If true, the variable of the current Raster will be appended to filename, if it actually exists.
• deflatelevel: Compression level: 0 (default) means no compression and 9 means maximum compression. Each chunk will be compressed individually.
• shuffle: If true, the shuffle filter is activated which can improve the compression ratio.
• checksum: The checksum method can be :fletcher32 or :nochecksum, the default.
• typename: The name of the NetCDF type required for vlen arrays (https://web.archive.org/save/https://www.unidata.ucar.edu/software/netcdf/netcdf-4/newdocs/netcdf-c/nc005fdef005fvlen.html)

GDAL Keywords

• force: false by default. If true it force writing to a file destructively, even if it already exists.
• driver: A GDAL driver name String or a GDAL driver retrieved via ArchGDAL.getdriver(drivername). By default driver is guessed from the filename extension.
• options::Dict{String,String}: A dictionary containing the dataset creation options passed to the driver. For example: Dict("COMPRESS" => "DEFLATE").

Valid options for each specific driver can be found at: https://gdal.org/drivers/raster/index.html

Source comments

R grd/grid files

Write a Raster to a .grd file with a .gri header file. Returns the base of filename with a .grd extension.

GDAL (tiff, and everything else)

Used if you write a Raster with a filename extension that no other backend can write. GDAL is the fallback, and writes a lot of file types, but is not guaranteed to work.

Returns filename.

source

# Base.write — Method.

Base.write(filename::AbstractString, ::Type{GRDsource}, s::AbstractRaster; kw...)

Write a Raster to a .grd file with a .gri header file.

This method is called automatically if you write a Raster with a .grd or .gri extension.

Keywords

• force: false by default. If true it force writing to a file destructively, even if it already exists.

If this method is called directly the extension of filename will be ignored.

Returns the base of filename with a .grd extension.

source

# Rasters.rplot — Method.

Rasters.rplot([position::GridPosition], raster; kw...)

raster may be a Raster (of 2 or 3 dimensions) or a RasterStack whose underlying rasters are 2 dimensional, or 3-dimensional with a singleton (length-1) third dimension.

Keywords

• plottype = Makie.Heatmap: The type of plot. Can be any Makie plot type which accepts a Raster; in practice, Heatmap, Contour, Contourf and Surface are the best bets.
• axistype = Makie.Axis: The type of axis. This can be an Axis, Axis3, LScene, or even a GeoAxis from GeoMakie.jl.
• X = XDim: The X dimension of the raster.
• Y = YDim: The Y dimension of the raster.
• Z = YDim: The Y dimension of the raster.
• draw_colorbar = true: Whether to draw a colorbar for the axis or not.
• colorbar_position = Makie.Right(): Indicates which side of the axis the colorbar should be placed on. Can be Makie.Top(), Makie.Bottom(), Makie.Left(), or Makie.Right().
• colorbar_padding = Makie.automatic: The amound of padding between the colorbar and its axis. If automatic, then this is set to the width of the colorbar.
• title = Makie.automatic: The titles of each plot. If automatic, these are set to the name of the band.
• xlabel = Makie.automatic: The x-label for the axis. If automatic, set to the dimension name of the X-dimension of the raster.
• ylabel = Makie.automatic: The y-label for the axis. If automatic, set to the dimension name of the Y-dimension of the raster.
• colorbarlabel = "": Usually nothing, but here if you need it. Sets the label on the colorbar.
• colormap = nothing: The colormap for the heatmap. This can be set to a vector of colormaps (symbols, strings, cgrads) if plotting a 3D raster or RasterStack.
• colorrange = Makie.automatic: The colormap for the heatmap. This can be set to a vector of (low, high) if plotting a 3D raster or RasterStack.
• nan_color = :transparent: The color which NaN values should take. Default to transparent.

source