2

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.

J.Con
  • 4,101
  • 4
  • 36
  • 64
Ali
  • 21
  • 1
  • use the github-only `albersusa` package and `points_elided()` https://github.com/hrbrmstr/albersusa/blob/master/R/points_elided.R – hrbrmstr May 20 '18 at 22:57
  • Thank you for your help. I have used the function points_elided it self because I was unable to install the package albersusa in my R version 3.4.2. But the function gave an error of "nrow(coords) > 0 is not TRUE". The error occurred when I started using maptools::elide function. And I noticed in the previous step that the sp::over wasn't able to subset the Alaska points. Please any help with that? – Ali May 21 '18 at 15:47

0 Answers0