Search code examples
rrasterr-rasterr-sprgdal

Need: clip/select/cut/subset (or just show) the values of a USA raster including the Contiguous US, Alaska and Hawaii

GIS people, I need to clip/select/cut/subset (or just show) the values of a USA raster including the mainland, Alaska and Hawaii. It is confusing to see the big map including some Islands or territories very far away. So, I have been trying to select/cut the raster to only include USA mainland, Alaska and Hawaii and then do the visualization. The code I have developed is as follow:

library(rgdal)
library(raster)

state <- getData("GADM", country="USA", level=1)
projection(state) <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

nestates <- c("Alabama","Arizona", "Arkansas", "California", "Colorado",  # Contiguous/Continental United States
              "Connecticut", "Delaware", "Florida", "Georgia", "Idaho",
              "Illinois", "Indiana", "Iowa", "Kansas", "Kentucky", 
              "Louisiana", "Maine", "Maryland", "Massachusetts", 
              "Michigan", "Minnesota", "Mississippi", "Missouri", 
              "Montana", "Nebraska", "Nevada", "New Hampshire", "New 
               Jersey", "New Mexico", "New York", "North Carolina",
               "North Dakota", "Ohio", "Oklahoma", "Oregon", 
               "Pennsylvania", "Rhode Island","South Carolina", 
               "South Dakota", "Tennessee", "Texas", "Utah", "Vermont",
               "Virginia", "Washington", "West Virginia", "Wisconsin", 
               "Wyoming",
               "Alaska", "Hawaii") # I tried excluding Hawaii too

# I believe the issue is with insular territories 

state.sub <- state[as.character(state@data$STATE_NAME) %in% nestates, ]

elevation <- raster("USA_1.tif")

elevation.sub <- crop(elevation, extent(state.sub))

elevation.sub <- mask(elevation.sub, state.sub) # Error in x@polygons[[i]] : subscript out of bounds

plot(elevation.sub)
plot(state.sub, add = TRUE)

The output, so far: enter image description here

The reproducible example:

Need something like this: enter image description here

I already tried this, this and this ones.

Any help is very much appreciated.

Solution

  • I have solved it in this way:

    us<-getData('GADM', country='USA', level=1)  #Get the County Shapefile for the US

nestates <-         c("Alabama","Arizona","Arkansas","California","Colorado","Connecticut","Delaware","District of Columbia",
          "Florida","Georgia","Idaho","Illinois","Indiana","Iowa","Kansas","Kentucky","Louisiana","Maine","Maryland",          
          "Massachusetts","Michigan","Minnesota","Mississippi","Missouri","Montana","Nebraska","Nevada","New Hampshire",       
          "New Jersey","New Mexico","New York","North Carolina","North Dakota","Ohio","Oklahoma","Oregon","Pennsylvania",    
          "Rhode Island","South Carolina","South Dakota","Tennessee","Texas","Utah","Vermont","Virginia","Washington",        
          "West Virginia","Wisconsin","Wyoming")
#"Alaska" polygons include the far away Islands   #"Hawaii"
# I followed these tutorials/Q&A
#http://data-analytics.net/wp-content/uploads/2014/09/geo2.html
#https://gis.stackexchange.com/questions/61243/clipping-a-raster-in-r/61278

ne = us[match(toupper(nestates),toupper(us$NAME_1)),]
raster_c <- crop(raster_1, extent(ne))

    I got this:

    enter image description here

    Hope it might help others