---
title: "GDAL Interoperability"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{GDAL Interoperability}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r setup}
library(vaster)
```

## Introduction

GDAL (Geospatial Data Abstraction Library) is the foundation of most geospatial
software. It uses specific conventions for representing raster grids that differ
from R's typical approach.

The **vaster** package provides functions to convert between vaster's
`dimension`/`extent` representation and GDAL's conventions, enabling seamless
integration with GDAL-based tools.

## GDAL GeoTransform

GDAL represents raster georeferencing using a 6-element *GeoTransform* array:

```
GT[0]: x-coordinate of the upper-left corner
GT[1]: pixel width (x resolution)
GT[2]: row rotation (typically 0)
GT[3]: y-coordinate of the upper-left corner
GT[4]: column rotation (typically 0)
GT[5]: pixel height (negative for north-up images)
```

For a typical north-up raster, only elements 0, 1, 3, and 5 are non-zero.

### Converting to GeoTransform

The `geo_transform0()` function converts vaster's dimension/extent to a GDAL
GeoTransform:

```{r}
dimension <- c(10, 5)
extent <- c(0, 100, 0, 50)

gt <- geo_transform0(dimension, extent)
gt
```

The GeoTransform tells us:

- Upper-left corner is at `(0, 50)` — note y is at the top
- Pixel width is `10` units
- Pixel height is `-10` units (negative because y decreases downward)

### Converting from GeoTransform

Given a GeoTransform and dimensions, `gt_dim_to_extent()` recovers the extent:

```{r}
gt <- c(0, 10, 0, 50, 0, -10)
dimension <- c(10, 5)

gt_dim_to_extent(gt, dimension)
```

And `extent_dim_to_gt()` does the reverse:

```{r}
extent_dim_to_gt(c(0, 100, 0, 50), c(10, 5))
```

## World Files

World files (`.tfw`, `.pgw`, `.jgw`, etc.) are a simpler georeferencing format
used by many applications. They contain 6 lines:

```
Line 1: pixel width
Line 2: row rotation
Line 3: column rotation
Line 4: pixel height (negative for north-up)
Line 5: x-coordinate of centre of upper-left pixel
Line 6: y-coordinate of centre of upper-left pixel
```

### Converting to World File Format

```{r}
dimension <- c(10, 5)
extent <- c(0, 100, 0, 50)

wf <- geo_world0(dimension, extent)
wf
```

Note the difference from GeoTransform: the world file specifies the **centre**
of the upper-left pixel, not its corner.

## GDAL RasterIO

GDAL's `RasterIO` function is the core mechanism for reading/writing raster
data. It uses a window specification:

```
(xoff, yoff, xsize, ysize)
```

Where:
- `xoff`: column offset (0-based, from left)
- `yoff`: row offset (0-based, from top)
- `xsize`: number of columns to read
- `ysize`: number of rows to read

### Converting Extent to RasterIO Window

The `rasterio0()` function converts an extent to RasterIO parameters:

```{r}
## Full grid
dimension <- c(360, 180)
extent <- c(-180, 180, -90, 90)

## Subregion to read (Australia)
subextent <- c(110, 155, -45, -10)

rio <- rasterio0(dimension, extent, subextent)
rio
```

This tells you:
- Start at column `r rio[1]` (0-based offset)
- Start at row `r rio[2]` (0-based offset)
- Read `r rio[3]` columns
- Read `r rio[4]` rows

### Using with vapour or stars

These parameters can be passed directly to GDAL-based R packages:

```{r, eval = FALSE}
## With vapour package
library(vapour)
data <- vapour_read_raster(
  "my_raster.tif",
  window = rasterio0(dimension, extent, subextent)
)

## With stars package
library(stars)
rio <- rasterio0(dimension, extent, subextent)
data <- stars::read_stars(
  "my_raster.tif",
  RasterIO = list(
    nXOff = rio[1] + 1,  # stars uses 1-based
    nYOff = rio[2] + 1,
    nXSize = rio[3],
    nYSize = rio[4]
  )
)
```

### Converting to sf-style RasterIO

The `rasterio_to_sfio()` function converts to the 1-based format used by
sf/stars:

```{r}
rio <- rasterio0(dimension, extent, subextent)
sfio <- rasterio_to_sfio(rio)
sfio
```

## VRT Files

GDAL Virtual Format (VRT) files are XML descriptions of raster datasets.
They're useful for creating virtual mosaics and defining subsets without
copying data.

### Extracting Extent from VRT

The `extent_vrt()` function parses a VRT file to extract its extent:

```{r, eval = FALSE}
## Parse VRT and get extent
vrt_file <- "path/to/mosaic.vrt"
ext <- extent_vrt(vrt_file)
ext
```

This is useful for understanding the coverage of a virtual mosaic before
reading data.

## Practical Example: Pre-computing Read Parameters

A common workflow is to pre-compute read parameters for efficient data access:

```{r}
## Known raster properties (e.g., from gdalinfo)
full_dimension <- c(43200, 21600)  # global 30 arc-second grid
full_extent <- c(-180, 180, -90, 90)

## Regions of interest
regions <- list(
  europe = c(-10, 40, 35, 70),
  australia = c(110, 155, -45, -10),
  amazon = c(-80, -45, -20, 5)
)

## Pre-compute RasterIO parameters for each region
read_params <- lapply(regions, function(roi) {
  aligned <- align_extent(roi, full_dimension, full_extent)
  crop_dim <- extent_dimension(aligned, full_dimension, full_extent)
  list(
    rasterio = rasterio0(full_dimension, full_extent, aligned),
    dimension = crop_dim,
    extent = aligned
  )
})

## Check Amazon parameters
read_params$amazon
```

This lets you plan I/O operations before touching the data.

## Practical Example: Creating Aligned Grids

When creating output grids compatible with GDAL tools:

```{r}
## Design a 1-degree global grid
target_res <- 1  # 1 degree cells
global_extent <- c(-180, 180, -90, 90)

## Compute dimensions from extent and resolution
dimension <- c(
  diff(global_extent[1:2]) / target_res,
  diff(global_extent[3:4]) / target_res
)
dimension

## Generate GeoTransform for GDAL
gt <- geo_transform0(dimension, global_extent)
gt

## Generate world file content
wf <- geo_world0(dimension, global_extent)
cat(wf, sep = "\n")
```

## Summary

| Function | Purpose |
|----------|---------|
| `geo_transform0()` | Convert to GDAL GeoTransform (6-element) |
| `geo_world0()` | Convert to world file format (6-element) |
| `gt_dim_to_extent()` | GeoTransform + dimension → extent |
| `extent_dim_to_gt()` | Extent + dimension → GeoTransform |
| `rasterio0()` | Compute GDAL RasterIO window parameters |
| `rasterio_to_sfio()` | Convert to sf/stars 1-based RasterIO |
| `extent_vrt()` | Extract extent from VRT file |

These functions bridge vaster's abstract grid logic with GDAL's concrete
representations, enabling you to plan operations and generate parameters
for GDAL-based tools.