CODEWITHSUNDEEP
rrasterr-rasterr-sprgdal
Diego Moya
Diego Moya

Reputation: 143

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.

Upvotes: 1

Views: 764

Answers (2)

lbusett
lbusett

Reputation: 5932

As an alternative, consider that package spData provides polygon boundaries for the US already splitted for contigouous US, Alaska, Hawaii and other islands that you can use directly. For example, using an elevation raster downloaded with raster::getData:

library(spData)
library(sf)
library(raster)

elev   <- raster::getData("alt", country="USA", level=1) 
#> returning a list of RasterLayer objects
usa    <- spData::us_states %>% sf::st_transform(4326)
hawaii <- spData::hawaii %>% sf::st_transform(4326)
alaska <- spData::alaska %>% sf::st_transform(4326)

# crop raster to area of interest and plot (Note that the data downloaded with `raster::getData`
# is split in four subdatasets, so in this case you need to select the correct one.)
usa_elev <- crop(elev[[1]], usa)
plot(usa_elev)
plot(sf::st_geometry(usa), add = TRUE)


alaska_elev <- crop(elev[[2]], alaska)
plot(alaska_elev)
plot(st_geometry(alaska), add = TRUE)


hawaii_elev <- crop(elev[[4]], hawaii)
plot(hawaii_elev)
plot(st_geometry(hawaii), add = TRUE)

Created on 2019-01-04 by the reprex package (v0.2.1)

Upvotes: 3

Diego Moya
Diego Moya

Reputation: 143

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

Upvotes: 0

Related Questions