Как проецировать пространственные точки на Аляске на карту США, которая масштабировала и преобразовала Аляску в R - PullRequest
Новая
       69

Как проецировать пространственные точки на Аляске на карту США, которая масштабировала и преобразовала Аляску в R

0 голосов
/ 21 мая 2018

Я пытаюсь построить карту США с помощью ALASKA, а затем спроецировать на нее набор точек.Я использовал опубликованные функции fixup и fix1 в следующей ссылке, чтобы переместить Аляску и Гавайи.

Перемещение Аляски и Гавайев на тематической карте США с помощью ggplot2

Затем спроецируйте мои точки со значениями широты и долготы как dput ниже 

my_points= structure(c(44.334567, 44.209571, 44.049845, 44.752622, 44.791511, 
44.391792, 44.540124, 46.075669, 46.007611, 45.836781, 46.595384, 
45.703999, 45.47594, 44.715117, 44.385954, 42.804004, 43.349842, 
43.252618, 42.891499, 56.212978, 55.391604, 55.395194, 60.476929, 
61.709743, 62.76729, 62.346439, 59.93483, 64.472359, 64.88569, 
-122.047007, -122.256733, -123.42621, -122.128965, -122.578974, 
-122.497582, -122.435915, -121.998698, -122.347319, -122.466208, 
-122.459556, -123.755405, -123.725121, -123.887335, -123.831778, 
-123.610907, -122.728942, -123.026172, -124.070652, -131.638358, 
-131.195583, -132.408627, -151.081668, -149.231938, -149.693379, 
-150.019201, -158.190006, -146.926265, -147.249648), .Dim = c(29L, 
2L))

Мой код: 

require(maptools)
require(rgdal)

fixup <- function(usa,alaskaFix,hawaiiFix){
  alaska=usa[usa$STATE_NAME=="Alaska",]
  alaska = fix1(alaska,alaskaFix)
  proj4string(alaska) <- proj4string(usa)
  hawaii = usa[usa$STATE_NAME=="Hawaii",]
  hawaii = fix1(hawaii,hawaiiFix)
  proj4string(hawaii) <- proj4string(usa)
  usa = usa[! usa$STATE_NAME %in% c("Alaska","Hawaii"),]
  usa = rbind(usa,alaska,hawaii)
  return(usa)
}

fix1 <- function(object,params){
  r=params[1];scale=params[2];shift=params[3:4]
  object = elide(object,rotate=r)
  size = max(apply(bbox(object),1,diff))/scale
  object = elide(object,scale=size)
  object = elide(object,shift=shift)
  object
}

## download the USA shapefile to the computer
    setwd(tempdir())
download.file("https://dl.dropbox.com/s/wl0z5rpygtowqbf/states_21basic.zip?dl=1", 
              "usmapdata.zip", 
              method = "curl")

 #This is a mirror of http://www.arcgis.com/home/item.html?
 #id=f7f805eb65eb4ab787a0a3e1116ca7e5
unzip("usmapdata.zip")

Тогда читайте в моем shapefile.Используйте rgdal: 

us = readOGR(dsn = "states_21basic",layer="states")

Теперь преобразуйте в равную область и запустите функцию исправления: 

usAEA = spTransform(us,CRS("+init=epsg:2163"))
usfix = fixup(usAEA,c(-35,1.5,-2800000,-2600000),c(-35,1,6800000,-1600000))
plot(usfix)

Параметры: вращения, масштабирование, сдвиг x и y для Аляски иГавайи соответственно и были получены методом проб и ошибок.Тщательно настройте их.Даже изменение параметра масштаба Гавайев на 0,99999 отослало его с планеты из-за больших чисел.

Поверните это обратно к широте: 

usfixLL = spTransform(usfix,CRS("+init=epsg:4326"))
plot(usfixLL)

Теперь, работая над моими точками, чтобы добавитьих на производимую карту 

require(sp)
my_NEW_S<-as.data.frame(my_points)
colnames(my_NEW_S)<-c("y","x")
coordinates(my_NEW_S)=~x+y #column names of the lat long cols
proj4string(my_NEW_S)<- CRS("+init=epsg:2163")

# add my points to the produced map
points(my_NEW_S,col="red",pch=16,cex=0.5)

Все точки совпадающих США проецируются на карту, КРОМЕ ТОГО связаны с Аляской и HI.Было бы очень признательно, если бы кто-нибудь помог мне с тем, как сделать мою проекцию очков Аляски такой же, как и перемещенную Аляску, или любое другое решение.Спасибо.

...