Нужно: обрезать / выбрать / вырезать / установить (или просто показать) значения растра США, включая смежные США, Аляску и Гавайи

Question

ГИС, мне нужно обрезать / выбрать / вырезать / подмножество (или просто показать) значения растра США, включая материк, Аляску и Гавайи.Это сбивает с толку видеть большую карту, включая некоторые острова или территории очень далеко.Итак, я пытался выбрать / вырезать растр, чтобы включить только материк США, Аляску и Гавайи, а затем сделать визуализацию.Код, который я разработал, выглядит следующим образом:

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)

Выходные данные:

Воспроизводимый Пример :

Нужно что-то вроде этого:

Я уже пробовал это , это и это Те.

Любая помощь очень ценится.

Answer 1

В качестве альтернативы рассмотрим, что пакет spData предоставляет границы полигонов для США, уже разделенных на смежные США, Аляску, Гавайи и другие острова, которые вы можете использовать напрямую.Например, с использованием растра высот, загруженного с помощью 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)

^{Создано в 2019-01-04 пакетом Представить (v0.2.1)}

Answer 2

Я решил это следующим образом:

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))

Я получил это:

Надеюсь, что это может помочь другим