#!mathjax#
library('mapmisc', quietly=TRUE)
library('terra')
if(!exists('fact'))
fact = 1
# uncomment for nicer images
# fact = 2
## [1] "/var/folders/1s/zkmc02qn4k18r6jdtbb459hc0000gn/T//RtmpmPTDOX"
## [1] FALSE
worldMap = unwrap(worldMap)
x = vect(cbind(130,15), crs=crsLL)
if(FALSE) {
okinawa = geocode("Okinawa, Japan")
hk = geocode("Hong Kong, China")
guam = geocode("Guam")
} else {
okinawa = unwrap(new("PackedSpatVector", type = "points", crs = crsLL,
coordinates = structure(c(127.911468664084, 26.4748948), dim = 1:2, dimnames = list(
NULL, c("x", "y"))), index = structure(c(1, 1, 0, 1), dim = c(1L,
4L), dimnames = list(NULL, c("geom", "part", "hole", "start"
))), attributes = structure(list(name = "Okinawa Island",
orig = "Okinawa, Japan", type = "island", category = "place",
importance = 0.544839401260831, place_id = 298674558L,
osm_type = "relation", osm_id = 4098955L, display_name = "沖縄本島, 沖縄県, 日本",
place_rank = 17L), class = "data.frame", row.names = c(NA, -1L))))
hk=unwrap(new("PackedSpatVector", type = "points", crs = crsLL,
coordinates = structure(c(114.1849161, 22.350627), dim = 1:2, dimnames = list(
NULL, c("x", "y"))), index = structure(c(1, 1, 0, 1), dim = c(1L,
4L), dimnames = list(NULL, c("geom", "part", "hole", "start"
))), attributes = structure(list(name = "Hong Kong", orig = "Hong Kong, China",
type = "administrative", category = "boundary", importance = 0.964301726010044,
place_id = 368735620L, osm_type = "relation", osm_id = 913110L,
display_name = "香港 Hong Kong, 中国", place_rank = 6L), class = "data.frame", row.names = c(NA, -1L))))
guam = unwrap(new("PackedSpatVector", type = "points", crs = crsLL,
coordinates = structure(c(144.757551029721, 13.4501257), dim = 1:2, dimnames = list(
NULL, c("x", "y"))), index = structure(c(1, 1, 0, 1), dim = c(1L,
4L), dimnames = list(NULL, c("geom", "part", "hole", "start"
))), attributes = structure(list(name = "Guam", orig = "Guam",
type = "island", category = "place", importance = 0.799784943797409,
place_id = 298225988L, osm_type = "relation", osm_id = 1152532L,
display_name = "Guam, Chalan Pago-Ordot Municipality, Guam, United States",
place_rank = 17L), class = "data.frame", row.names = c(NA, -1L))))
}
myProj = tpers(x=vect(c(x,hk)), hKm = 2*1000, tilt =-15, axis='seu')
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Reprojection failed, err = 2050, further errors will be suppressed on
## the transform object. (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Reprojection failed, err = 2050, further errors will be suppressed on
## the transform object. (GDAL error 1)
map.new(attributes(myProj)$regionLL, buffer=2)
plot(attributes(myProj)$crop, add=TRUE, col='red')
plot(project(worldMap, crsLL), add=TRUE)
text(guam, 'guam')
text(hk, 'hk')
text(okinawa, 'Okinawa')
xProj = project(guam, myProj)
yProj = project(hk, myProj)
okinawaProj = project(okinawa, myProj)
myMap = openmap(myProj, zoom=2,fact=2)
map.new(myProj, lwd=3)
plot(myMap, add=TRUE)
plot(wrapPoly(worldMap, myProj), add=TRUE)
myProjDateLine = tpers(x=c(175,75), azi=40, hKm = 2*1000, tilt =-15, axis='seu')
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Reprojection failed, err = 2050, further errors will be suppressed on
## the transform object. (GDAL error 1)
myMapDateLine = openmap(myProjDateLine, zoom=2, fact=1)
map.new(myProjDateLine)
plot(myMapDateLine, add=TRUE)
plot(wrapPoly(worldMap, myProjDateLine), add=TRUE)
myProjDateLine = tpers(x=c(-165,0), azi=-70, hKm = 30*1000, tilt =-15, axis='seu')
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Reprojection failed, err = 2050, further errors will be suppressed on
## the transform object. (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Reprojection failed, err = 2050, further errors will be suppressed on
## the transform object. (GDAL error 1)
myMapDateLine = openmap(myProjDateLine, zoom=2, fact=1, verbose=TRUE)
## http://tile.openstreetmap.org/
## http://tile.openstreetmap.org/
## tile /var/folders/1s/zkmc02qn4k18r6jdtbb459hc0000gn/T//RtmpmPTDOX/org/2/0/0.png cached
## tile /var/folders/1s/zkmc02qn4k18r6jdtbb459hc0000gn/T//RtmpmPTDOX/org/2/0/1.png cached
## downloading http://tile.openstreetmap.org/2/0/2.png
## downloading http://tile.openstreetmap.org/2/0/3.png
## tile /var/folders/1s/zkmc02qn4k18r6jdtbb459hc0000gn/T//RtmpmPTDOX/org/2/1/1.png cached
## downloading http://tile.openstreetmap.org/2/1/2.png
## tile /var/folders/1s/zkmc02qn4k18r6jdtbb459hc0000gn/T//RtmpmPTDOX/org/2/3/0.png cached
## tile /var/folders/1s/zkmc02qn4k18r6jdtbb459hc0000gn/T//RtmpmPTDOX/org/2/3/1.png cached
## tile /var/folders/1s/zkmc02qn4k18r6jdtbb459hc0000gn/T//RtmpmPTDOX/org/2/3/2.png cached
## downloading http://tile.openstreetmap.org/2/3/3.png
## reprojecting
## reprojecting: 3 cycles:1 2 2
map.new(myProjDateLine)
plot(myMapDateLine, add=TRUE)
plot(wrapPoly(worldMap, myProjDateLine), add=TRUE)
myProjDateLine = tpers(x=c(175,-85), azi=40, hKm = 20*1000, tilt =-15, axis='seu')
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Point outside of projection domain (GDAL error 1)
## Warning: Reprojection failed, err = 2050, further errors will be suppressed on
## the transform object. (GDAL error 1)
myMapDateLine = openmap(myProjDateLine, zoom=2)
map.new(myProjDateLine)
plot(myMapDateLine, add=TRUE)
plot(wrapPoly(worldMap, myProjDateLine), add=TRUE)
map.new(myProj, bg='black', buffer=2*100*1000)
plot(myMap, add=TRUE)
plot(wrapPoly(worldMap, myProj), add=TRUE)
gridlinesWrap(myProj, col='orange',
norths=seq(-40,40,by=10), easts=seq(-150,180,by=10),
lty=2)
points(xProj, col='red', pch=16, cex=2)
text(xProj,
as.character(xProj$name),
pos=2, col='red')
points(okinawaProj, col='red', pch=16, cex=2)
text(okinawaProj,
as.character(okinawaProj$name),
pos=2, col='red')
points(yProj, col='red', pch=16, cex=2)
text(yProj,
as.character(yProj$name),
pos=2, col='red')