comparisons with GeoBox

[mdsumner]

---

title: "Untitled"
format:
html:
code-fold: true
jupyter: python3

import odc.geo
from odc.geo.geobox import GeoBox
g = GeoBox.from_bbox(
   (-2_000_000, -5_000_000,
   2_250_000, -1_000_000),
   "epsg:3577", resolution=1000)
g.coordinates['y'].values

ex <- c(-2000000.0,2250000.0, -5000000.0, -1000000.0)
dm <- diff(ex)[c(1, 3)] %/% 1000
idx <- matrix(c(1, 3, 2, 3, 2, 4, 1, 4, 1, 3), ncol = 2L, byrow = TRUE)
xxyy <- matrix(ex[t(cbind(idx[-nrow(idx), ], idx[-1, ]))], ncol = 4, byrow = TRUE)
xxyy[,1:2] <- reproj::reproj_xy(xxyy[,1:2, drop = FALSE], "OGC:CRS84", source = "EPSG:3577")
xxyy[,3:4] <- reproj::reproj_xy(xxyy[,3:4, drop = FALSE], "OGC:CRS84", source = "EPSG:3577")

plot(rbind(xxyy[,1:2], xxyy[,3:4]), ylim = c(-50, -5), asp = 1/cos(mean(ex[3:4] * pi/180)))

for (i in 1:nrow(xxyy)) {
  ll <- bigcurve:::bisect(matrix(xxyy[i, , drop = FALSE], nrow = 2, byrow = T), "EPSG:3577", dist = min(diff(ex)[c(1, 3)]/50000))
  lines(ll)
  points(ll)
}

maps::map(add = TRUE)

library(reticulate)
odc.geo <- import("odc.geo")
g <- odc.geo$geobox$GeoBox$from_bbox(ex[c(1, 3, 2, 4)], crs = "EPSG:3577", resolution = 1000)
boundary <- g$boundary(100L)
## convert back through the affine
x <- boundary[,1] * 1000 + ex[1]
y <- boundary[,2] * -1000 + ex[4]
points(reproj::reproj_xy(cbind(x, y), "OGC:CRS84", source = "EPSG:3577"), pch = 19, col = "red")


lines(reproj::reproj_xy(vaster::vaster_boundary(dm %/% 20, ex), "OGC:CRS84", source = "EPSG:3577"), lwd = 2)

[mdsumner]

this is better, but still my crappy R needs a cleanup - I understand what GeoBox is doing now, and it's match by new vaster_boundary()

ex <- c(-2000000.0,2250000.0, -5000000.0, -1000000.0)
dm <- diff(ex)[c(1, 3)] %/% 1000
idx <- matrix(c(1, 3, 2, 3, 2, 4, 1, 4, 1, 3), ncol = 2L, byrow = TRUE)
xxyy <- matrix(ex[t(cbind(idx[-nrow(idx), ], idx[-1, ]))], ncol = 4, byrow = TRUE)
xxyy[,1:2] <- reproj::reproj_xy(xxyy[,1:2, drop = FALSE], "OGC:CRS84", source = "EPSG:3577")
xxyy[,3:4] <- reproj::reproj_xy(xxyy[,3:4, drop = FALSE], "OGC:CRS84", source = "EPSG:3577")

plot(rbind(xxyy[,1:2], xxyy[,3:4]), ylim = c(-50, -5), asp = 1/cos(mean(ex[3:4] * pi/180)))

## controlledburn is the best great circle by fewest points _in_that_crs_
for (i in 1:nrow(xxyy)) {
  ll <- bigcurve:::bisect(matrix(xxyy[i, , drop = FALSE], nrow = 2, byrow = T), "EPSG:3577", dist = min(diff(ex)[c(1, 3)]/5000))
  lines(ll)
  points(ll)
}

maps::map(add = TRUE)

library(reticulate)
odc.geo <- import("odc.geo")
g <- odc.geo$geobox$GeoBox$from_bbox(ex[c(1, 3, 2, 4)], crs = "EPSG:3577", resolution = 1000)
boundary <- g$boundary(100L)
## convert back through the affine
x <- boundary[,1] * 1000 + ex[1]
y <- boundary[,2] * -1000 + ex[4]
points(reproj::reproj_xy(cbind(x, y), "OGC:CRS84", source = "EPSG:3577"), pch = 19, col = "red")


lines(reproj::reproj_xy(vaster::vaster_boundary(dm %/% 20, ex), "OGC:CRS84", source = "EPSG:3577"), lwd = 2)

aa <- reproj::reproj_xy(cbind(ex[c(1, 2)], ex[3]), "OGC:CRS84", source = "EPSG:3577")
lines(geosphere::gcIntermediate(aa[1, ], aa[2, ]), col = "green")
bb <- reproj::reproj_xy(cbind(ex[c(1, 2)], ex[4]), "OGC:CRS84", source = "EPSG:3577")
lines(geosphere::gcIntermediate(bb[1, ], bb[2, ]), col = "green")