2m_pixels_at_2m_tall.md

library(arrow)
#> 
#> Attaching package: 'arrow'
#> The following object is masked from 'package:utils':
#> 
#>     timestamp
library(terra)
#> terra 1.7.71
#> 
#> Attaching package: 'terra'
#> The following object is masked from 'package:arrow':
#> 
#>     buffer
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:terra':
#> 
#>     intersect, union
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
## enter an address
addressname <- c("2862",  "Lyell Highway", "HAYES")
#addressname <- c("381",  "Elizabeth Street", "North Hobart")


addressname <- toupper(addressname)


srcs <- read_parquet("https://github.com/mdsumner/listmap.sources/raw/main/data-raw/listmap-sources-vsi.parquet") |> 
  distinct(dsn, .keep_all = TRUE)

## 2m dem
dsn <- "/vsicurl/https://s3.us-west-2.amazonaws.com/us-west-2.opendata.source.coop/alexgleith/tasmania-dem-2m/Tasmania_Statewide_2m_DEM_14-08-2021.tif"

## index for the address source
id1<- grep("ADDRESS", srcs$dsn)
base1 <- srcs$dsn[id1]
file1 <- srcs$vsilist[id1]
dsnaddress <- sprintf("%s/%s", base1, file1)
layeraddress <- vector_layers(dsnaddress)[1]

ST_NO <- as.integer(addressname[1])
STREET_TYPE <- tail(strsplit(addressname[2], "\\s+")[[1]], 1)
STREET <- trimws(gsub(sprintf("%s$", STREET_TYPE), "", addressname[2]))

LOCALITY <- addressname[3]

## get address point
address  <- terra::vect(dsnaddress, query = sprintf('SELECT * FROM "%s" WHERE STREET = \'%s\' AND ST_NO_FROM = %i AND LOCALITY = \'%s\'" ', 
                                                    layeraddress,  STREET, ST_NO , LOCALITY))
#> Warning: Did not find end-of-string character (GDAL error 1)

#> Warning: Did not find end-of-string character (GDAL error 1)


## read all the municipalities
munic <- grep("LIST_LOCALITY_POSTCODE", srcs$dsn, value = TRUE)
## municipality is always the second layer, we only need its NAME
municipality <- do.call(rbind, lapply(munic, \(.x) vect(.x, vector_layers(.x)[2])))
municipality <- municipality[which(relate(municipality, address, "intersects")), ]

## now get the parcel source from the address LOCALITY
parcelsource <- grep(sprintf("PARCELS_%s", toupper(gsub("\\s+", "_", municipality$NAME))), srcs$dsn)

## index for the parcel source

base <- srcs$dsn[parcelsource]
file <- srcs$vsilist[parcelsource]

dsnparcel <- sprintf("%s/%s", base, file)
layerparcel <- vector_layers(dsnparcel)


parcel <- terra::vect(dsnparcel, query = sprintf("SELECT * FROM %s WHERE PID = %s", layerparcel, address$PID[1]))

dem <- project(rast(dsn), rast(align(ext(as.vector(ext(parcel))), 2), res = 2, crs = crs(parcel)), by_util = TRUE)
plot(dem)
plot(parcel, add = TRUE)
plot(address, add = TRUE)
text(geom(address)[, c("x", "y"), drop = FALSE], pos = 1, lab = paste(address$ST_NO_FROM, address$STREET, address$LOCALITY))

^{Created on 2024-04-07 with reprex v2.0.2}