This vignette explains the data model of `stars`

objects, with illustrations.

With a very simple file created from a \(4 \times 5\) matrix

```
suppressPackageStartupMessages(library(stars))
m = matrix(1:20, nrow = 5, ncol = 4)
dim(m) = c(x = 5, y = 4) # named dim
(s = st_as_stars(m))
## stars object with 2 dimensions and 1 attribute
## attribute(s):
## A1
## Min. : 1.00
## 1st Qu.: 5.75
## Median :10.50
## Mean :10.50
## 3rd Qu.:15.25
## Max. :20.00
## dimension(s):
## from to offset delta refsys point values
## x 1 5 0 1 NA NA NULL
## y 1 4 0 1 NA NA NULL
```

we see that

- the rows (5) are mapped to the first dimension, the x-coordinate
- the columns (4) are mapped to the second dimension, the y-coordinate
- the
`from`

and`to`

fields of the dimensions are redundant, as they also are present in the array dimension:

```
dim(s[[1]])
## x y
## 5 4
```

- offset and delta specify how increasing row/column index maps to x and y coordinate values

When we plot this object, using the `image`

method for `stars`

objects, we see

`image(s, text_values = TRUE, axes = TRUE)`

Where it becomes clear that \((0,0)\) is the origin of the grid (grid corner), and \(1\) the coordinate value increase from one index (row, col) to the next. It means that consecutive matrix columns represent grid lines, going from south to north. Grids defined this way are **regular**: grid cell size is constant everywhere.

Most grid data comes with y coordinates (grid rows) going from North to South (top to bottom); this is established using a negative `delta`

. We see that the grid origing \((0,0)\) did not change:

```
attr(s, "dimensions")$y$delta = -1
attr(s, "dimensions")$x$geotransform[6] = -1
attr(s, "dimensions")$y$geotransform[6] = -1
image(s, text_values = TRUE, axes = TRUE)
```

An example is the GeoTIFF carried in the package, which has a negative `delta`

for the `y`

-coordinate:

```
tif = system.file("tif/L7_ETMs.tif", package = "stars")
st_dimensions(read_stars(tif))["y"]
## from to offset delta refsys point values
## y 1 352 9120761 -28.5 +proj=utm +zone=25 +south... FALSE NULL
```

`x`

and `y`

coordinate grids carry a `geotransform`

field, used to compute \(x\) and \(y\) coordinates from grid index \(i\) and \(j\):

\[x = gt_1 + (i-1) gt_2 + (j-1) gt_3\]

\[y = gt_4 + (i-1) gt_5 + (j-1) gt_6\] We can rotate grids by setting \(gt_3\) and \(gt_5\) to a constant, non-zero value:

```
attr(s, "dimensions")$x$geotransform[c(3,5)] = 0.1
attr(s, "dimensions")$y$geotransform[c(3,5)] = 0.1
# FIXME: use image()
plot(st_as_sf(s, as_points = FALSE), axes = TRUE, nbreaks = 20)
```

The rotation angle, in degrees, is

```
atan2(0.1, 1) * 180 / pi
## [1] 5.710593
```

Sheared grids are obtained when the two rotation coefficients, \(gt_3\) and \(gt_5\), are unequal:

```
attr(s, "dimensions")$x$geotransform[c(3,5)] = c(0.1, 0.2)
attr(s, "dimensions")$y$geotransform[c(3,5)] = c(0.1, 0.2)
plot(st_as_sf(s, as_points = FALSE), axes = TRUE, nbreaks = 20)
```