Rasters.jl

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.

Mean over the time dimension:

julia

using Rasters, Statistics, RasterDataSources

A = Raster(WorldClim{BioClim}, 5)

╭──────────────────────────────────╮
│ 2160×1080 Raster{Float32,2} bio5 │
├──────────────────────────────────┴───────────────────────────────────── dims ┐
  ↓ X Projected{Float64} LinRange{Float64}(-180.0, 179.83333333333331, 2160) ForwardOrdered Regular Intervals{Start},
  → Y Projected{Float64} LinRange{Float64}(89.83333333333333, -90.0, 1080) ReverseOrdered Regular Intervals{Start}
├──────────────────────────────────────────────────────────────────── metadata ┤
  Metadata{Rasters.GDALsource} of Dict{String, Any} with 4 entries:
  "units"    => ""
  "offset"   => 0.0
  "filepath" => "./WorldClim/BioClim/wc2.1_10m_bio_5.tif"
  "scale"    => 1.0
├────────────────────────────────────────────────────────────────────── raster ┤
  extent: Extent(X = (-180.0, 179.99999999999997), Y = (-90.0, 90.0))
  missingval: -3.4f38
  crs: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
└──────────────────────────────────────────────────────────────────────────────┘
    ↓ →    89.8333  89.6667  89.5     89.3333  …  -89.6667  -89.8333  -90.0
 -180.0    -3.4f38  -3.4f38  -3.4f38  -3.4f38     -15.399   -13.805   -14.046
    ⋮                                          ⋱                        ⋮
  179.833  -3.4f38  -3.4f38  -3.4f38  -3.4f38  …  -17.1478  -15.4243  -15.701

Then we do the mean over the X dimension

julia

mean(A, dims=X) # Ti if time were available would also be possible

╭───────────────────────────────╮
│ 1×1080 Raster{Float32,2} bio5 │
├───────────────────────────────┴──────────────────────────────────────── dims ┐
  ↓ X Projected{Float64} -180.0:360.0:-180.0 ForwardOrdered Regular Intervals{Start},
  → Y Projected{Float64} LinRange{Float64}(89.83333333333333, -90.0, 1080) ReverseOrdered Regular Intervals{Start}
├──────────────────────────────────────────────────────────────────── metadata ┤
  Metadata{Rasters.GDALsource} of Dict{String, Any} with 4 entries:
  "units"    => ""
  "offset"   => 0.0
  "filepath" => "./WorldClim/BioClim/wc2.1_10m_bio_5.tif"
  "scale"    => 1.0
├────────────────────────────────────────────────────────────────────── raster ┤
  extent: Extent(X = (-180.0, 180.0), Y = (-90.0, 90.0))
  missingval: -3.4f38
  crs: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
└──────────────────────────────────────────────────────────────────────────────┘
    ↓ →   89.8333   89.6667   89.5   89.3333  …  -89.6667  -89.8333  -90.0
 -180.0  -Inf      -Inf      -Inf   -Inf         -23.0768  -22.9373  -22.0094

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.

Broadcasting

For a disk-based array A, this will only be applied when indexing occurs or when we read the array.

julia

A .*= 2

╭──────────────────────────────────╮
│ 2160×1080 Raster{Float32,2} bio5 │
├──────────────────────────────────┴───────────────────────────────────── dims ┐
  ↓ X Projected{Float64} LinRange{Float64}(-180.0, 179.83333333333331, 2160) ForwardOrdered Regular Intervals{Start},
  → Y Projected{Float64} LinRange{Float64}(89.83333333333333, -90.0, 1080) ReverseOrdered Regular Intervals{Start}
├──────────────────────────────────────────────────────────────────── metadata ┤
  Metadata{Rasters.GDALsource} of Dict{String, Any} with 4 entries:
  "units"    => ""
  "offset"   => 0.0
  "filepath" => "./WorldClim/BioClim/wc2.1_10m_bio_5.tif"
  "scale"    => 1.0
├────────────────────────────────────────────────────────────────────── raster ┤
  extent: Extent(X = (-180.0, 179.99999999999997), Y = (-90.0, 90.0))
  missingval: -3.4f38
  crs: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
└──────────────────────────────────────────────────────────────────────────────┘
    ↓ →     89.8333   89.6667   89.5  …  -89.6667  -89.8333  -90.0
 -180.0    -Inf      -Inf      -Inf      -30.798   -27.61    -28.092
    ⋮                                 ⋱                        ⋮
  179.833  -Inf      -Inf      -Inf   …  -34.2955  -30.8485  -31.402

To broadcast directly to disk, we need to open the file in write mode:

julia

open(Raster(filename); write=true) do O
   O .*= 2
end

To broadcast over a RasterStack use map, which applies a function to the raster layers of the stack.

julia

newstack = map(stack) do raster
   raster .* 2
end

Modifying object properties

rebuild can be used to modify the fields of an object, generating a new object (but possibly holding the same arrays or files).

If you know that a file had an incorrectly specified missing value, you can do:

rebuild

julia

A = Raster(WorldClim{BioClim}, 5)
rebuild(A; missingval=-9999)

╭──────────────────────────────────╮
│ 2160×1080 Raster{Float32,2} bio5 │
├──────────────────────────────────┴───────────────────────────────────── dims ┐
  ↓ X Projected{Float64} LinRange{Float64}(-180.0, 179.83333333333331, 2160) ForwardOrdered Regular Intervals{Start},
  → Y Projected{Float64} LinRange{Float64}(89.83333333333333, -90.0, 1080) ReverseOrdered Regular Intervals{Start}
├──────────────────────────────────────────────────────────────────── metadata ┤
  Metadata{Rasters.GDALsource} of Dict{String, Any} with 4 entries:
  "units"    => ""
  "offset"   => 0.0
  "filepath" => "./WorldClim/BioClim/wc2.1_10m_bio_5.tif"
  "scale"    => 1.0
├────────────────────────────────────────────────────────────────────── raster ┤
  extent: Extent(X = (-180.0, 179.99999999999997), Y = (-90.0, 90.0))
  missingval: -9999.0f0
  crs: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
└──────────────────────────────────────────────────────────────────────────────┘
    ↓ →    89.8333  89.6667  89.5     89.3333  …  -89.6667  -89.8333  -90.0
 -180.0    -3.4f38  -3.4f38  -3.4f38  -3.4f38     -15.399   -13.805   -14.046
    ⋮                                          ⋱                        ⋮
  179.833  -3.4f38  -3.4f38  -3.4f38  -3.4f38  …  -17.1478  -15.4243  -15.701

(replace_missing will actually replace the current values).

Or if you need to change the name of the layer:

Then use rebuild as

julia

B = rebuild(A; name=:temperature)
B.name

:temperature

replace_missing

julia

replace_missing(A, missingval=-9999)

╭──────────────────────────────────╮
│ 2160×1080 Raster{Float32,2} bio5 │
├──────────────────────────────────┴───────────────────────────────────── dims ┐
  ↓ X Projected{Float64} LinRange{Float64}(-180.0, 179.83333333333331, 2160) ForwardOrdered Regular Intervals{Start},
  → Y Projected{Float64} LinRange{Float64}(89.83333333333333, -90.0, 1080) ReverseOrdered Regular Intervals{Start}
├──────────────────────────────────────────────────────────────────── metadata ┤
  Metadata{Rasters.GDALsource} of Dict{String, Any} with 4 entries:
  "units"    => ""
  "offset"   => 0.0
  "filepath" => "./WorldClim/BioClim/wc2.1_10m_bio_5.tif"
  "scale"    => 1.0
├────────────────────────────────────────────────────────────────────── raster ┤
  extent: Extent(X = (-180.0, 179.99999999999997), Y = (-90.0, 90.0))
  missingval: -9999.0f0
  crs: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
└──────────────────────────────────────────────────────────────────────────────┘
    ↓ →       89.8333     89.6667     89.5  …  -89.6667  -89.8333  -90.0
 -180.0    -9999.0     -9999.0     -9999.0     -15.399   -13.805   -14.046
    ⋮                                       ⋱                        ⋮
  179.833  -9999.0     -9999.0     -9999.0  …  -17.1478  -15.4243  -15.701

set

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:

julia

using Rasters: Points
set(A, X => Points)

╭──────────────────────────────────╮
│ 2160×1080 Raster{Float32,2} bio5 │
├──────────────────────────────────┴───────────────────────────────────── dims ┐
  ↓ X Projected{Float64} LinRange{Float64}(-180.0, 179.83333333333331, 2160) ForwardOrdered Regular Points,
  → Y Projected{Float64} LinRange{Float64}(89.83333333333333, -90.0, 1080) ReverseOrdered Regular Intervals{Start}
├──────────────────────────────────────────────────────────────────── metadata ┤
  Metadata{Rasters.GDALsource} of Dict{String, Any} with 4 entries:
  "units"    => ""
  "offset"   => 0.0
  "filepath" => "./WorldClim/BioClim/wc2.1_10m_bio_5.tif"
  "scale"    => 1.0
├────────────────────────────────────────────────────────────────────── raster ┤
  extent: Extent(X = (-180.0, 179.83333333333331), Y = (-90.0, 90.0))
  missingval: -3.4f38
  crs: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
└──────────────────────────────────────────────────────────────────────────────┘
    ↓ →    89.8333  89.6667  89.5     89.3333  …  -89.6667  -89.8333  -90.0
 -180.0    -3.4f38  -3.4f38  -3.4f38  -3.4f38     -15.399   -13.805   -14.046
    ⋮                                          ⋱                        ⋮
  179.833  -3.4f38  -3.4f38  -3.4f38  -3.4f38  …  -17.1478  -15.4243  -15.701

We can also reassign dimensions, here X becomes Z:

julia

set(A, X => Z)

╭──────────────────────────────────╮
│ 2160×1080 Raster{Float32,2} bio5 │
├──────────────────────────────────┴───────────────────────────────────── dims ┐
  ↓ Z Projected{Float64} LinRange{Float64}(-180.0, 179.83333333333331, 2160) ForwardOrdered Regular Intervals{Start},
  → Y Projected{Float64} LinRange{Float64}(89.83333333333333, -90.0, 1080) ReverseOrdered Regular Intervals{Start}
├──────────────────────────────────────────────────────────────────── metadata ┤
  Metadata{Rasters.GDALsource} of Dict{String, Any} with 4 entries:
  "units"    => ""
  "offset"   => 0.0
  "filepath" => "./WorldClim/BioClim/wc2.1_10m_bio_5.tif"
  "scale"    => 1.0
├────────────────────────────────────────────────────────────────────── raster ┤
  extent: Extent(Z = (-180.0, 179.99999999999997), Y = (-90.0, 90.0))
  missingval: -3.4f38
  crs: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
└──────────────────────────────────────────────────────────────────────────────┘
    ↓ →    89.8333  89.6667  89.5     89.3333  …  -89.6667  -89.8333  -90.0
 -180.0    -3.4f38  -3.4f38  -3.4f38  -3.4f38     -15.399   -13.805   -14.046
    ⋮                                          ⋱                        ⋮
  179.833  -3.4f38  -3.4f38  -3.4f38  -3.4f38  …  -17.1478  -15.4243  -15.701

setcrs(A, crs) and setmappedcrs(A, crs) will set the crs value/s of an object to any GeoFormat from GeoFormatTypes.jl.

Mean over the time dimension: ​

Broadcasting ​

Modifying object properties ​

rebuild ​

replace_missing ​

set ​

Mean over the time dimension:

Broadcasting

Modifying object properties

rebuild

replace_missing

set