R: построение шейп-файлов

Question

Я очень плохо знаком с R и сейчас пытаюсь создать карту из шейп-файла. Однако я застрял на том, как заполнить карту цветом на основе заданного параметра в моем наборе данных ... Я был бы признателен за любую помощь! Спасибо!

Мой код выглядит следующим образом:

#Read in data and set variables
parameter_results<- readRDS("param_results_2014.RDS")
NJ<- readOGR(dsn="V:/lum/WM&S/BEAR (Bureau of Environmental Analysis and Restoration)/Envpln/Hourly Employees/JohnDoe/Rwork/2014IR/Maps/shapefiles",layer="2014_NJ_Integrated_Report_AU")
#join common field
NJ_merge<- merge(NJ,parameter_results,by.x="AU_NUM",by.y="Waterbody")
colors <- c("#91D79E","#FFFF73","#FF7F7F","#FF7F7F")


#Following lines create plot with scale and arrow
plotexpression<-plot(NJ_merge, xlim=c(200000.732,905000.646), ylim=c(-5812.321,900000.543),main = "Figure 2.4A: Assessment Results for Trout Aquatic\nLife Use, Spatial Extent",col=NJ_merge$watsup)

prettymap(plotexpression,oma=c(3,3,4,3),drawbox = TRUE,scale.plotunit="mi",drawscale = TRUE,scale.pos = "bottomright",drawarrow = TRUE,arrow.scale = .5,box.lwd = 1, arrow.cols = c("black","black"),arrow.text.col = "black")


legend(640000, 400000, legend=c("Fully Supported", "Insufficient Data", "TMDL Waters (Not Supported)", "Not Supporting","Not Applicable"), 
       fill=c(colors), bty="n", title = "Aquatic Life - Trout\nDesignated Use 2014\nAssessment",text.font = 2, cex=0.6,title.adj=0.2,title.col=1)

Результаты параметра набора данных:

# A tibble: 958 x 89
   WMA   Waterbody  Name    `Biological (Cau~ `Biological Trou~ DO    `DO Trout` Temperature
   <chr> <chr>      <chr>   <chr>             <chr>             <chr> <chr>      <chr>      
 1 15    020403020~ Abseco~ Attaining         Not Applicable    Atta~ Not Appli~ Attaining  
 2 15    020403020~ Abseco~ Insufficient Inf~ Not Applicable    Non ~ Not Appli~ Attaining  
 3 15    020403020~ Abseco~ Attaining         Not Applicable    Insu~ Not Appli~ Insufficie~
 4 15    020403020~ Abseco~ Attaining         Not Applicable    Atta~ Not Appli~ Attaining  
 5 14    020403011~ Albert~ Non Attaining     Not Applicable    Atta~ Not Appli~ Attaining  
 6 11    020401052~ Alexau~ Attaining         Attaining         Insu~ Attaining  Insufficie~
 7 11    020401052~ Alexau~ Attaining         Attaining         Insu~ Attaining  Insufficie~
 8 17    020402060~ Allowa~ Non Attaining     Not Applicable    Atta~ Not Appli~ Attaining  
 9 17    020402060~ Allowa~ Insufficient Inf~ Not Applicable    Atta~ Not Appli~ Attaining  
10 17    020402060~ Allowa~ Insufficient Inf~ Not Applicable    Insu~ Not Appli~ Insufficie~

NJ_merge выглядит следующим образом:

str(NJ_merge,2)
Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
  ..@ data       :'data.frame': 958 obs. of  101 variables:
  ..@ polygons   :List of 958
  ..@ plotOrder  : int [1:958] 950 844 853 421 687 329 334 721 251 321 ...
  ..@ bbox       : num [1:2, 1:2] 190378 10574 659480 919549
  .. ..- attr(*, "dimnames")=List of 2
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot

Сюжет, который я сейчас получаю:

Что мне не хватает, чтобы заполнить области на карте цветами, которые есть в легенде ???

Answer 1

Это очень просто сделать с пакетом sf. Используя включенные данные шейп-файла для Северной Каролины, я создаю групповую переменную и отображаю ее с помощью либо включенного метода plot.sf, либо geom_sf, включенного в версию для разработки ggplot2.

library(tidyverse)
library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.2.3, proj.4 4.9.3
set.seed(100)
nc <- system.file("shape/nc.shp", package="sf") %>%
  st_read() %>%
  mutate(
    group = sample.int(5, 100, replace = TRUE),
    group = parse_factor(group, levels = 1:5)
    )
#> Reading layer `nc' from data source `C:\Users\Calum You\Documents\R\win-library\3.4\sf\shape\nc.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 100 features and 14 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> epsg (SRID):    4267
#> proj4string:    +proj=longlat +datum=NAD27 +no_defs

plot(nc[, "group"])

ggplot(nc, aes(fill = group)) +
  theme_minimal() +
  geom_sf() +
  scale_fill_brewer(type = "qual")

Создано в 2018-04-24 пакетом представ. (v0.2.0).