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.
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:
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.
An AbstractSampledLookup, 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.
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:
julia
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.
A generic AbstractRaster for spatial/raster array data. It can hold either memory-backed arrays or, if lazy=true, a FileArray, which stores the String path to an unopened file.
If lazy=true, the file 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 will not load data from disk; they will be 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 retrieve the, alphabetically first, named layer if Raster is used on a multi-layered file like a NetCDF. If instead RasterStack is used to read the multi-layered file, by default, all variables will be added to the stack.
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 incorrect, 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
A RasterSeries is an array of Rasters or RasterStacks, along some dimension(s).
Existing RasterRasterStack 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.
We can load a RasterSeries with a DateTime lookup:
julia
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 instead.
Just using the directory will also work, unless there are other files mixed in it:
julia
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.
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 incorrect, 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.
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 aggregated 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.
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 aggregated 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.
Create a mask array of Bool values, from another Raster. AbstractRasterStack or AbstractRasterSeries are also accepted.
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
a Raster or one or multiple geometries. Geometries can be a GeoInterface.jl AbstractGeometry, a nested Vector of AbstractGeometry, or a Tables.jl compatible object containing a :geometry column or points and values columns, in which case geometrycolumn must be specified.
Raster / RasterStack Keywords
invert: invert the mask, so that areas no missing in with are masked, and areas missing in with are masked.
missingval: The missing value of the source array, with default missingval(raster).
Keywords
alllayers: if true a mask is taken for all layers, otherwise only the first layer is used. Defaults to true
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.
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.
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 excessive memory use. Caution should also be used: threaded should not be used in in-place functions writing 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.
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.
Gives the approximate area of each gridcell of x. By assuming the earth is a sphere, it approximates the true size to about 0.1%, depending on latitude.
Run using ArchGDAL to make this method fully available.
method can be Spherical(; radius) (the default) or Planar().
Spherical will compute cell area on the sphere, by transforming all points back to long-lat. You can specify the radius by the radius keyword argument here. By default, this is 6371008.8, the mean radius of the Earth.
Planar will compute cell area in the plane of the CRS you have chosen. Be warned that this will likely be incorrect for non-equal-area projections.
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.
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.
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.
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.
julia
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("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.
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.
If lazy, concatenate lazily. The default is to concatenate lazily for lazy Rasters and eagerly otherwise.
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.
Calculate the area of a raster covered by GeoInterface.jl compatible geometry 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 practice, 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 excessive memory use. Caution should also be used: threaded should not be used in in-place functions writing 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.
Calculate the area of a raster covered by GeoInterface.jl compatible geometry 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 practice, 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 excessive memory use. Caution should also be used: threaded should not be used in in-place functions writing 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.
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.
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.
crop(x; to, touches=false, [geometrycolumn])
+crop(xs...; to)
Crop one or multiple AbstractRaster or AbstractRasterStackx 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 allocating new memory.
Keywords
to: the object to crop to. This can be a Raster or one or multiple geometries. Geometries can be a GeoInterface.jl AbstractGeometry, a nested Vector of AbstractGeometry, or a Tables.jl compatible object containing a :geometry column or points and values columns, in which case geometrycolumn must be specified. 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.
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.
As crop is lazy, filename and suffix keywords are not used.
Example
Crop to another raster:
julia
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("build/nz_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.
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.
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.
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 wrapper 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.
julia
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("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.
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.
data: a GeoInterface.jl AbstractGeometry, a nested Vector of AbstractGeometry, or a Tables.jl compatible object containing a :geometry column or points and values columns, in which case geometrycolumn must be specified.
Keywords
geometry: include :geometry in returned NamedTuple, true by default.
index: include :index of the CartesianIndex in returned NamedTuple, false by default.
name: a Symbol or Tuple of Symbol corresponding to layer/s of a RasterStack to extract. All layers by default.
skipmissing: skip missing points automatically.
atol: a tolerance for floating point lookup values for when the Lookup contains Points. atol is ignored for Intervals.
– 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.
Example
Here we extract points matching the occurrence of the Mountain Pygmy Possum, Burramis parvus. This could be used to fit a species distribution model.
julia
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.
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.
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).
invert: invert the mask, so that areas no missing in with are masked, and areas missing in with are masked.
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.
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.
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).
invert: invert the mask, so that areas no missing in with are masked, and areas missing in with are masked.
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.
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.
Example
Mask an unmasked AWAP layer with a masked WorldClim layer, by first resampling the 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.
Create a mask array of missing and true values, from another Raster. AbstractRasterStack or AbstractRasterSeries are also accepted-
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.
Arguments
obj: a Raster or one or multiple geometries. Geometries can be a GeoInterface.jl AbstractGeometry, a nested Vector of AbstractGeometry, or a Tables.jl compatible object containing a :geometry column or points and values columns, in which case geometrycolumn must be specified.
Keywords
alllayers: if true a mask is taken for all layers, otherwise only the first layer is used. Defaults to true
invert: invert the mask, so that areas no missing in with are masked, and areas missing in with are masked.
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.
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.
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.
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.
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.
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.
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.
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.
Rasterize a GeoInterface.jl compatable geometry or feature, or a Tables.jl table with a :geometry column of GeoInterface.jl objects, or points columns specified by geometrycolumn
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, a nested Vector of AbstractGeometry, or a Tables.jl compatible object containing a :geometry column or points and values columns, in which case geometrycolumn must be specified.
Keywords
These are detected automatically from data where possible.
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.
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 when 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 excessive memory use. Caution should also be used: threaded should not be used in in-place functions writing 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, minimumcount 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.
julia
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("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.
Rasterize the geometries in data into the Raster or RasterStackdest, 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, a nested Vector of AbstractGeometry, or a Tables.jl compatible object containing a :geometry column or points and values columns, in which case geometrycolumn must be specified.
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 when 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 excessive memory use. Caution should also be used: threaded should not be used in in-place functions writing 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
julia
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("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.
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 aligned: 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.
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:
: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:
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.
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 MixedGeoFormat type
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.
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.
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.
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.
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 DimTuple, Extent, a Raster or one or multiple geometries. Geometries can be a GeoInterface.jl AbstractGeometry, a nested Vector of AbstractGeometry, or a Tables.jl compatible object containing a :geometry column or points and values columns, in which case geometrycolumn must be specified.
Keywords
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.
These can be used when of is or contains (a) GeoInterface.jl compatible object(s):
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..
skipmissing: wether to apply f to the result of skipmissing(A) or not. If truef will be passed an iterator over the values, which loses all spatial information. if falsef will be passes a masked Raster or RasterStack, and will be responsible for handling missing values itself. The default value is true.
Example
julia
using Rasters, RasterDataSources, ArchGDAL, Shapefile, DataFrames, Downloads, Statistics, Dates
+
+# Download a borders shapefile
+ne_url = "https://github.com/nvkelso/natural-earth-vector/raw/master/10m_cultural/ne_10m_admin_0_countries"
+shp_url, dbf_url = ne_url * ".shp", ne_url * ".dbf"
+shp_name, dbf_name = "country_borders.shp", "country_borders.dbf"
+isfile(shp_name) || Downloads.download(shp_url, shp_name)
+isfile(dbf_url) || Downloads.download(dbf_url, dbf_name)
+
+# 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
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.
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.
open is used to open any lazy=trueAbstractRaster 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.
julia
# 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.
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 don't 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.
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 don't 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.
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").
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.
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 don't 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.
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").
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.
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 amount 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.
Most base methods work as for regular julia Arrays, such as reverse and rotations like rotl90. Base, statistics and linear algebra methods like mean that take a dims argument can also use the dimension name.
broadcast works lazily from disk when lazy=true, and is only applied when data is directly indexed. Adding a dot to any function will use broadcast over a Raster just like an Array.
set can be used to modify the nested properties of an objects dimensions, that are more difficult to change with rebuild. set works on the principle that dimension properties can only be in one specific field, so we generally don't have to specify which one it is. set will also try to update anything affected by a change you make.
This will set the X axis to specify points, instead of intervals:
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.
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:
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.
An AbstractSampledLookup, 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.
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:
julia
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.
A generic AbstractRaster for spatial/raster array data. It can hold either memory-backed arrays or, if lazy=true, a FileArray, which stores the String path to an unopened file.
If lazy=true, the file 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 will not load data from disk; they will be 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 retrieve the, alphabetically first, named layer if Raster is used on a multi-layered file like a NetCDF. If instead RasterStack is used to read the multi-layered file, by default, all variables will be added to the stack.
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 incorrect, 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
A RasterSeries is an array of Rasters or RasterStacks, along some dimension(s).
Existing RasterRasterStack 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.
We can load a RasterSeries with a DateTime lookup:
julia
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 instead.
Just using the directory will also work, unless there are other files mixed in it:
julia
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.
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 incorrect, 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.
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 aggregated 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.
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 aggregated 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.
Create a mask array of Bool values, from another Raster. AbstractRasterStack or AbstractRasterSeries are also accepted.
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
a Raster or one or multiple geometries. Geometries can be a GeoInterface.jl AbstractGeometry, a nested Vector of AbstractGeometry, or a Tables.jl compatible object containing a :geometry column or points and values columns, in which case geometrycolumn must be specified.
Raster / RasterStack Keywords
invert: invert the mask, so that areas no missing in with are masked, and areas missing in with are masked.
missingval: The missing value of the source array, with default missingval(raster).
Keywords
alllayers: if true a mask is taken for all layers, otherwise only the first layer is used. Defaults to true
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.
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.
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 excessive memory use. Caution should also be used: threaded should not be used in in-place functions writing 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.
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.
Gives the approximate area of each gridcell of x. By assuming the earth is a sphere, it approximates the true size to about 0.1%, depending on latitude.
Run using ArchGDAL to make this method fully available.
method can be Spherical(; radius) (the default) or Planar().
Spherical will compute cell area on the sphere, by transforming all points back to long-lat. You can specify the radius by the radius keyword argument here. By default, this is 6371008.8, the mean radius of the Earth.
Planar will compute cell area in the plane of the CRS you have chosen. Be warned that this will likely be incorrect for non-equal-area projections.
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.
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.
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.
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.
julia
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("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.
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.
If lazy, concatenate lazily. The default is to concatenate lazily for lazy Rasters and eagerly otherwise.
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.
Calculate the area of a raster covered by GeoInterface.jl compatible geometry 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 practice, 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 excessive memory use. Caution should also be used: threaded should not be used in in-place functions writing 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.
Calculate the area of a raster covered by GeoInterface.jl compatible geometry 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 practice, 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 excessive memory use. Caution should also be used: threaded should not be used in in-place functions writing 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.
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.
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.
crop(x; to, touches=false, [geometrycolumn])
+crop(xs...; to)
Crop one or multiple AbstractRaster or AbstractRasterStackx 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 allocating new memory.
Keywords
to: the object to crop to. This can be a Raster or one or multiple geometries. Geometries can be a GeoInterface.jl AbstractGeometry, a nested Vector of AbstractGeometry, or a Tables.jl compatible object containing a :geometry column or points and values columns, in which case geometrycolumn must be specified. 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.
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.
As crop is lazy, filename and suffix keywords are not used.
Example
Crop to another raster:
julia
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("build/nz_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.
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.
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.
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 wrapper 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.
julia
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("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.
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.
data: a GeoInterface.jl AbstractGeometry, a nested Vector of AbstractGeometry, or a Tables.jl compatible object containing a :geometry column or points and values columns, in which case geometrycolumn must be specified.
Keywords
geometry: include :geometry in returned NamedTuple, true by default.
index: include :index of the CartesianIndex in returned NamedTuple, false by default.
name: a Symbol or Tuple of Symbol corresponding to layer/s of a RasterStack to extract. All layers by default.
skipmissing: skip missing points automatically.
atol: a tolerance for floating point lookup values for when the Lookup contains Points. atol is ignored for Intervals.
– 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.
Example
Here we extract points matching the occurrence of the Mountain Pygmy Possum, Burramis parvus. This could be used to fit a species distribution model.
julia
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.
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.
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).
invert: invert the mask, so that areas no missing in with are masked, and areas missing in with are masked.
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.
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.
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).
invert: invert the mask, so that areas no missing in with are masked, and areas missing in with are masked.
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.
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.
Example
Mask an unmasked AWAP layer with a masked WorldClim layer, by first resampling the 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.
Create a mask array of missing and true values, from another Raster. AbstractRasterStack or AbstractRasterSeries are also accepted-
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.
Arguments
obj: a Raster or one or multiple geometries. Geometries can be a GeoInterface.jl AbstractGeometry, a nested Vector of AbstractGeometry, or a Tables.jl compatible object containing a :geometry column or points and values columns, in which case geometrycolumn must be specified.
Keywords
alllayers: if true a mask is taken for all layers, otherwise only the first layer is used. Defaults to true
invert: invert the mask, so that areas no missing in with are masked, and areas missing in with are masked.
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.
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.
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.
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.
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.
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.
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.
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.
Rasterize a GeoInterface.jl compatable geometry or feature, or a Tables.jl table with a :geometry column of GeoInterface.jl objects, or points columns specified by geometrycolumn
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, a nested Vector of AbstractGeometry, or a Tables.jl compatible object containing a :geometry column or points and values columns, in which case geometrycolumn must be specified.
Keywords
These are detected automatically from data where possible.
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.
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 when 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 excessive memory use. Caution should also be used: threaded should not be used in in-place functions writing 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, minimumcount 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.
julia
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("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.
Rasterize the geometries in data into the Raster or RasterStackdest, 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, a nested Vector of AbstractGeometry, or a Tables.jl compatible object containing a :geometry column or points and values columns, in which case geometrycolumn must be specified.
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 when 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 excessive memory use. Caution should also be used: threaded should not be used in in-place functions writing 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
julia
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("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.
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 aligned: 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.
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:
: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:
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.
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 MixedGeoFormat type
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.
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.
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.
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.
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 DimTuple, Extent, a Raster or one or multiple geometries. Geometries can be a GeoInterface.jl AbstractGeometry, a nested Vector of AbstractGeometry, or a Tables.jl compatible object containing a :geometry column or points and values columns, in which case geometrycolumn must be specified.
Keywords
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.
These can be used when of is or contains (a) GeoInterface.jl compatible object(s):
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..
skipmissing: wether to apply f to the result of skipmissing(A) or not. If truef will be passed an iterator over the values, which loses all spatial information. if falsef will be passes a masked Raster or RasterStack, and will be responsible for handling missing values itself. The default value is true.
Example
julia
using Rasters, RasterDataSources, ArchGDAL, Shapefile, DataFrames, Downloads, Statistics, Dates
+
+# Download a borders shapefile
+ne_url = "https://github.com/nvkelso/natural-earth-vector/raw/master/10m_cultural/ne_10m_admin_0_countries"
+shp_url, dbf_url = ne_url * ".shp", ne_url * ".dbf"
+shp_name, dbf_name = "country_borders.shp", "country_borders.dbf"
+isfile(shp_name) || Downloads.download(shp_url, shp_name)
+isfile(dbf_url) || Downloads.download(dbf_url, dbf_name)
+
+# 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
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.
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.
open is used to open any lazy=trueAbstractRaster 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.
julia
# 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.
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 don't 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.
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 don't 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.
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").
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.
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 don't 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.
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").
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.
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 amount 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.
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.
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:
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.
An AbstractSampledLookup, 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.
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:
julia
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.
A generic AbstractRaster for spatial/raster array data. It can hold either memory-backed arrays or, if lazy=true, a FileArray, which stores the String path to an unopened file.
If lazy=true, the file 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 will not load data from disk; they will be 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 retrieve the, alphabetically first, named layer if Raster is used on a multi-layered file like a NetCDF. If instead RasterStack is used to read the multi-layered file, by default, all variables will be added to the stack.
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 incorrect, 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
A RasterSeries is an array of Rasters or RasterStacks, along some dimension(s).
Existing RasterRasterStack 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.
We can load a RasterSeries with a DateTime lookup:
julia
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 instead.
Just using the directory will also work, unless there are other files mixed in it:
julia
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.
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 incorrect, 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.
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 aggregated 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.
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 aggregated 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.
Create a mask array of Bool values, from another Raster. AbstractRasterStack or AbstractRasterSeries are also accepted.
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
a Raster or one or multiple geometries. Geometries can be a GeoInterface.jl AbstractGeometry, a nested Vector of AbstractGeometry, or a Tables.jl compatible object containing a :geometry column or points and values columns, in which case geometrycolumn must be specified.
Raster / RasterStack Keywords
invert: invert the mask, so that areas no missing in with are masked, and areas missing in with are masked.
missingval: The missing value of the source array, with default missingval(raster).
Keywords
alllayers: if true a mask is taken for all layers, otherwise only the first layer is used. Defaults to true
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.
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.
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 excessive memory use. Caution should also be used: threaded should not be used in in-place functions writing 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.
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.
Gives the approximate area of each gridcell of x. By assuming the earth is a sphere, it approximates the true size to about 0.1%, depending on latitude.
Run using ArchGDAL to make this method fully available.
method can be Spherical(; radius) (the default) or Planar().
Spherical will compute cell area on the sphere, by transforming all points back to long-lat. You can specify the radius by the radius keyword argument here. By default, this is 6371008.8, the mean radius of the Earth.
Planar will compute cell area in the plane of the CRS you have chosen. Be warned that this will likely be incorrect for non-equal-area projections.
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.
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.
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.
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.
julia
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("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.
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.
If lazy, concatenate lazily. The default is to concatenate lazily for lazy Rasters and eagerly otherwise.
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.
Calculate the area of a raster covered by GeoInterface.jl compatible geometry 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 practice, 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 excessive memory use. Caution should also be used: threaded should not be used in in-place functions writing 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.
Calculate the area of a raster covered by GeoInterface.jl compatible geometry 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 practice, 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 excessive memory use. Caution should also be used: threaded should not be used in in-place functions writing 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.
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.
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.
crop(x; to, touches=false, [geometrycolumn])
+crop(xs...; to)
Crop one or multiple AbstractRaster or AbstractRasterStackx 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 allocating new memory.
Keywords
to: the object to crop to. This can be a Raster or one or multiple geometries. Geometries can be a GeoInterface.jl AbstractGeometry, a nested Vector of AbstractGeometry, or a Tables.jl compatible object containing a :geometry column or points and values columns, in which case geometrycolumn must be specified. 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.
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.
As crop is lazy, filename and suffix keywords are not used.
Example
Crop to another raster:
julia
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("build/nz_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.
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.
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.
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 wrapper 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.
julia
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("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.
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.
data: a GeoInterface.jl AbstractGeometry, a nested Vector of AbstractGeometry, or a Tables.jl compatible object containing a :geometry column or points and values columns, in which case geometrycolumn must be specified.
Keywords
geometry: include :geometry in returned NamedTuple, true by default.
index: include :index of the CartesianIndex in returned NamedTuple, false by default.
name: a Symbol or Tuple of Symbol corresponding to layer/s of a RasterStack to extract. All layers by default.
skipmissing: skip missing points automatically.
atol: a tolerance for floating point lookup values for when the Lookup contains Points. atol is ignored for Intervals.
– 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.
Example
Here we extract points matching the occurrence of the Mountain Pygmy Possum, Burramis parvus. This could be used to fit a species distribution model.
julia
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.
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.
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).
invert: invert the mask, so that areas no missing in with are masked, and areas missing in with are masked.
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.
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.
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).
invert: invert the mask, so that areas no missing in with are masked, and areas missing in with are masked.
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.
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.
Example
Mask an unmasked AWAP layer with a masked WorldClim layer, by first resampling the 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.
Create a mask array of missing and true values, from another Raster. AbstractRasterStack or AbstractRasterSeries are also accepted-
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.
Arguments
obj: a Raster or one or multiple geometries. Geometries can be a GeoInterface.jl AbstractGeometry, a nested Vector of AbstractGeometry, or a Tables.jl compatible object containing a :geometry column or points and values columns, in which case geometrycolumn must be specified.
Keywords
alllayers: if true a mask is taken for all layers, otherwise only the first layer is used. Defaults to true
invert: invert the mask, so that areas no missing in with are masked, and areas missing in with are masked.
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.
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.
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.
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.
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.
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.
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.
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.
Rasterize a GeoInterface.jl compatable geometry or feature, or a Tables.jl table with a :geometry column of GeoInterface.jl objects, or points columns specified by geometrycolumn
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, a nested Vector of AbstractGeometry, or a Tables.jl compatible object containing a :geometry column or points and values columns, in which case geometrycolumn must be specified.
Keywords
These are detected automatically from data where possible.
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.
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 when 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 excessive memory use. Caution should also be used: threaded should not be used in in-place functions writing 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, minimumcount 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.
julia
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("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.
Rasterize the geometries in data into the Raster or RasterStackdest, 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, a nested Vector of AbstractGeometry, or a Tables.jl compatible object containing a :geometry column or points and values columns, in which case geometrycolumn must be specified.
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 when 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 excessive memory use. Caution should also be used: threaded should not be used in in-place functions writing 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
julia
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("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.
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 aligned: 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.
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:
: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:
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.
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 MixedGeoFormat type
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.
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.
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.
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.
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 DimTuple, Extent, a Raster or one or multiple geometries. Geometries can be a GeoInterface.jl AbstractGeometry, a nested Vector of AbstractGeometry, or a Tables.jl compatible object containing a :geometry column or points and values columns, in which case geometrycolumn must be specified.
Keywords
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.
These can be used when of is or contains (a) GeoInterface.jl compatible object(s):
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..
skipmissing: wether to apply f to the result of skipmissing(A) or not. If truef will be passed an iterator over the values, which loses all spatial information. if falsef will be passes a masked Raster or RasterStack, and will be responsible for handling missing values itself. The default value is true.
Example
julia
using Rasters, RasterDataSources, ArchGDAL, Shapefile, DataFrames, Downloads, Statistics, Dates
+
+# Download a borders shapefile
+ne_url = "https://github.com/nvkelso/natural-earth-vector/raw/master/10m_cultural/ne_10m_admin_0_countries"
+shp_url, dbf_url = ne_url * ".shp", ne_url * ".dbf"
+shp_name, dbf_name = "country_borders.shp", "country_borders.dbf"
+isfile(shp_name) || Downloads.download(shp_url, shp_name)
+isfile(dbf_url) || Downloads.download(dbf_url, dbf_name)
+
+# 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
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.
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.
open is used to open any lazy=trueAbstractRaster 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.
julia
# 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.
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 don't 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.
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 don't 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.
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").
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.
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 don't 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.
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").
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.
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 amount 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.
',5))])])}const Hs=l(x,[["render",Ms]]);export{Ys as __pageData,Hs as default};
diff --git a/previews/PR807/assets/app.hIvj5Ijy.js b/previews/PR807/assets/app.hIvj5Ijy.js
new file mode 100644
index 00000000..2eab900d
--- /dev/null
+++ b/previews/PR807/assets/app.hIvj5Ijy.js
@@ -0,0 +1 @@
+import{R as p}from"./chunks/theme.CJjspAw-.js";import{R as o,a6 as u,a7 as c,a8 as l,a9 as f,aa as d,ab as m,ac as h,ad as g,ae as A,af as v,d as P,u as R,v as w,s as y,ag as C,ah as b,ai as E,a4 as S}from"./chunks/framework.Cl7EIXwS.js";function i(e){if(e.extends){const a=i(e.extends);return{...a,...e,async enhanceApp(t){a.enhanceApp&&await a.enhanceApp(t),e.enhanceApp&&await e.enhanceApp(t)}}}return e}const s=i(p),T=P({name:"VitePressApp",setup(){const{site:e,lang:a,dir:t}=R();return w(()=>{y(()=>{document.documentElement.lang=a.value,document.documentElement.dir=t.value})}),e.value.router.prefetchLinks&&C(),b(),E(),s.setup&&s.setup(),()=>S(s.Layout)}});async function D(){globalThis.__VITEPRESS__=!0;const e=j(),a=_();a.provide(c,e);const t=l(e.route);return a.provide(f,t),a.component("Content",d),a.component("ClientOnly",m),Object.defineProperties(a.config.globalProperties,{$frontmatter:{get(){return t.frontmatter.value}},$params:{get(){return t.page.value.params}}}),s.enhanceApp&&await s.enhanceApp({app:a,router:e,siteData:h}),{app:a,router:e,data:t}}function _(){return g(T)}function j(){let e=o,a;return A(t=>{let n=v(t),r=null;return n&&(e&&(a=n),(e||a===n)&&(n=n.replace(/\.js$/,".lean.js")),r=import(n)),o&&(e=!1),r},s.NotFound)}o&&D().then(({app:e,router:a,data:t})=>{a.go().then(()=>{u(a.route,t.site),e.mount("#app")})});export{D as createApp};
diff --git a/previews/PR807/assets/argentina_crop_example.DlKBBk5m.png b/previews/PR807/assets/argentina_crop_example.DlKBBk5m.png
new file mode 100644
index 00000000..26a5ae55
Binary files /dev/null and b/previews/PR807/assets/argentina_crop_example.DlKBBk5m.png differ
diff --git a/previews/PR807/assets/array_operations.md.CRu0fcQR.js b/previews/PR807/assets/array_operations.md.CRu0fcQR.js
new file mode 100644
index 00000000..8e20eedf
--- /dev/null
+++ b/previews/PR807/assets/array_operations.md.CRu0fcQR.js
@@ -0,0 +1,141 @@
+import{_ as a,c as n,a5 as t,o as e}from"./chunks/framework.Cl7EIXwS.js";const h=JSON.parse('{"title":"","description":"","frontmatter":{},"headers":[],"relativePath":"array_operations.md","filePath":"array_operations.md","lastUpdated":null}'),i={name:"array_operations.md"};function p(l,s,o,u,r,c){return e(),n("div",null,s[0]||(s[0]=[t(`
Most base methods work as for regular julia Arrays, such as reverse and rotations like rotl90. Base, statistics and linear algebra methods like mean that take a dims argument can also use the dimension name.
broadcast works lazily from disk when lazy=true, and is only applied when data is directly indexed. Adding a dot to any function will use broadcast over a Raster just like an Array.
set can be used to modify the nested properties of an objects dimensions, that are more difficult to change with rebuild. set works on the principle that dimension properties can only be in one specific field, so we generally don't have to specify which one it is. set will also try to update anything affected by a change you make.
This will set the X axis to specify points, instead of intervals:
setcrs(A, crs) and setmappedcrs(A, crs) will set the crs value/s of an object to any GeoFormat from GeoFormatTypes.jl.
`,39)]))}const q=a(i,[["render",p]]);export{h as __pageData,q as default};
diff --git a/previews/PR807/assets/array_operations.md.CRu0fcQR.lean.js b/previews/PR807/assets/array_operations.md.CRu0fcQR.lean.js
new file mode 100644
index 00000000..8e20eedf
--- /dev/null
+++ b/previews/PR807/assets/array_operations.md.CRu0fcQR.lean.js
@@ -0,0 +1,141 @@
+import{_ as a,c as n,a5 as t,o as e}from"./chunks/framework.Cl7EIXwS.js";const h=JSON.parse('{"title":"","description":"","frontmatter":{},"headers":[],"relativePath":"array_operations.md","filePath":"array_operations.md","lastUpdated":null}'),i={name:"array_operations.md"};function p(l,s,o,u,r,c){return e(),n("div",null,s[0]||(s[0]=[t(`
Most base methods work as for regular julia Arrays, such as reverse and rotations like rotl90. Base, statistics and linear algebra methods like mean that take a dims argument can also use the dimension name.
broadcast works lazily from disk when lazy=true, and is only applied when data is directly indexed. Adding a dot to any function will use broadcast over a Raster just like an Array.
set can be used to modify the nested properties of an objects dimensions, that are more difficult to change with rebuild. set works on the principle that dimension properties can only be in one specific field, so we generally don't have to specify which one it is. set will also try to update anything affected by a change you make.
This will set the X axis to specify points, instead of intervals:
Rasters.jl uses a number of backends to load raster data. Raster, RasterStack and RasterSeries will detect which backend to use for you, automatically.
R GRD files can be loaded natively, using Julias MMap - which means they are very fast, but are not compressed. They are always 3 dimensional, and have Y, X and Band dimensions.
NetCDF .nc files are loaded using NCDatasets.jl. Layers from files can be loaded as Raster("filename.nc"; name=:layername). Without name the first layer is used. RasterStack("filename.nc") will use all netcdf variables in the file that are not dimensions as layers.
NetCDF layers can have arbitrary dimensions. Known, common dimension names are converted to X, YZ, and Ti, otherwise Dim{:layername} is used. Layers in the same file may also have different dimensions.
NetCDF files still have issues loading directly from disk for some operations. Using read(ncstack) may help.
All files GDAL can access, such as .tiff and .asc files, can be loaded, using ArchGDAL.jl. These are generally best loaded as Raster("filename.tif"), but can be loaded as RasterStack("filename.tif"; layersfrom=Band), taking layers from the Band dimension, which is also the default.
The Soil Moisture Active-Passive dataset provides global layers of soil moisture, temperature and other related data, in a custom HDF5 format. Layers are always 2 dimensional, with Y and X dimensions.
These can be loaded as multi-layered RasterStack("filename.h5"). Individual layers can be loaded as Raster("filename.h5"; name=:layername), without name the first layer is used.
julia
using Rasters
Missing docstring.
Missing docstring for smapseries. Check Documenter's build log for details.
Files can be written to disk in all formats other than SMAP HDF5 using write("filename.ext", A). See the docs for write. They can (with some caveats) be written to different formats than they were loaded in as, providing file-type conversion for spatial data.
Some metadata may be lost in formats that store little metadata, or where metadata conversion has not been completely implemented.
RasterDataSources.jl standardises the download of common raster data sources, with a focus on datasets used in ecology and the environmental sciences. RasterDataSources.jl is tightly integrated into Rasters.jl, so that datsets and keywords can be used directly to download and load data as a Raster, RasterStack, or RasterSeries.
See the docs for Raster, RasterStack and RasterSeries, and the docs for RasterDataSources.getraster for syntax to specify various data sources.
',23)]))}const k=a(r,[["render",n]]);export{m as __pageData,k as default};
diff --git a/previews/PR807/assets/data_sources.md.Dcl5DOui.lean.js b/previews/PR807/assets/data_sources.md.Dcl5DOui.lean.js
new file mode 100644
index 00000000..2f4406ef
--- /dev/null
+++ b/previews/PR807/assets/data_sources.md.Dcl5DOui.lean.js
@@ -0,0 +1,4 @@
+import{_ as a,c as s,a5 as t,o as i}from"./chunks/framework.Cl7EIXwS.js";const o="/Rasters.jl/previews/PR807/assets/itwpihm.CBMwZ7qg.png",m=JSON.parse('{"title":"Data sources","description":"","frontmatter":{},"headers":[],"relativePath":"data_sources.md","filePath":"data_sources.md","lastUpdated":null}'),r={name:"data_sources.md"};function n(d,e,l,c,h,p){return i(),s("div",null,e[0]||(e[0]=[t(`
Rasters.jl uses a number of backends to load raster data. Raster, RasterStack and RasterSeries will detect which backend to use for you, automatically.
R GRD files can be loaded natively, using Julias MMap - which means they are very fast, but are not compressed. They are always 3 dimensional, and have Y, X and Band dimensions.
NetCDF .nc files are loaded using NCDatasets.jl. Layers from files can be loaded as Raster("filename.nc"; name=:layername). Without name the first layer is used. RasterStack("filename.nc") will use all netcdf variables in the file that are not dimensions as layers.
NetCDF layers can have arbitrary dimensions. Known, common dimension names are converted to X, YZ, and Ti, otherwise Dim{:layername} is used. Layers in the same file may also have different dimensions.
NetCDF files still have issues loading directly from disk for some operations. Using read(ncstack) may help.
All files GDAL can access, such as .tiff and .asc files, can be loaded, using ArchGDAL.jl. These are generally best loaded as Raster("filename.tif"), but can be loaded as RasterStack("filename.tif"; layersfrom=Band), taking layers from the Band dimension, which is also the default.
The Soil Moisture Active-Passive dataset provides global layers of soil moisture, temperature and other related data, in a custom HDF5 format. Layers are always 2 dimensional, with Y and X dimensions.
These can be loaded as multi-layered RasterStack("filename.h5"). Individual layers can be loaded as Raster("filename.h5"; name=:layername), without name the first layer is used.
julia
using Rasters
Missing docstring.
Missing docstring for smapseries. Check Documenter's build log for details.
Files can be written to disk in all formats other than SMAP HDF5 using write("filename.ext", A). See the docs for write. They can (with some caveats) be written to different formats than they were loaded in as, providing file-type conversion for spatial data.
Some metadata may be lost in formats that store little metadata, or where metadata conversion has not been completely implemented.
RasterDataSources.jl standardises the download of common raster data sources, with a focus on datasets used in ecology and the environmental sciences. RasterDataSources.jl is tightly integrated into Rasters.jl, so that datsets and keywords can be used directly to download and load data as a Raster, RasterStack, or RasterSeries.
See the docs for Raster, RasterStack and RasterSeries, and the docs for RasterDataSources.getraster for syntax to specify various data sources.
',23)]))}const k=a(r,[["render",n]]);export{m as __pageData,k as default};
diff --git a/previews/PR807/assets/edcrsgk.2YhfOv1i.png b/previews/PR807/assets/edcrsgk.2YhfOv1i.png
new file mode 100644
index 00000000..68177e08
Binary files /dev/null and b/previews/PR807/assets/edcrsgk.2YhfOv1i.png differ
diff --git a/previews/PR807/assets/eqdfnkz.8o5RdXOt.png b/previews/PR807/assets/eqdfnkz.8o5RdXOt.png
new file mode 100644
index 00000000..243878be
Binary files /dev/null and b/previews/PR807/assets/eqdfnkz.8o5RdXOt.png differ
diff --git a/previews/PR807/assets/extend_example.DNJ4wwKN.png b/previews/PR807/assets/extend_example.DNJ4wwKN.png
new file mode 100644
index 00000000..6e71c96d
Binary files /dev/null and b/previews/PR807/assets/extend_example.DNJ4wwKN.png differ
diff --git a/previews/PR807/assets/fzzzwbh.ALj2aDbZ.png b/previews/PR807/assets/fzzzwbh.ALj2aDbZ.png
new file mode 100644
index 00000000..f4cee49b
Binary files /dev/null and b/previews/PR807/assets/fzzzwbh.ALj2aDbZ.png differ
diff --git a/previews/PR807/assets/gbif_wflow.md.B-xemxXG.js b/previews/PR807/assets/gbif_wflow.md.B-xemxXG.js
new file mode 100644
index 00000000..c34cf7f3
--- /dev/null
+++ b/previews/PR807/assets/gbif_wflow.md.B-xemxXG.js
@@ -0,0 +1,67 @@
+import{_ as a,c as i,a5 as n,o as p}from"./chunks/framework.Cl7EIXwS.js";const t="/Rasters.jl/previews/PR807/assets/olovemv.CVwZE_Z9.png",g=JSON.parse('{"title":"","description":"","frontmatter":{},"headers":[],"relativePath":"gbif_wflow.md","filePath":"gbif_wflow.md","lastUpdated":null}'),l={name:"gbif_wflow.md"};function e(h,s,k,d,r,o){return p(),i("div",null,s[0]||(s[0]=[n(`
Load occurrences for the Mountain Pygmy Possum using GBIF.jl
Then extract predictor variables and write to CSV.
julia
using CSV
+predictors = collect(extract(se_aus, coords))
+CSV.write("burramys_parvus_predictors.csv", predictors)
"burramys_parvus_predictors.csv"
Or convert them to a DataFrame:
julia
using DataFrames
+df = DataFrame(predictors)
+df[1:5,:]
`,22)]))}const E=a(l,[["render",e]]);export{g as __pageData,E as default};
diff --git a/previews/PR807/assets/gbif_wflow.md.B-xemxXG.lean.js b/previews/PR807/assets/gbif_wflow.md.B-xemxXG.lean.js
new file mode 100644
index 00000000..c34cf7f3
--- /dev/null
+++ b/previews/PR807/assets/gbif_wflow.md.B-xemxXG.lean.js
@@ -0,0 +1,67 @@
+import{_ as a,c as i,a5 as n,o as p}from"./chunks/framework.Cl7EIXwS.js";const t="/Rasters.jl/previews/PR807/assets/olovemv.CVwZE_Z9.png",g=JSON.parse('{"title":"","description":"","frontmatter":{},"headers":[],"relativePath":"gbif_wflow.md","filePath":"gbif_wflow.md","lastUpdated":null}'),l={name:"gbif_wflow.md"};function e(h,s,k,d,r,o){return p(),i("div",null,s[0]||(s[0]=[n(`
Load occurrences for the Mountain Pygmy Possum using GBIF.jl
Then extract predictor variables and write to CSV.
julia
using CSV
+predictors = collect(extract(se_aus, coords))
+CSV.write("burramys_parvus_predictors.csv", predictors)
"burramys_parvus_predictors.csv"
Or convert them to a DataFrame:
julia
using DataFrames
+df = DataFrame(predictors)
+df[1:5,:]
`,22)]))}const E=a(l,[["render",e]]);export{g as __pageData,E as default};
diff --git a/previews/PR807/assets/get_started.md.CRLYYECF.js b/previews/PR807/assets/get_started.md.CRLYYECF.js
new file mode 100644
index 00000000..8e762bd2
--- /dev/null
+++ b/previews/PR807/assets/get_started.md.CRLYYECF.js
@@ -0,0 +1,95 @@
+import{_ as a,c as i,a5 as n,o as e}from"./chunks/framework.Cl7EIXwS.js";const t="/Rasters.jl/previews/PR807/assets/vjnjunt.DojgcnZi.png",g=JSON.parse('{"title":"Quick start","description":"","frontmatter":{},"headers":[],"relativePath":"get_started.md","filePath":"get_started.md","lastUpdated":null}'),p={name:"get_started.md"};function l(h,s,d,k,r,o){return e(),i("div",null,s[0]||(s[0]=[n(`
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.
julia
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
More options are available, like Near, Contains and Where.
Dimensions
Rasters uses X, Y, and Z dimensions from DimensionalData to represent spatial directions like longitude, latitude and the vertical dimension, and subset data with them. Ti is used for time, and Band represent bands. Other dimensions can have arbitrary names, but will be treated generically. See DimensionalData for more details on how they work.
Lookup Arrays
These specify properties of the index associated with e.g. the X and Y dimension. Rasters.jl defines additional lookup arrays: Projected to handle dimensions with projections, and Mapped where the projection is mapped to another projection like EPSG(4326). Mapped is largely designed to handle NetCDF dimensions, especially with Explicit spans.
Regular getindex (e.g. A[1:100, :]) and view work on all objects just as with an Array. view is always lazy, and reads from disk are deferred until getindex is used. DimensionalData.jlDimensions and Selectors are the other way to subset an object, making use of the objects index to find values at e.g. certain X/Y coordinates. The available selectors are listed here:
Selectors
Description
At(x)
get the index exactly matching the passed in value(s).
Near(x)
get the closest index to the passed in value(s).
Where(f::Function)
filter the array axis by a function of the dimension index values.
a..b/Between(a, b)
get all indices between two values, excluding the high value.
Contains(x)
get indices where the value x falls within an interval.
Info
Use the .. selector to take a view of madagascar:
julia
using Rasters, RasterDataSources
+const RS = Rasters
+using CairoMakie
+CairoMakie.activate!()
+
+A = Raster(WorldClim{BioClim}, 5)
+madagascar = view(A, X(43.25 .. 50.48), Y(-25.61 .. -12.04))
+# Note the space between .. -12
+Makie.plot(madagascar)
',33)]))}const E=a(p,[["render",l]]);export{g as __pageData,E as default};
diff --git a/previews/PR807/assets/get_started.md.CRLYYECF.lean.js b/previews/PR807/assets/get_started.md.CRLYYECF.lean.js
new file mode 100644
index 00000000..8e762bd2
--- /dev/null
+++ b/previews/PR807/assets/get_started.md.CRLYYECF.lean.js
@@ -0,0 +1,95 @@
+import{_ as a,c as i,a5 as n,o as e}from"./chunks/framework.Cl7EIXwS.js";const t="/Rasters.jl/previews/PR807/assets/vjnjunt.DojgcnZi.png",g=JSON.parse('{"title":"Quick start","description":"","frontmatter":{},"headers":[],"relativePath":"get_started.md","filePath":"get_started.md","lastUpdated":null}'),p={name:"get_started.md"};function l(h,s,d,k,r,o){return e(),i("div",null,s[0]||(s[0]=[n(`
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.
julia
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
More options are available, like Near, Contains and Where.
Dimensions
Rasters uses X, Y, and Z dimensions from DimensionalData to represent spatial directions like longitude, latitude and the vertical dimension, and subset data with them. Ti is used for time, and Band represent bands. Other dimensions can have arbitrary names, but will be treated generically. See DimensionalData for more details on how they work.
Lookup Arrays
These specify properties of the index associated with e.g. the X and Y dimension. Rasters.jl defines additional lookup arrays: Projected to handle dimensions with projections, and Mapped where the projection is mapped to another projection like EPSG(4326). Mapped is largely designed to handle NetCDF dimensions, especially with Explicit spans.
Regular getindex (e.g. A[1:100, :]) and view work on all objects just as with an Array. view is always lazy, and reads from disk are deferred until getindex is used. DimensionalData.jlDimensions and Selectors are the other way to subset an object, making use of the objects index to find values at e.g. certain X/Y coordinates. The available selectors are listed here:
Selectors
Description
At(x)
get the index exactly matching the passed in value(s).
Near(x)
get the closest index to the passed in value(s).
Where(f::Function)
filter the array axis by a function of the dimension index values.
a..b/Between(a, b)
get all indices between two values, excluding the high value.
Contains(x)
get indices where the value x falls within an interval.
Info
Use the .. selector to take a view of madagascar:
julia
using Rasters, RasterDataSources
+const RS = Rasters
+using CairoMakie
+CairoMakie.activate!()
+
+A = Raster(WorldClim{BioClim}, 5)
+madagascar = view(A, X(43.25 .. 50.48), Y(-25.61 .. -12.04))
+# Note the space between .. -12
+Makie.plot(madagascar)
',33)]))}const E=a(p,[["render",l]]);export{g as __pageData,E as default};
diff --git a/previews/PR807/assets/huuwpka.CVxd1wxn.png b/previews/PR807/assets/huuwpka.CVxd1wxn.png
new file mode 100644
index 00000000..43cf8b0b
Binary files /dev/null and b/previews/PR807/assets/huuwpka.CVxd1wxn.png differ
diff --git a/previews/PR807/assets/index.md.CYSmGW74.js b/previews/PR807/assets/index.md.CYSmGW74.js
new file mode 100644
index 00000000..d5980064
--- /dev/null
+++ b/previews/PR807/assets/index.md.CYSmGW74.js
@@ -0,0 +1 @@
+import{_ as a,c as t,a5 as s,o as i}from"./chunks/framework.Cl7EIXwS.js";const h=JSON.parse('{"title":"","description":"","frontmatter":{"layout":"home","hero":{"name":"Rasters.jl","tagline":"Manipulating rasterized spatial data","image":{"src":"/logo.png","alt":"Rasters"},"actions":[{"theme":"brand","text":"Get Started","link":"/get_started"},{"theme":"alt","text":"View on Github","link":"https://github.com/rafaqz/Rasters.jl"},{"theme":"alt","text":"API Reference","link":"/api"}]},"features":[{"title":"Rasters.jl","details":"Defines common types and methods for reading, writing and manipulating rasterized spatial data.","link":"/markdown-examples"},{"title":"Data Formats","details":"These currently include raster arrays like GeoTIF and NetCDF, R grd files, multi-layered stacks, and multi-file series of arrays and stacks."}]},"headers":[],"relativePath":"index.md","filePath":"index.md","lastUpdated":null}'),o={name:"index.md"};function n(l,e,r,d,c,p){return i(),t("div",null,e[0]||(e[0]=[s('
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.
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:
julia
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: not yet finished.
Other functionality in extensions:
Raster data downloads, like Worldclim{Climate}: using RasterDataSources
Makie plots: using Makie
Coordinate transformations for gdal rasters: using CoordinateTransformations
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, that's great! but its not always possible).
To make an issue we can fix quickly (or at all) there are three key steps:
Include the file in an accessible place on web without authentication 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.
Add a minimum working example to the issue template that first downloads the file, then runs the function that triggers the bug.
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.
',14)]))}const g=a(o,[["render",n]]);export{h as __pageData,g as default};
diff --git a/previews/PR807/assets/index.md.CYSmGW74.lean.js b/previews/PR807/assets/index.md.CYSmGW74.lean.js
new file mode 100644
index 00000000..d5980064
--- /dev/null
+++ b/previews/PR807/assets/index.md.CYSmGW74.lean.js
@@ -0,0 +1 @@
+import{_ as a,c as t,a5 as s,o as i}from"./chunks/framework.Cl7EIXwS.js";const h=JSON.parse('{"title":"","description":"","frontmatter":{"layout":"home","hero":{"name":"Rasters.jl","tagline":"Manipulating rasterized spatial data","image":{"src":"/logo.png","alt":"Rasters"},"actions":[{"theme":"brand","text":"Get Started","link":"/get_started"},{"theme":"alt","text":"View on Github","link":"https://github.com/rafaqz/Rasters.jl"},{"theme":"alt","text":"API Reference","link":"/api"}]},"features":[{"title":"Rasters.jl","details":"Defines common types and methods for reading, writing and manipulating rasterized spatial data.","link":"/markdown-examples"},{"title":"Data Formats","details":"These currently include raster arrays like GeoTIF and NetCDF, R grd files, multi-layered stacks, and multi-file series of arrays and stacks."}]},"headers":[],"relativePath":"index.md","filePath":"index.md","lastUpdated":null}'),o={name:"index.md"};function n(l,e,r,d,c,p){return i(),t("div",null,e[0]||(e[0]=[s('
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.
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:
julia
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: not yet finished.
Other functionality in extensions:
Raster data downloads, like Worldclim{Climate}: using RasterDataSources
Makie plots: using Makie
Coordinate transformations for gdal rasters: using CoordinateTransformations
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, that's great! but its not always possible).
To make an issue we can fix quickly (or at all) there are three key steps:
Include the file in an accessible place on web without authentication 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.
Add a minimum working example to the issue template that first downloads the file, then runs the function that triggers the bug.
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.
',14)]))}const g=a(o,[["render",n]]);export{h as __pageData,g as default};
diff --git a/previews/PR807/assets/indonesia_rasterized.CAASrLmh.png b/previews/PR807/assets/indonesia_rasterized.CAASrLmh.png
new file mode 100644
index 00000000..5773fea7
Binary files /dev/null and b/previews/PR807/assets/indonesia_rasterized.CAASrLmh.png differ
diff --git a/previews/PR807/assets/inter-italic-cyrillic-ext.r48I6akx.woff2 b/previews/PR807/assets/inter-italic-cyrillic-ext.r48I6akx.woff2
new file mode 100644
index 00000000..b6b603d5
Binary files /dev/null and b/previews/PR807/assets/inter-italic-cyrillic-ext.r48I6akx.woff2 differ
diff --git a/previews/PR807/assets/inter-italic-cyrillic.By2_1cv3.woff2 b/previews/PR807/assets/inter-italic-cyrillic.By2_1cv3.woff2
new file mode 100644
index 00000000..def40a4f
Binary files /dev/null and b/previews/PR807/assets/inter-italic-cyrillic.By2_1cv3.woff2 differ
diff --git a/previews/PR807/assets/inter-italic-greek-ext.1u6EdAuj.woff2 b/previews/PR807/assets/inter-italic-greek-ext.1u6EdAuj.woff2
new file mode 100644
index 00000000..e070c3d3
Binary files /dev/null and b/previews/PR807/assets/inter-italic-greek-ext.1u6EdAuj.woff2 differ
diff --git a/previews/PR807/assets/inter-italic-greek.DJ8dCoTZ.woff2 b/previews/PR807/assets/inter-italic-greek.DJ8dCoTZ.woff2
new file mode 100644
index 00000000..a3c16ca4
Binary files /dev/null and b/previews/PR807/assets/inter-italic-greek.DJ8dCoTZ.woff2 differ
diff --git a/previews/PR807/assets/inter-italic-latin-ext.CN1xVJS-.woff2 b/previews/PR807/assets/inter-italic-latin-ext.CN1xVJS-.woff2
new file mode 100644
index 00000000..2210a899
Binary files /dev/null and b/previews/PR807/assets/inter-italic-latin-ext.CN1xVJS-.woff2 differ
diff --git a/previews/PR807/assets/inter-italic-latin.C2AdPX0b.woff2 b/previews/PR807/assets/inter-italic-latin.C2AdPX0b.woff2
new file mode 100644
index 00000000..790d62dc
Binary files /dev/null and b/previews/PR807/assets/inter-italic-latin.C2AdPX0b.woff2 differ
diff --git a/previews/PR807/assets/inter-italic-vietnamese.BSbpV94h.woff2 b/previews/PR807/assets/inter-italic-vietnamese.BSbpV94h.woff2
new file mode 100644
index 00000000..1eec0775
Binary files /dev/null and b/previews/PR807/assets/inter-italic-vietnamese.BSbpV94h.woff2 differ
diff --git a/previews/PR807/assets/inter-roman-cyrillic-ext.BBPuwvHQ.woff2 b/previews/PR807/assets/inter-roman-cyrillic-ext.BBPuwvHQ.woff2
new file mode 100644
index 00000000..2cfe6153
Binary files /dev/null and b/previews/PR807/assets/inter-roman-cyrillic-ext.BBPuwvHQ.woff2 differ
diff --git a/previews/PR807/assets/inter-roman-cyrillic.C5lxZ8CY.woff2 b/previews/PR807/assets/inter-roman-cyrillic.C5lxZ8CY.woff2
new file mode 100644
index 00000000..e3886dd1
Binary files /dev/null and b/previews/PR807/assets/inter-roman-cyrillic.C5lxZ8CY.woff2 differ
diff --git a/previews/PR807/assets/inter-roman-greek-ext.CqjqNYQ-.woff2 b/previews/PR807/assets/inter-roman-greek-ext.CqjqNYQ-.woff2
new file mode 100644
index 00000000..36d67487
Binary files /dev/null and b/previews/PR807/assets/inter-roman-greek-ext.CqjqNYQ-.woff2 differ
diff --git a/previews/PR807/assets/inter-roman-greek.BBVDIX6e.woff2 b/previews/PR807/assets/inter-roman-greek.BBVDIX6e.woff2
new file mode 100644
index 00000000..2bed1e85
Binary files /dev/null and b/previews/PR807/assets/inter-roman-greek.BBVDIX6e.woff2 differ
diff --git a/previews/PR807/assets/inter-roman-latin-ext.4ZJIpNVo.woff2 b/previews/PR807/assets/inter-roman-latin-ext.4ZJIpNVo.woff2
new file mode 100644
index 00000000..9a8d1e2b
Binary files /dev/null and b/previews/PR807/assets/inter-roman-latin-ext.4ZJIpNVo.woff2 differ
diff --git a/previews/PR807/assets/inter-roman-latin.Di8DUHzh.woff2 b/previews/PR807/assets/inter-roman-latin.Di8DUHzh.woff2
new file mode 100644
index 00000000..07d3c53a
Binary files /dev/null and b/previews/PR807/assets/inter-roman-latin.Di8DUHzh.woff2 differ
diff --git a/previews/PR807/assets/inter-roman-vietnamese.BjW4sHH5.woff2 b/previews/PR807/assets/inter-roman-vietnamese.BjW4sHH5.woff2
new file mode 100644
index 00000000..57bdc22a
Binary files /dev/null and b/previews/PR807/assets/inter-roman-vietnamese.BjW4sHH5.woff2 differ
diff --git a/previews/PR807/assets/itwpihm.CBMwZ7qg.png b/previews/PR807/assets/itwpihm.CBMwZ7qg.png
new file mode 100644
index 00000000..288a96b7
Binary files /dev/null and b/previews/PR807/assets/itwpihm.CBMwZ7qg.png differ
diff --git a/previews/PR807/assets/maoowux.CO0mgiME.png b/previews/PR807/assets/maoowux.CO0mgiME.png
new file mode 100644
index 00000000..5d70e1e6
Binary files /dev/null and b/previews/PR807/assets/maoowux.CO0mgiME.png differ
diff --git a/previews/PR807/assets/methods.md.tEJSj4eC.js b/previews/PR807/assets/methods.md.tEJSj4eC.js
new file mode 100644
index 00000000..d1cd47d2
--- /dev/null
+++ b/previews/PR807/assets/methods.md.tEJSj4eC.js
@@ -0,0 +1 @@
+import{_ as e,c as a,a5 as s,o as r}from"./chunks/framework.Cl7EIXwS.js";const p=JSON.parse('{"title":"","description":"","frontmatter":{},"headers":[],"relativePath":"methods.md","filePath":"methods.md","lastUpdated":null}'),o={name:"methods.md"};function l(d,t,i,n,c,h){return r(),a("div",null,t[0]||(t[0]=[s('
Methods that change the resolution or extent of an object
Click through to the function documentation for more in-depth descriptions and examples.
Note that most regular Julia methods, such as replace, work as for a standard Array. These additional methods are commonly required in GIS applications.
',10)]))}const g=e(o,[["render",l]]);export{p as __pageData,g as default};
diff --git a/previews/PR807/assets/methods.md.tEJSj4eC.lean.js b/previews/PR807/assets/methods.md.tEJSj4eC.lean.js
new file mode 100644
index 00000000..d1cd47d2
--- /dev/null
+++ b/previews/PR807/assets/methods.md.tEJSj4eC.lean.js
@@ -0,0 +1 @@
+import{_ as e,c as a,a5 as s,o as r}from"./chunks/framework.Cl7EIXwS.js";const p=JSON.parse('{"title":"","description":"","frontmatter":{},"headers":[],"relativePath":"methods.md","filePath":"methods.md","lastUpdated":null}'),o={name:"methods.md"};function l(d,t,i,n,c,h){return r(),a("div",null,t[0]||(t[0]=[s('
Methods that change the resolution or extent of an object
Click through to the function documentation for more in-depth descriptions and examples.
Note that most regular Julia methods, such as replace, work as for a standard Array. These additional methods are commonly required in GIS applications.
',10)]))}const g=e(o,[["render",l]]);export{p as __pageData,g as default};
diff --git a/previews/PR807/assets/mosaic_bang_example.ptHNiUCT.png b/previews/PR807/assets/mosaic_bang_example.ptHNiUCT.png
new file mode 100644
index 00000000..974dd8e4
Binary files /dev/null and b/previews/PR807/assets/mosaic_bang_example.ptHNiUCT.png differ
diff --git a/previews/PR807/assets/mosaic_example_africa.Dpr9JnNl.png b/previews/PR807/assets/mosaic_example_africa.Dpr9JnNl.png
new file mode 100644
index 00000000..95db5ce0
Binary files /dev/null and b/previews/PR807/assets/mosaic_example_africa.Dpr9JnNl.png differ
diff --git a/previews/PR807/assets/mosaic_example_aus.3EfcKnQU.png b/previews/PR807/assets/mosaic_example_aus.3EfcKnQU.png
new file mode 100644
index 00000000..694a77d9
Binary files /dev/null and b/previews/PR807/assets/mosaic_example_aus.3EfcKnQU.png differ
diff --git a/previews/PR807/assets/mosaic_example_combined.XY5Q_nfP.png b/previews/PR807/assets/mosaic_example_combined.XY5Q_nfP.png
new file mode 100644
index 00000000..94005ed4
Binary files /dev/null and b/previews/PR807/assets/mosaic_example_combined.XY5Q_nfP.png differ
diff --git a/previews/PR807/assets/motceku.oyWWmvC1.png b/previews/PR807/assets/motceku.oyWWmvC1.png
new file mode 100644
index 00000000..40438a0c
Binary files /dev/null and b/previews/PR807/assets/motceku.oyWWmvC1.png differ
diff --git a/previews/PR807/assets/nuekexd.BRgQ_uqM.png b/previews/PR807/assets/nuekexd.BRgQ_uqM.png
new file mode 100644
index 00000000..fc0327dd
Binary files /dev/null and b/previews/PR807/assets/nuekexd.BRgQ_uqM.png differ
diff --git a/previews/PR807/assets/nz_crop_example.CeBIxUDy.png b/previews/PR807/assets/nz_crop_example.CeBIxUDy.png
new file mode 100644
index 00000000..41a0304f
Binary files /dev/null and b/previews/PR807/assets/nz_crop_example.CeBIxUDy.png differ
diff --git a/previews/PR807/assets/olovemv.CVwZE_Z9.png b/previews/PR807/assets/olovemv.CVwZE_Z9.png
new file mode 100644
index 00000000..0f46062e
Binary files /dev/null and b/previews/PR807/assets/olovemv.CVwZE_Z9.png differ
diff --git a/previews/PR807/assets/ouepwku.BPZOvOxi.png b/previews/PR807/assets/ouepwku.BPZOvOxi.png
new file mode 100644
index 00000000..3f00b7e3
Binary files /dev/null and b/previews/PR807/assets/ouepwku.BPZOvOxi.png differ
diff --git a/previews/PR807/assets/plot_makie.md.DYZzsleT.js b/previews/PR807/assets/plot_makie.md.DYZzsleT.js
new file mode 100644
index 00000000..7199f260
--- /dev/null
+++ b/previews/PR807/assets/plot_makie.md.DYZzsleT.js
@@ -0,0 +1,78 @@
+import{_ as n,c as e,a5 as a,j as i,a as l,G as h,B as p,o as k}from"./chunks/framework.Cl7EIXwS.js";const r="/Rasters.jl/previews/PR807/assets/fzzzwbh.ALj2aDbZ.png",d="/Rasters.jl/previews/PR807/assets/vxnkhbp.Dw_eTAB7.png",o="/Rasters.jl/previews/PR807/assets/swkcnzv.B2QjG_n4.png",E="/Rasters.jl/previews/PR807/assets/rplot.CwrU8Sen.mp4",g="/Rasters.jl/previews/PR807/assets/aus_trim.B4Z7jnS4.png",D=JSON.parse('{"title":"reset theme","description":"","frontmatter":{},"headers":[],"relativePath":"plot_makie.md","filePath":"plot_makie.md","lastUpdated":null}'),c={name:"plot_makie.md"},y={class:"jldocstring custom-block",open:""};function F(u,s,C,m,b,f){const t=p("Badge");return k(),e("div",null,[s[3]||(s[3]=a(`
Plotting in Makie works somewhat differently than Plots, since the recipe system is different. You can pass a 2-D raster to any surface-like function (heatmap, contour, contourf, or even surface for a 3D plot) with ease.
fig, ax, _ = plot(A)
+contour(fig[1, 2], A)
+ax = Axis(fig[2, 1]; aspect = DataAspect())
+contourf!(ax, A)
+surface(fig[2, 2], A) # even a 3D plot work!
+fig
Rasters.rplot should support Observable input out of the box, but the dimensions of that input must remain the same - i.e., the element names of a RasterStack must remain the same.
julia
Makie.set_theme!(Rasters.theme_rasters())
+# \`stack\` is the WorldClim climate data for January
+stack_obs = Observable(stack)
+fig = Rasters.rplot(stack_obs;
+ Colorbar=(; height=Relative(0.75), width=5)
+)
+record(fig, "rplot.mp4", 1:12; framerate = 3) do i
+ stack_obs[] = RasterStack(WorldClim{Climate}; month = i)
+end
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 amount 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.
# colorbar attributes
+colormap = :batlow
+flipaxis = false
+tickalign=1
+width = 13
+ticksize = 13
+# figure
+with_theme(theme_dark()) do
+ fig = Figure(; size=(600, 600), backgroundcolor=:transparent)
+ axs = [Axis(fig[i,j], xlabel = "lon", ylabel = "lat",
+ backgroundcolor=:transparent) for i in 1:2 for j in 1:2]
+ plt = [Makie.heatmap!(axs[i], aus[l]; colormap) for (i, l) in enumerate(layers)]
+ for (i, l) in enumerate(layers) axs[i].title = string(l) end
+ hidexdecorations!.(axs[1:2]; grid=false, ticks=false)
+ hideydecorations!.(axs[[2,4]]; grid=false, ticks=false)
+ Colorbar(fig[1, 0], plt[1]; flipaxis, tickalign, width, ticksize)
+ Colorbar(fig[1, 3], plt[2]; tickalign, width, ticksize)
+ Colorbar(fig[2, 0], plt[3]; flipaxis, tickalign, width, ticksize)
+ Colorbar(fig[2, 3], plt[4]; tickalign, width, ticksize)
+ colgap!(fig.layout, 5)
+ rowgap!(fig.layout, 5)
+ Label(fig[0, :], "RasterStack of EarthEnv HabitatHeterogeneity layers, trimmed to Australia")
+ fig
+end
+save("aus_trim.png", current_figure());
CairoMakie.Screen{IMAGE}
',9))])}const v=n(c,[["render",F]]);export{D as __pageData,v as default};
diff --git a/previews/PR807/assets/plot_makie.md.DYZzsleT.lean.js b/previews/PR807/assets/plot_makie.md.DYZzsleT.lean.js
new file mode 100644
index 00000000..7199f260
--- /dev/null
+++ b/previews/PR807/assets/plot_makie.md.DYZzsleT.lean.js
@@ -0,0 +1,78 @@
+import{_ as n,c as e,a5 as a,j as i,a as l,G as h,B as p,o as k}from"./chunks/framework.Cl7EIXwS.js";const r="/Rasters.jl/previews/PR807/assets/fzzzwbh.ALj2aDbZ.png",d="/Rasters.jl/previews/PR807/assets/vxnkhbp.Dw_eTAB7.png",o="/Rasters.jl/previews/PR807/assets/swkcnzv.B2QjG_n4.png",E="/Rasters.jl/previews/PR807/assets/rplot.CwrU8Sen.mp4",g="/Rasters.jl/previews/PR807/assets/aus_trim.B4Z7jnS4.png",D=JSON.parse('{"title":"reset theme","description":"","frontmatter":{},"headers":[],"relativePath":"plot_makie.md","filePath":"plot_makie.md","lastUpdated":null}'),c={name:"plot_makie.md"},y={class:"jldocstring custom-block",open:""};function F(u,s,C,m,b,f){const t=p("Badge");return k(),e("div",null,[s[3]||(s[3]=a(`
Plotting in Makie works somewhat differently than Plots, since the recipe system is different. You can pass a 2-D raster to any surface-like function (heatmap, contour, contourf, or even surface for a 3D plot) with ease.
fig, ax, _ = plot(A)
+contour(fig[1, 2], A)
+ax = Axis(fig[2, 1]; aspect = DataAspect())
+contourf!(ax, A)
+surface(fig[2, 2], A) # even a 3D plot work!
+fig
Rasters.rplot should support Observable input out of the box, but the dimensions of that input must remain the same - i.e., the element names of a RasterStack must remain the same.
julia
Makie.set_theme!(Rasters.theme_rasters())
+# \`stack\` is the WorldClim climate data for January
+stack_obs = Observable(stack)
+fig = Rasters.rplot(stack_obs;
+ Colorbar=(; height=Relative(0.75), width=5)
+)
+record(fig, "rplot.mp4", 1:12; framerate = 3) do i
+ stack_obs[] = RasterStack(WorldClim{Climate}; month = i)
+end
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 amount 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.
# colorbar attributes
+colormap = :batlow
+flipaxis = false
+tickalign=1
+width = 13
+ticksize = 13
+# figure
+with_theme(theme_dark()) do
+ fig = Figure(; size=(600, 600), backgroundcolor=:transparent)
+ axs = [Axis(fig[i,j], xlabel = "lon", ylabel = "lat",
+ backgroundcolor=:transparent) for i in 1:2 for j in 1:2]
+ plt = [Makie.heatmap!(axs[i], aus[l]; colormap) for (i, l) in enumerate(layers)]
+ for (i, l) in enumerate(layers) axs[i].title = string(l) end
+ hidexdecorations!.(axs[1:2]; grid=false, ticks=false)
+ hideydecorations!.(axs[[2,4]]; grid=false, ticks=false)
+ Colorbar(fig[1, 0], plt[1]; flipaxis, tickalign, width, ticksize)
+ Colorbar(fig[1, 3], plt[2]; tickalign, width, ticksize)
+ Colorbar(fig[2, 0], plt[3]; flipaxis, tickalign, width, ticksize)
+ Colorbar(fig[2, 3], plt[4]; tickalign, width, ticksize)
+ colgap!(fig.layout, 5)
+ rowgap!(fig.layout, 5)
+ Label(fig[0, :], "RasterStack of EarthEnv HabitatHeterogeneity layers, trimmed to Australia")
+ fig
+end
+save("aus_trim.png", current_figure());
CairoMakie.Screen{IMAGE}
',9))])}const v=n(c,[["render",F]]);export{D as __pageData,v as default};
diff --git a/previews/PR807/assets/plotting.md.n1eGQxVQ.js b/previews/PR807/assets/plotting.md.n1eGQxVQ.js
new file mode 100644
index 00000000..288b22f5
--- /dev/null
+++ b/previews/PR807/assets/plotting.md.n1eGQxVQ.js
@@ -0,0 +1,135 @@
+import{_ as a,c as i,a5 as n,o as t}from"./chunks/framework.Cl7EIXwS.js";const e="/Rasters.jl/previews/PR807/assets/motceku.oyWWmvC1.png",p="/Rasters.jl/previews/PR807/assets/nuekexd.BRgQ_uqM.png",l="/Rasters.jl/previews/PR807/assets/edcrsgk.2YhfOv1i.png",h="/Rasters.jl/previews/PR807/assets/maoowux.CO0mgiME.png",k="/Rasters.jl/previews/PR807/assets/qeatmhs.BsE3Kr89.png",o="/Rasters.jl/previews/PR807/assets/huuwpka.CVxd1wxn.png",d="/Rasters.jl/previews/PR807/assets/ouepwku.BPZOvOxi.png",r="/Rasters.jl/previews/PR807/assets/xqojlxo.Djo3H46f.png",c="/Rasters.jl/previews/PR807/assets/uedqoev.CUldC_pQ.png",g="/Rasters.jl/previews/PR807/assets/eqdfnkz.8o5RdXOt.png",E="/Rasters.jl/previews/PR807/assets/spxvigw.D9tQTCLW.png",u="/Rasters.jl/previews/PR807/assets/ryavibo.H6rBHY04.png",D=JSON.parse('{"title":"","description":"","frontmatter":{},"headers":[],"relativePath":"plotting.md","filePath":"plotting.md","lastUpdated":null}'),y={name:"plotting.md"};function m(F,s,C,v,b,q){return t(),i("div",null,s[0]||(s[0]=[n(`
Plots.jl and Makie.jl are fully supported by Rasters.jl, with recipes for plotting Raster and RasterStack provided. plot will plot a heatmap with axes matching dimension values. If mappedcrs is used, converted values will be shown on axes instead of the underlying crs values. contourf will similarly plot a filled contour plot.
Pixel resolution is limited to allow loading very large files quickly. max_res specifies the maximum pixel resolution to show on the longest axis of the array. It can be set manually to change the resolution (e.g. for large or high-quality plots):
For Makie, plot functions in a similar way. plot will only accept two-dimensional rasters. You can invoke contour, contourf, heatmap, surface or any Makie plotting function which supports surface-like data on a 2D raster.
To obtain tiled plots for 3D rasters and RasterStacks, use the function Rasters.rplot([gridposition], raster; kw_args...). This is an unexported function, since we're not sure how the API will change going forward.
╭──────────────────────────────────────────────────╮
+│ 180×170×24 Raster{Union{Missing, Float32},3} tos │
+├──────────────────────────────────────────────────┴───────────────────── dims ┐
+ ↓ X Mapped{Float64} [1.0, 3.0, …, 357.0, 359.0] ForwardOrdered Regular Intervals{Center},
+ → Y Mapped{Float64} [-79.5, -78.5, …, 88.5, 89.5] ForwardOrdered Regular Intervals{Center},
+ ↗ Ti Sampled{DateTime360Day} [DateTime360Day(2001-01-16T00:00:00), …, DateTime360Day(2002-12-16T00:00:00)] ForwardOrdered Explicit Intervals{Center}
+├──────────────────────────────────────────────────────────────────── metadata ┤
+ Metadata{Rasters.NCDsource} of Dict{String, Any} with 9 entries:
+ "units" => "K"
+ "missing_value" => 1.0f20
+ "original_units" => "degC"
+ "cell_methods" => "time: mean (interval: 30 minutes)"
+ "history" => " At 16:37:23 on 01/11/2005: CMOR altered the data in t…
+ "long_name" => "Sea Surface Temperature"
+ "standard_name" => "sea_surface_temperature"
+ "_FillValue" => 1.0f20
+ "original_name" => "sosstsst"
+├────────────────────────────────────────────────────────────────────── raster ┤
+ extent: Extent(X = (0.0, 360.0), Y = (-80.0, 90.0), Ti = (DateTime360Day(2001-01-01T00:00:00), DateTime360Day(2003-01-01T00:00:00)))
+ missingval: missing
+ crs: EPSG:4326
+ mappedcrs: EPSG:4326
+└──────────────────────────────────────────────────────────────────────────────┘
+[:, :, 1]
+ ⋮ ⋱
Objects with Dimensions other than X and Y will produce multi-pane plots. Here we plot every third month in the first year in one plot:
julia
A[Ti=1:3:12] |> plot
Now plot the ocean temperatures around the Americas in the first month of 2001. Notice we are using lat/lon coordinates and date/time instead of regular indices. The time dimension uses DateTime360Day, so we need to load CFTime.jl to index it with Near.
Now get the mean over the timespan, then save it to disk, and plot it as a filled contour.
Other plot functions and sliced objects that have only one X/Y/Z dimension fall back to generic DimensionalData.jl plotting, which will still correctly label plot axes.
julia
using Statistics
+# Take the mean
+mean_tos = mean(A; dims=Ti)
╭─────────────────────────────────────────────────╮
+│ 180×170×1 Raster{Union{Missing, Float32},3} tos │
+├─────────────────────────────────────────────────┴────────────────────── dims ┐
+ ↓ X Mapped{Float64} [1.0, 3.0, …, 357.0, 359.0] ForwardOrdered Regular Intervals{Center},
+ → Y Mapped{Float64} [-79.5, -78.5, …, 88.5, 89.5] ForwardOrdered Regular Intervals{Center},
+ ↗ Ti Sampled{DateTime360Day} DateTime360Day(2002-01-16T00:00:00):Millisecond(2592000000):DateTime360Day(2002-01-16T00:00:00) ForwardOrdered Explicit Intervals{Center}
+├──────────────────────────────────────────────────────────────────── metadata ┤
+ Metadata{Rasters.NCDsource} of Dict{String, Any} with 9 entries:
+ "units" => "K"
+ "missing_value" => 1.0f20
+ "original_units" => "degC"
+ "cell_methods" => "time: mean (interval: 30 minutes)"
+ "history" => " At 16:37:23 on 01/11/2005: CMOR altered the data in t…
+ "long_name" => "Sea Surface Temperature"
+ "standard_name" => "sea_surface_temperature"
+ "_FillValue" => 1.0f20
+ "original_name" => "sosstsst"
+├────────────────────────────────────────────────────────────────────── raster ┤
+ extent: Extent(X = (0.0, 360.0), Y = (-80.0, 90.0), Ti = (DateTime360Day(2001-01-01T00:00:00), DateTime360Day(2003-01-01T00:00:00)))
+ missingval: missing
+ crs: EPSG:4326
+ mappedcrs: EPSG:4326
+└──────────────────────────────────────────────────────────────────────────────┘
+[:, :, 1]
+ ⋮ ⋱
Plotting recipes in DimensionalData.jl are the fallback for Rasters.jl when the object doesn't have 2 X/Y/Z dimensions, or a non-spatial plot command is used. So (as a random example) we could plot a transect of ocean surface temperature at 20 degree latitude :
Rasters.jl provides a range of other methods that are being added to over time. Where applicable these methods read and write lazily to and from disk-based arrays of common raster file types. These methods also work for entire RasterStacks and RasterSeries using the same syntax.
`,78)]))}const f=a(y,[["render",m]]);export{D as __pageData,f as default};
diff --git a/previews/PR807/assets/plotting.md.n1eGQxVQ.lean.js b/previews/PR807/assets/plotting.md.n1eGQxVQ.lean.js
new file mode 100644
index 00000000..288b22f5
--- /dev/null
+++ b/previews/PR807/assets/plotting.md.n1eGQxVQ.lean.js
@@ -0,0 +1,135 @@
+import{_ as a,c as i,a5 as n,o as t}from"./chunks/framework.Cl7EIXwS.js";const e="/Rasters.jl/previews/PR807/assets/motceku.oyWWmvC1.png",p="/Rasters.jl/previews/PR807/assets/nuekexd.BRgQ_uqM.png",l="/Rasters.jl/previews/PR807/assets/edcrsgk.2YhfOv1i.png",h="/Rasters.jl/previews/PR807/assets/maoowux.CO0mgiME.png",k="/Rasters.jl/previews/PR807/assets/qeatmhs.BsE3Kr89.png",o="/Rasters.jl/previews/PR807/assets/huuwpka.CVxd1wxn.png",d="/Rasters.jl/previews/PR807/assets/ouepwku.BPZOvOxi.png",r="/Rasters.jl/previews/PR807/assets/xqojlxo.Djo3H46f.png",c="/Rasters.jl/previews/PR807/assets/uedqoev.CUldC_pQ.png",g="/Rasters.jl/previews/PR807/assets/eqdfnkz.8o5RdXOt.png",E="/Rasters.jl/previews/PR807/assets/spxvigw.D9tQTCLW.png",u="/Rasters.jl/previews/PR807/assets/ryavibo.H6rBHY04.png",D=JSON.parse('{"title":"","description":"","frontmatter":{},"headers":[],"relativePath":"plotting.md","filePath":"plotting.md","lastUpdated":null}'),y={name:"plotting.md"};function m(F,s,C,v,b,q){return t(),i("div",null,s[0]||(s[0]=[n(`
Plots.jl and Makie.jl are fully supported by Rasters.jl, with recipes for plotting Raster and RasterStack provided. plot will plot a heatmap with axes matching dimension values. If mappedcrs is used, converted values will be shown on axes instead of the underlying crs values. contourf will similarly plot a filled contour plot.
Pixel resolution is limited to allow loading very large files quickly. max_res specifies the maximum pixel resolution to show on the longest axis of the array. It can be set manually to change the resolution (e.g. for large or high-quality plots):
For Makie, plot functions in a similar way. plot will only accept two-dimensional rasters. You can invoke contour, contourf, heatmap, surface or any Makie plotting function which supports surface-like data on a 2D raster.
To obtain tiled plots for 3D rasters and RasterStacks, use the function Rasters.rplot([gridposition], raster; kw_args...). This is an unexported function, since we're not sure how the API will change going forward.
╭──────────────────────────────────────────────────╮
+│ 180×170×24 Raster{Union{Missing, Float32},3} tos │
+├──────────────────────────────────────────────────┴───────────────────── dims ┐
+ ↓ X Mapped{Float64} [1.0, 3.0, …, 357.0, 359.0] ForwardOrdered Regular Intervals{Center},
+ → Y Mapped{Float64} [-79.5, -78.5, …, 88.5, 89.5] ForwardOrdered Regular Intervals{Center},
+ ↗ Ti Sampled{DateTime360Day} [DateTime360Day(2001-01-16T00:00:00), …, DateTime360Day(2002-12-16T00:00:00)] ForwardOrdered Explicit Intervals{Center}
+├──────────────────────────────────────────────────────────────────── metadata ┤
+ Metadata{Rasters.NCDsource} of Dict{String, Any} with 9 entries:
+ "units" => "K"
+ "missing_value" => 1.0f20
+ "original_units" => "degC"
+ "cell_methods" => "time: mean (interval: 30 minutes)"
+ "history" => " At 16:37:23 on 01/11/2005: CMOR altered the data in t…
+ "long_name" => "Sea Surface Temperature"
+ "standard_name" => "sea_surface_temperature"
+ "_FillValue" => 1.0f20
+ "original_name" => "sosstsst"
+├────────────────────────────────────────────────────────────────────── raster ┤
+ extent: Extent(X = (0.0, 360.0), Y = (-80.0, 90.0), Ti = (DateTime360Day(2001-01-01T00:00:00), DateTime360Day(2003-01-01T00:00:00)))
+ missingval: missing
+ crs: EPSG:4326
+ mappedcrs: EPSG:4326
+└──────────────────────────────────────────────────────────────────────────────┘
+[:, :, 1]
+ ⋮ ⋱
Objects with Dimensions other than X and Y will produce multi-pane plots. Here we plot every third month in the first year in one plot:
julia
A[Ti=1:3:12] |> plot
Now plot the ocean temperatures around the Americas in the first month of 2001. Notice we are using lat/lon coordinates and date/time instead of regular indices. The time dimension uses DateTime360Day, so we need to load CFTime.jl to index it with Near.
Now get the mean over the timespan, then save it to disk, and plot it as a filled contour.
Other plot functions and sliced objects that have only one X/Y/Z dimension fall back to generic DimensionalData.jl plotting, which will still correctly label plot axes.
julia
using Statistics
+# Take the mean
+mean_tos = mean(A; dims=Ti)
╭─────────────────────────────────────────────────╮
+│ 180×170×1 Raster{Union{Missing, Float32},3} tos │
+├─────────────────────────────────────────────────┴────────────────────── dims ┐
+ ↓ X Mapped{Float64} [1.0, 3.0, …, 357.0, 359.0] ForwardOrdered Regular Intervals{Center},
+ → Y Mapped{Float64} [-79.5, -78.5, …, 88.5, 89.5] ForwardOrdered Regular Intervals{Center},
+ ↗ Ti Sampled{DateTime360Day} DateTime360Day(2002-01-16T00:00:00):Millisecond(2592000000):DateTime360Day(2002-01-16T00:00:00) ForwardOrdered Explicit Intervals{Center}
+├──────────────────────────────────────────────────────────────────── metadata ┤
+ Metadata{Rasters.NCDsource} of Dict{String, Any} with 9 entries:
+ "units" => "K"
+ "missing_value" => 1.0f20
+ "original_units" => "degC"
+ "cell_methods" => "time: mean (interval: 30 minutes)"
+ "history" => " At 16:37:23 on 01/11/2005: CMOR altered the data in t…
+ "long_name" => "Sea Surface Temperature"
+ "standard_name" => "sea_surface_temperature"
+ "_FillValue" => 1.0f20
+ "original_name" => "sosstsst"
+├────────────────────────────────────────────────────────────────────── raster ┤
+ extent: Extent(X = (0.0, 360.0), Y = (-80.0, 90.0), Ti = (DateTime360Day(2001-01-01T00:00:00), DateTime360Day(2003-01-01T00:00:00)))
+ missingval: missing
+ crs: EPSG:4326
+ mappedcrs: EPSG:4326
+└──────────────────────────────────────────────────────────────────────────────┘
+[:, :, 1]
+ ⋮ ⋱
Plotting recipes in DimensionalData.jl are the fallback for Rasters.jl when the object doesn't have 2 X/Y/Z dimensions, or a non-spatial plot command is used. So (as a random example) we could plot a transect of ocean surface temperature at 20 degree latitude :
Rasters.jl provides a range of other methods that are being added to over time. Where applicable these methods read and write lazily to and from disk-based arrays of common raster file types. These methods also work for entire RasterStacks and RasterSeries using the same syntax.
Rasters.jl uses a number of backends to load raster data. Raster, RasterStack and RasterSeries will detect which backend to use for you, automatically.
R GRD files can be loaded natively, using Julias MMap - which means they are very fast, but are not compressed. They are always 3 dimensional, and have Y, X and Band dimensions.
NetCDF .nc files are loaded using NCDatasets.jl. Layers from files can be loaded as Raster("filename.nc"; name=:layername). Without name the first layer is used. RasterStack("filename.nc") will use all netcdf variables in the file that are not dimensions as layers.
NetCDF layers can have arbitrary dimensions. Known, common dimension names are converted to X, YZ, and Ti, otherwise Dim{:layername} is used. Layers in the same file may also have different dimensions.
NetCDF files still have issues loading directly from disk for some operations. Using read(ncstack) may help.
All files GDAL can access, such as .tiff and .asc files, can be loaded, using ArchGDAL.jl. These are generally best loaded as Raster("filename.tif"), but can be loaded as RasterStack("filename.tif"; layersfrom=Band), taking layers from the Band dimension, which is also the default.
The Soil Moisture Active-Passive dataset provides global layers of soil moisture, temperature and other related data, in a custom HDF5 format. Layers are always 2 dimensional, with Y and X dimensions.
These can be loaded as multi-layered RasterStack("filename.h5"). Individual layers can be loaded as Raster("filename.h5"; name=:layername), without name the first layer is used.
julia
using Rasters
Missing docstring.
Missing docstring for smapseries. Check Documenter's build log for details.
Files can be written to disk in all formats other than SMAP HDF5 using write("filename.ext", A). See the docs for write. They can (with some caveats) be written to different formats than they were loaded in as, providing file-type conversion for spatial data.
Some metadata may be lost in formats that store little metadata, or where metadata conversion has not been completely implemented.
RasterDataSources.jl standardises the download of common raster data sources, with a focus on datasets used in ecology and the environmental sciences. RasterDataSources.jl is tightly integrated into Rasters.jl, so that datsets and keywords can be used directly to download and load data as a Raster, RasterStack, or RasterSeries.
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.
julia
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
More options are available, like Near, Contains and Where.
Dimensions
Rasters uses X, Y, and Z dimensions from DimensionalData to represent spatial directions like longitude, latitude and the vertical dimension, and subset data with them. Ti is used for time, and Band represent bands. Other dimensions can have arbitrary names, but will be treated generically. See DimensionalData for more details on how they work.
Lookup Arrays
These specify properties of the index associated with e.g. the X and Y dimension. Rasters.jl defines additional lookup arrays: Projected to handle dimensions with projections, and Mapped where the projection is mapped to another projection like EPSG(4326). Mapped is largely designed to handle NetCDF dimensions, especially with Explicit spans.
Regular getindex (e.g. A[1:100, :]) and view work on all objects just as with an Array. view is always lazy, and reads from disk are deferred until getindex is used. DimensionalData.jlDimensions and Selectors are the other way to subset an object, making use of the objects index to find values at e.g. certain X/Y coordinates. The available selectors are listed here:
Selectors
Description
At(x)
get the index exactly matching the passed in value(s).
Near(x)
get the closest index to the passed in value(s).
Where(f::Function)
filter the array axis by a function of the dimension index values.
a..b/Between(a, b)
get all indices between two values, excluding the high value.
Contains(x)
get indices where the value x falls within an interval.
Info
Use the .. selector to take a view of madagascar:
julia
using Rasters, RasterDataSources
+const RS = Rasters
+using CairoMakie
+CairoMakie.activate!()
+
+A = Raster(WorldClim{BioClim}, 5)
+madagascar = view(A, X(43.25 .. 50.48), Y(-25.61 .. -12.04))
+# Note the space between .. -12
+Makie.plot(madagascar)
+
+
+
+
\ No newline at end of file
diff --git a/previews/PR807/hashmap.json b/previews/PR807/hashmap.json
new file mode 100644
index 00000000..cd80ee2c
--- /dev/null
+++ b/previews/PR807/hashmap.json
@@ -0,0 +1 @@
+{"api.md":"OKz0SH8-","array_operations.md":"CRu0fcQR","data_sources.md":"Dcl5DOui","gbif_wflow.md":"B-xemxXG","get_started.md":"CRLYYECF","index.md":"CYSmGW74","methods.md":"tEJSj4eC","plot_makie.md":"DYZzsleT","plotting.md":"n1eGQxVQ"}
diff --git a/previews/PR807/index.html b/previews/PR807/index.html
new file mode 100644
index 00000000..2ce2c359
--- /dev/null
+++ b/previews/PR807/index.html
@@ -0,0 +1,25 @@
+
+
+
+
+
+ Rasters.jl
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
These currently include raster arrays like GeoTIF and NetCDF, R grd files, multi-layered stacks, and multi-file series of arrays and stacks.
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.
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:
julia
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: not yet finished.
Other functionality in extensions:
Raster data downloads, like Worldclim{Climate}: using RasterDataSources
Makie plots: using Makie
Coordinate transformations for gdal rasters: using CoordinateTransformations
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, that's great! but its not always possible).
To make an issue we can fix quickly (or at all) there are three key steps:
Include the file in an accessible place on web without authentication 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.
Add a minimum working example to the issue template that first downloads the file, then runs the function that triggers the bug.
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.
+
+
+
+
\ No newline at end of file
diff --git a/previews/PR807/logo.png b/previews/PR807/logo.png
new file mode 100644
index 00000000..4b53fffb
Binary files /dev/null and b/previews/PR807/logo.png differ
diff --git a/previews/PR807/methods.html b/previews/PR807/methods.html
new file mode 100644
index 00000000..9fab7f85
--- /dev/null
+++ b/previews/PR807/methods.html
@@ -0,0 +1,25 @@
+
+
+
+
+
+ Rasters.jl
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Note that most regular Julia methods, such as replace, work as for a standard Array. These additional methods are commonly required in GIS applications.
Plotting in Makie works somewhat differently than Plots, since the recipe system is different. You can pass a 2-D raster to any surface-like function (heatmap, contour, contourf, or even surface for a 3D plot) with ease.
fig, ax, _ = plot(A)
+contour(fig[1, 2], A)
+ax = Axis(fig[2, 1]; aspect = DataAspect())
+contourf!(ax, A)
+surface(fig[2, 2], A) # even a 3D plot work!
+fig
Rasters.rplot should support Observable input out of the box, but the dimensions of that input must remain the same - i.e., the element names of a RasterStack must remain the same.
julia
Makie.set_theme!(Rasters.theme_rasters())
+# `stack` is the WorldClim climate data for January
+stack_obs = Observable(stack)
+fig = Rasters.rplot(stack_obs;
+ Colorbar=(; height=Relative(0.75), width=5)
+)
+record(fig, "rplot.mp4", 1:12; framerate = 3) do i
+ stack_obs[] = RasterStack(WorldClim{Climate}; month = i)
+end
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 amount 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.
Plots.jl and Makie.jl are fully supported by Rasters.jl, with recipes for plotting Raster and RasterStack provided. plot will plot a heatmap with axes matching dimension values. If mappedcrs is used, converted values will be shown on axes instead of the underlying crs values. contourf will similarly plot a filled contour plot.
Pixel resolution is limited to allow loading very large files quickly. max_res specifies the maximum pixel resolution to show on the longest axis of the array. It can be set manually to change the resolution (e.g. for large or high-quality plots):
For Makie, plot functions in a similar way. plot will only accept two-dimensional rasters. You can invoke contour, contourf, heatmap, surface or any Makie plotting function which supports surface-like data on a 2D raster.
To obtain tiled plots for 3D rasters and RasterStacks, use the function Rasters.rplot([gridposition], raster; kw_args...). This is an unexported function, since we're not sure how the API will change going forward.
╭──────────────────────────────────────────────────╮
+│ 180×170×24 Raster{Union{Missing, Float32},3} tos │
+├──────────────────────────────────────────────────┴───────────────────── dims ┐
+ ↓ X Mapped{Float64} [1.0, 3.0, …, 357.0, 359.0] ForwardOrdered Regular Intervals{Center},
+ → Y Mapped{Float64} [-79.5, -78.5, …, 88.5, 89.5] ForwardOrdered Regular Intervals{Center},
+ ↗ Ti Sampled{DateTime360Day} [DateTime360Day(2001-01-16T00:00:00), …, DateTime360Day(2002-12-16T00:00:00)] ForwardOrdered Explicit Intervals{Center}
+├──────────────────────────────────────────────────────────────────── metadata ┤
+ Metadata{Rasters.NCDsource} of Dict{String, Any} with 9 entries:
+ "units" => "K"
+ "missing_value" => 1.0f20
+ "original_units" => "degC"
+ "cell_methods" => "time: mean (interval: 30 minutes)"
+ "history" => " At 16:37:23 on 01/11/2005: CMOR altered the data in t…
+ "long_name" => "Sea Surface Temperature"
+ "standard_name" => "sea_surface_temperature"
+ "_FillValue" => 1.0f20
+ "original_name" => "sosstsst"
+├────────────────────────────────────────────────────────────────────── raster ┤
+ extent: Extent(X = (0.0, 360.0), Y = (-80.0, 90.0), Ti = (DateTime360Day(2001-01-01T00:00:00), DateTime360Day(2003-01-01T00:00:00)))
+ missingval: missing
+ crs: EPSG:4326
+ mappedcrs: EPSG:4326
+└──────────────────────────────────────────────────────────────────────────────┘
+[:, :, 1]
+ ⋮ ⋱
Objects with Dimensions other than X and Y will produce multi-pane plots. Here we plot every third month in the first year in one plot:
julia
A[Ti=1:3:12] |> plot
Now plot the ocean temperatures around the Americas in the first month of 2001. Notice we are using lat/lon coordinates and date/time instead of regular indices. The time dimension uses DateTime360Day, so we need to load CFTime.jl to index it with Near.
Now get the mean over the timespan, then save it to disk, and plot it as a filled contour.
Other plot functions and sliced objects that have only one X/Y/Z dimension fall back to generic DimensionalData.jl plotting, which will still correctly label plot axes.
julia
using Statistics
+# Take the mean
+mean_tos = mean(A; dims=Ti)
╭─────────────────────────────────────────────────╮
+│ 180×170×1 Raster{Union{Missing, Float32},3} tos │
+├─────────────────────────────────────────────────┴────────────────────── dims ┐
+ ↓ X Mapped{Float64} [1.0, 3.0, …, 357.0, 359.0] ForwardOrdered Regular Intervals{Center},
+ → Y Mapped{Float64} [-79.5, -78.5, …, 88.5, 89.5] ForwardOrdered Regular Intervals{Center},
+ ↗ Ti Sampled{DateTime360Day} DateTime360Day(2002-01-16T00:00:00):Millisecond(2592000000):DateTime360Day(2002-01-16T00:00:00) ForwardOrdered Explicit Intervals{Center}
+├──────────────────────────────────────────────────────────────────── metadata ┤
+ Metadata{Rasters.NCDsource} of Dict{String, Any} with 9 entries:
+ "units" => "K"
+ "missing_value" => 1.0f20
+ "original_units" => "degC"
+ "cell_methods" => "time: mean (interval: 30 minutes)"
+ "history" => " At 16:37:23 on 01/11/2005: CMOR altered the data in t…
+ "long_name" => "Sea Surface Temperature"
+ "standard_name" => "sea_surface_temperature"
+ "_FillValue" => 1.0f20
+ "original_name" => "sosstsst"
+├────────────────────────────────────────────────────────────────────── raster ┤
+ extent: Extent(X = (0.0, 360.0), Y = (-80.0, 90.0), Ti = (DateTime360Day(2001-01-01T00:00:00), DateTime360Day(2003-01-01T00:00:00)))
+ missingval: missing
+ crs: EPSG:4326
+ mappedcrs: EPSG:4326
+└──────────────────────────────────────────────────────────────────────────────┘
+[:, :, 1]
+ ⋮ ⋱
Plotting recipes in DimensionalData.jl are the fallback for Rasters.jl when the object doesn't have 2 X/Y/Z dimensions, or a non-spatial plot command is used. So (as a random example) we could plot a transect of ocean surface temperature at 20 degree latitude :
Rasters.jl provides a range of other methods that are being added to over time. Where applicable these methods read and write lazily to and from disk-based arrays of common raster file types. These methods also work for entire RasterStacks and RasterSeries using the same syntax.
+
+
+
+
\ No newline at end of file
diff --git a/previews/PR807/siteinfo.js b/previews/PR807/siteinfo.js
new file mode 100644
index 00000000..ea3998f2
--- /dev/null
+++ b/previews/PR807/siteinfo.js
@@ -0,0 +1 @@
+var DOCUMENTER_CURRENT_VERSION = "previews/PR807";
