Использование R: Работа с геопространственными данными (и ggplot2)

Загрузка библиотек и данных

# load libraries
library(ggplot2)
library(sp)
library(rgdal)
library(rgeos)

# create a local directory for the data
localDir <- "R_GIS_data"
if (!file.exists(localDir)) {
  dir.create(localDir)
}

# download and unzip the data
url <- "ftp://www.ecy.wa.gov/gis_a/inlandWaters/wria.zip"
file <- paste(localDir, basename(url), sep='/')
if (!file.exists(file)) {
  download.file(url, file)
  unzip(file,exdir=localDir)
}

# create a layer name for the shapefiles (text before file extension)
layerName <- "WRIA_poly"

# read data into a SpatialPolygonsDataFrame object
dataProjected <- readOGR(dsn=localDir, layer=layerName)

Преобразование данных

Пока что мы не сделали ничего радикально отличного от предыдущего, но для того, чтобы подготовить данные для построения графика в ggplot, нам нужно будет сделать пару манипуляций со структурой данных. ggplot2будет работать только с data.frameобъектом, поэтому наш объект класса SpatialPolygonsDataFrameне будет подходящим для построения. Давайте напишем некоторый код и обсудим, почему этот вид преобразования необходим.

# add to data a new column termed "id" composed of the rownames of data
dataProjected@data$id <- rownames(dataProjected@data)

# create a data.frame from our spatial object
watershedPoints <- fortify(dataProjected, region = "id")

# merge the "fortified" data with the data from our spatial object
watershedDF <- merge(watershedPoints, dataProjected@data, by = "id")

# NOTE : If we so choose, we could have loaded the plyr library to use the
#      : join() function. For those familiar with SQL, this may be a more
#      : intuitive way to understand the merging of two data.frames. An
#      : equivalent SQL statement might look something like this:
#      : SELECT *
#      : FROM dataProjected@data
#      : INNER JOIN watershedPoints
#      : ON dataProjected@data$id = watershedPoints$id

# library(plyr)
# watershedDF <- join(watershedPoints, dataProjected@data, by = "id")

Что означает весь этот код и зачем он нам нужен? Давайте пройдемся по этой строке за строкой.

dataProjected@data$id <- rownames(dataProjected@data)

Здесь мы добавляем к данным дополнительный столбец с именем «id». Этот столбец будет содержать имена строк, поэтому мы определим явную связь между данными и полигонами, связанными с этими данными.

watershedPoints <- fortify(dataProjected, region = "id")

Fortify? Что это вообще значит? Быстрый поиск в Интернете даст полезную документацию. (См. Документацию fortify.sp ). В основном, fortifyвозьмите два аргумента:, modelкоторый будет состоять из SpatialPolygonsDataFrameобъекта, который мы хотим преобразовать, и regionимя переменной, по которой нужно разделить регионы. Если все идет по плану, происходит какая-то магия, и мы получаем data.frame, как мы и хотели … ну, не совсем. Если вы проверите это data.frame, вы заметите, что в нем отсутствует какая-то важная информация. Не волнуйтесь! Используя отношения, которые мы создали ранее, мы можем объединить эти два набора данных с помощью следующей команды.

watershedDF <- merge(watershedPoints, dataProjected@data, by = "id")

И альт! Теперь, когда мы создали , data.frameчто ggplot2любит, мы можем начать построение. Прежде чем мы начнем рисовать, давайте взглянем на это новое, data.frameчто мы создали.

head(watershedDF)