generating aerial images with gdalio

test_gdalio.R

library(sf)
library(terra)
library(gdalio)
library(geobr)

### Choose either a small or large area

i = 49 # small area

i = 1066 # large area



###### 1. Get census tract data --------------------------------

# read local table data with all census tracts
all_cts <- geobr::read_census_tract(code_tract = "SE")

all_cts$area <- st_area(all_cts) |> as.numeric()

# subset census tract
temp_ct <- all_cts[i,]

# fix eventual topoly errors
temp_ct <- sf::st_make_valid(temp_ct)


###### 2. gdalio --------------------------------

## {terra}
gdalio_terra <- function(dsn, ..., band_output_type = "numeric") {
  v <- gdalio_data(dsn, ..., band_output_type  = band_output_type)
  g <- gdalio_get_default_grid()
  r <- terra::rast(terra::ext(g$extent), nrows = g$dimension[2], ncols = g$dimension[1], crs = g$projection)
  if (length(v) > 1) terra::nlyr(r) <- length(v)
  terra::setValues(r, do.call(cbind, v))
}

virtualearth_imagery <- tempfile(fileext = ".xml")

writeLines('<GDAL_WMS>
  <Service name="VirtualEarth">
    <ServerUrl>http://a${server_num}.ortho.tiles.virtualearth.net/tiles/a${quadkey}.jpeg?g=90</ServerUrl>
  </Service>
  <MaxConnections>4</MaxConnections>
  <Cache/>
</GDAL_WMS>', virtualearth_imagery)


### calculate area extension

# bounding box
bb <- sf::st_bbox(temp_ct)

# width length
point12 <- st_point(c(bb[1], bb[2])) |> st_sfc(crs = st_crs(temp_ct)) 
point32 <- st_point(c(bb[3], bb[2])) |> st_sfc(crs = st_crs(temp_ct)) 
x_dist <- st_distance(point12, point32) |> as.numeric()

# height length
point12 <- st_point(c(bb[1], bb[2])) |> st_sfc(crs = st_crs(temp_ct)) 
point14 <- st_point(c(bb[1], bb[4])) |> st_sfc(crs = st_crs(temp_ct)) 
y_dist <- st_distance(point12, point14) |> as.numeric()

# extension
ext <- (max(x_dist, y_dist) * 1.05) |> round()
ext <- ifelse(ext < 2000, 2000, ext)

### pick dimensions (resolution)
ct_area <- sf::st_area(temp_ct) |> as.numeric()
dim <- ifelse(ct_area > 5000000, 5000, 9000)

# centroid coordinates
coords <- st_centroid(temp_ct) |> st_coordinates()

my_zoom <- paste0("+proj=laea +lon_0=", coords[1], " +lat_0=", coords[2])


### prepare grid query
grid1 <- list(extent = c(-1, 1, -1, 1) * ext,
              dimension = c(dim, dim), 
              projection = my_zoom)

gdalio_set_default_grid(grid1)

img <- gdalio_terra(virtualearth_imagery, bands = 1:3)
# terra::plotRGB(img)


###### 3. mask, crop and save image --------------------------------

# reproject census tract
temp_ct2 <- st_transform(temp_ct, crs= st_crs(img))

# mask raster with census tract
temp_plot <- terra::mask(img, temp_ct2)

# crop with buffer of 80 meters
buff <- sf::st_buffer(temp_ct2, dist = 50)
temp_plot <- terra::crop(temp_plot, buff)
# terra::plotRGB(temp_plot)

# image proportions
r <- nrow(temp_plot) /ncol(temp_plot)

# save image to tempfile
tempd <- tempdir()
image_file <- paste0(tempd, '/image_',i,'.png')

png(image_file, res = 500, width = 15, height = 15*r, units = 'cm') 
raster::plotRGB(temp_plot)
dev.off()

# open image locally
utils::browseURL(image_file)

Here's a little bit simpler, I use the lonlat centroid to generate the projected version, and just use the extent of that no need for sf/distance etc. There's also no need to write the xml to file (I learnt relatively recently) so I've also changed that.

I've set it to small dimensions so it runs fast, and turned off the write to file - just so you can see immediately what it'sdoing and then you can choose the resolution you want. It could be set to give a specific resolution (i.e. pixel size) by simply dividiing the as.integer(diff(ext)[c(1, 3)]) by the width/height you want in metres of each pix.

So, you could even use a constant projection for all of Brazil rather than by the local region, but this will work globally.

library(sf)
library(terra)
library(gdalio)
library(geobr)

### Choose either a small or large area

#i = 49 # small area

i = 1066 # large area

###### 1. Get census tract data --------------------------------

# read local table data with all census tracts
all_cts <- geobr::read_census_tract(code_tract = "SE")

all_cts$area <- st_area(all_cts) |> as.numeric()

# subset census tract
temp_ct <- all_cts[i,]

# fix eventual topoly errors
temp_ct <- sf::st_make_valid(temp_ct)


###### 2. gdalio --------------------------------

## {terra}
gdalio_terra <- function(dsn, ..., band_output_type = "numeric") {
  v <- gdalio_data(dsn, ..., band_output_type  = band_output_type)
  g <- gdalio_get_default_grid()
  r <- terra::rast(terra::ext(g$extent), nrows = g$dimension[2], ncols = g$dimension[1], crs = g$projection)
  if (length(v) > 1) terra::nlyr(r) <- length(v)
  terra::setValues(r, do.call(cbind, v))
}

virtualearth_imagery <- '<GDAL_WMS>
  <Service name="VirtualEarth">
    <ServerUrl>http://a${server_num}.ortho.tiles.virtualearth.net/tiles/a${quadkey}.jpeg?g=90</ServerUrl>
  </Service>
  <MaxConnections>4</MaxConnections>
  <Cache/>
</GDAL_WMS>'


### calculate area extension
# centroid coordinates
coords <- st_centroid(temp_ct) |> st_coordinates()

prj <- paste0("+proj=laea +lon_0=", coords[1], " +lat_0=", coords[2])
temp_ct_prj <- sf::st_transform(temp_ct, prj)
# bounding box
ext <- sf::st_bbox(temp_ct_prj)[c(1, 3, 2, 4)]

### pick dimensions (resolution)
ct_area <- sf::st_area(temp_ct_prj) |> as.numeric()
dim <- ifelse(ct_area > 5000000, 5000, 9000)/8


### prepare grid query
grid1 <- list(extent = ext,
              dimension = c(dim, dim), 
              projection = prj)

gdalio_set_default_grid(grid1)

img <- gdalio_terra(virtualearth_imagery, bands = 1:3)
# terra::plotRGB(img)


###### 3. mask, crop and save image --------------------------------

# mask raster with census tract
temp_plot <- terra::mask(img, temp_ct_prj)

# crop with buffer of 80 meters
buff <- sf::st_buffer(temp_ct_prj, dist = 50)
temp_plot <- terra::crop(temp_plot, buff)
terra::plotRGB(temp_plot)

# # image proportions
# r <- nrow(temp_plot) /ncol(temp_plot)
# 
# # save image to tempfile
# tempd <- tempdir()
# image_file <- paste0(tempd, '/image_',i,'.png')
# 
# png(image_file, res = 500, width = 15, height = 15*r, units = 'cm') 
# raster::plotRGB(temp_plot)
# dev.off()
# 
# # open image locally
# utils::browseURL(image_file)