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)
title: "Untitled"
format:
html:
code-fold: true
jupyter: python3