# "Raster" syntax in GeoStats.jl

**URL:** <https://discourse.julialang.org/t/raster-syntax-in-geostats-jl/119214>\
**Category:** Geo\
**Tags:** indexing, array, geo, rasters\
**Created:** [September 9, 2024, 11:07am UTC](https://discourse.julialang.org/t/raster-syntax-in-geostats-jl/119214 "2024-09-09T11:07:09Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [September 9, 2024, 11:07am UTC](https://discourse.julialang.org/t/raster-syntax-in-geostats-jl/119214/1 "2024-09-09T11:07:09Z")

</div>

Given that most people in the `geo` community are not aware of the “raster” syntax in GeoStats.jl, I collect a few examples here to facilitate the transition. This post may become part of the official project documentation in the future.

# The “raster” model

A “raster” is nothing more than a Julia array that is georeferenced over a geospatial grid. The “raster” model is popular in the scientific community for two main reasons:

1. it provides convenient syntax to access grid elements within “rectangular” regions,  
including syntax to name array dimensions as “X”, “Y”, … and extract rectangular regions based on coordinate values.

2. it can rely on DiskArrays.jl and similar packages for loading large datasets lazily, without consuming all the available RAM.

In the following sections I illustrate how `GeoTable`s over `Grid`s share these benefits.

## 1. Convenient syntax

### 1.1 n-dimensional array syntax

Any `GeoTable` over a `Grid` domain supports the n-dimensional array syntax in the **row selector**.

Consider the following “raster” from the NaturalEarth dataset:

```julia
using GeoStats
using GeoArtifacts

import GLMakie as Mke

raster = NaturalEarth.naturalearth1("water") |> Upscale(10, 5)

raster |> viewer

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/5/b/5b1b52b321eb410725593c4f077969eba19bbfcc.jpeg)

We can check the size of the underlying grid with:

```julia
julia> size(raster.geometry)
(1620, 1620)

```

This dataset is stored with `LatLon` coordinates. For convenience, we always store the “horizontal” coordinate along the first axis of the grid (longitude in this case). We can slice the Earth as follows:

```julia
raster[(1:800, :), :] |> viewer

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/e/f/ef157dc8568fc3011518de84a572eaf06b3cf8c6.jpeg)

```julia
raster[(800:1620, :), :] |> viewer

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/8/e/8ead5b65bac7a160cc0d6e6d8a69b3f862a3b578.jpeg)

```julia
raster[(:, 1:800), :] |> viewer

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/5/f/5f1209d09c816a59fe63ca945160b74c1748f6aa.jpeg)

```julia
raster[(:, 500:800), :] |> viewer

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/2/5/258f6d30e5c568460434039ca65583c1f02642f9.jpeg)

We can also slice after 2D projections, which is more common in the literature:

```julia
projec = raster |> Proj(Robinson)

projec |> viewer

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/c/a/ca77a785e8923c9fba8d14ac4104843ef8b992df.jpeg)

```julia
projec[(300:800, 600:1400), :] |> viewer

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/3/a/3ad4a0d17e47680a68d349245b53f974bff5a414.jpeg)

Throughout all these examples we used the `:` **column selector** to retain all columns of the geotable. We can retain a subset of the columns with usual dataframe syntax:

```julia
raster[(1:800, 500:800), ["TEMPERATURE", "PRESSURE"]]

```

### 1.2 Slice transform syntax

To select a subset of the dataset from actual coordinate ranges, we can use the `Slice` transform. To slice the `x` and `y` coordinates of the projected dataset with values in kilometers, we do:

```julia
using Unitful: km

projec |> Slice(x=(0km, 10000km), y=(0km, 5000km)) |> viewer

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/8/d/8d198de099f9b8e1bc42afcb3e6f98f740160620.jpeg)

Notice that the **result is no longer a “raster”** because the `Robinson` projection deforms the graticule, and we can’t simply rely on the sorted “x” and “y” coordinates of the _first_ row and column to find indices across the other rows and columns. The following example makes this clearer:

```julia
projec |> Slice(x=(0km, 20000km)) |> viewer

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/d/d/ddf185c93e27d071c6849989edfe40651bdd0e44.jpeg)

The “first row” of the result has more “columns” than the row of the Equator. This result **cannot** be represented with the “raster” model. It **can** , however, be represented with the `GeoTable` model over a domain view.

The names `x` and `y` can be replaced by `lon` and `lat` if the `crs` of the domain is geographic. Below are examples of subgrids over the ellipsoid using lat/lon ranges of known places:

```julia
using Unitful: °

# Brazil
raster |> Slice(lat=(-60°, 20°), lon=(-100°, -20°)) |> viewer

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/c/4/c403844919c5e39325e01a20c7acbf3a5e118d54.jpeg)

```julia
# south pole
raster |> Slice(lat=(-90°, -60°)) |> viewer

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/e/3/e330cf75d796b690d55b815c5838c64997d2afaa.jpeg)

## 2. Lazy loading

The `Grid` itself is lazy, it doesn’t allocate:

```julia
julia> Base.summarysize(raster.geometry)
192

```

Remember that a n-dimensional array is simply a flat memory buffer with an additional tuple containing the number of elements along each dimension. In Julia, the `vec` and `reshape` operations are lazy:

```julia
julia> x = rand(1000, 1000);

julia> @allocated vec(x)
80

julia> @allocated reshape(vec(x), 1000, 1000)
176

```

If you have a DiskArray or any other array type that satisfies your memory requirements, simply `georef` it over a `Grid` to obtain a “lazy” `GeoTable`:

```julia
raster = georef((var1 = vec(array1), var2 = vec(array2)), grid)

```

If you load your data with GeoIO.jl, it will take care of these low-level memory details for you.