I am trying to plot the USA map with ALASKA and then project a set of points on it. I used the posted functions fixup
and fix1
in the following link to relocate Alaska & Hawaii.
Relocating Alaska and Hawaii on thematic map of the USA with ggplot2
Then project my points that have latitude and longitude values as dput
below
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))
My code is:
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")
Then read in my shapefile
. Use rgdal
:
us = readOGR(dsn = "states_21basic",layer="states")
Now transform to equal-area, and run the fixup function:
usAEA = spTransform(us,CRS("+init=epsg:2163"))
usfix = fixup(usAEA,c(-35,1.5,-2800000,-2600000),c(-35,1,6800000,-1600000))
plot(usfix)
The parameters are rotations, scaling, x and y shift for Alaska and Hawaii respectively, and were obtained by trial and error. Tweak them carefully. Even changing Hawaii's scale parameter to 0.99999 sent it off the planet because of the large numbers involved.
Turn this back to lat-long:
usfixLL = spTransform(usfix,CRS("+init=epsg:4326"))
plot(usfixLL)
Now, working on my points to add them to the produced map
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)
All the points of conterminous US are projected on the map EXCEPT the ones are related to Alaska and HI. It would be greatly appreciated if someone helps me with how to make my Alaska points projection as same as the relocated Alaska, or any other solution. Thank you.