Skip to content

start in using vaster to work with the index #4

@mdsumner

Description

@mdsumner
shp <- sf::read_sf("/vsicurl/https://github.com/wmgeolab/geoBoundaries/raw/main/releaseData/CGAZ/geoBoundariesCGAZ_ADM0.gpkg")
r <- controlledburn:::burn_polygon(shp, ext <- c(-180, 180, -90, 90), dm <- c(360, 180) * 4000)
## our index is triplets of start,end,line where the polygon edge was detected - 
## this essentially an rle by scanline of start,end polygon coverage
index <- matrix(unlist(r, use.names = F), ncol = 4L, byrow = TRUE) + 1 ## plus one because 0-index internally
## plot just the start and ends of each scanline detected
xy0 <- vaster::xy_from_cell(dm, ext, vaster::cell_from_row_col(dm, t(index[,c(3, 3)]), t(index[,1:2])))
plot(xy0, pch = ".")


xy1 <- xy0[seq(1, nrow(xy0), by = 2), ]
xy2 <- xy0[seq(2, nrow(xy0), by = 2), ]
segments(xy1[,1], xy1[,2], xy2[,1], xy2[,2], col = palr::d_pal(rep(index[,4], each = 1)))


plot(xy0, pch = ".", xlim = c(140, 155), ylim = c(-20, 0))
segments(xy1[,1], xy1[,2], xy2[,1], xy2[,2], col = palr::d_pal(rep(index[,4], each = 1)))

plot(xy0, pch = ".", xlim = c(151, 151.2), ylim = c(-8.8, -8.3))
segments(xy1[,1], xy1[,2], xy2[,1], xy2[,2], col = palr::d_pal(rep(index[,4], each = 1)))
plot(shp, add = T)

cr <- do.call(rbind, apply(index, 1,\(.x) cbind(seq(.x[1], .x[2]), .x[3], .x[4])))
xy <- vaster::xy_from_cell(dm, ext, vaster::cell_from_row_col(dm, cr[,2], cr[,1]))
plot(xy, pch = 19, cex = .3, col = palr::d_pal(cr[,3]))


## how would we use cells?

## 1) area_from_row determine area (easy for longlat)
## 2) group up values (like area weighting)
## 3) values from other datasets (the one we rasterized from its definition)
## 4) efficient plotting of the scan lines
## 5) efficient conversion to polygon regions from raster
## 6) ability to crop
## 7) ability to farm out the job (but that's easy, just loop over chunks)
area_from_row <- function(dm, ext, row) {
  SCALE <- 1852 * 60
  lat <- vaster::y_from_row(dm, ext, row)
  size <- vaster::x_res(dm, ext) * cos(lat * pi/180) * vaster::x_res(dm, ext) * SCALE
  size
}


a <- area_from_row(dm, ext, cr[,2])
plot(xy, pch = 19, cex = .3, col = palr::d_pal(a))


Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions