R-TIGER Shapefile问题中的等值线图

Question

对R映射有一个问题,特别是在R中的等值区映射周围.

我有一个分配给一个和一些相关数据的邮政编码数据集(数据集在这里).

我的最终数据格式是:区域ID,ZIP,概率值,客户计数,区域概率和区域客户总数.我试图通过在地图上绘制面积概率和面积客户总数来呈现此数据.我试图通过使用人口普查TIGER Shapefiles来做到这一点,但我猜R不能处理整个国家.

我对统计功能感到满意,现在我将所有我的映射从第三方GIS专注的应用程序转移到我在R中完成所有映射.有没有人有任何关于如何在R内实现这一点的指针？

为了更加详细一点,这里是R停止工作的地方 -

shapes <- readShapeSpatial("tl_2013_us_zcta510.shp")

(其中shp文件是人口普查/ TIGER)形状文件.

编辑 - 提供更多详细信息.我试图首先阅读TIGER shapefile,希望将这个空间数据集与我的数据结合起来并最终绘制.我在尝试读取形状文件时最初遇到问题.下面是输出的代码

require(maptools)
shapes<-readShapeSpatial("tl_2013_us_zcta510.shp")

Error: cannot allocate vector of size 317 Kb

Answer 1

有几个关于使用R制作地图的示例和教程,但大多数是非常一般的,不幸的是,大多数地图项目都有细微差别,会产生难以理解的问题.你的就是一个很好的例子.

我遇到的最大问题是美国人口普查局的邮政编码区域形状文件对于整个美国来说是巨大的:~800MB.使用readOGR(...)R SpatialPolygonDataFrame对象加载时大约是913MB.尝试处理这个大小的文件(例如,使用转换为数据框fortify(...)),至少在我的系统上,导致了如上所述的错误.因此,解决方案是根据实际存储在数据中的邮政编码对文件进行子集化.

这张地图:

是使用以下代码从您的数据.

library(rgdal)
library(ggplot2)
library(stringr)
library(RColorBrewer)

setwd("<directory containing shapfiles and sample data>")

data     <- read.csv("Sample.csv",header=T) # your sample data, downloaded as csv
data$ZIP <- str_pad(data$ZIP,5,"left","0") # convert ZIP to char(5) w/leading zeros

zips     <- readOGR(dsn=".","tl_2013_us_zcta510") # import zip code polygon shapefile
map      <- zips[zips$ZCTA5CE10 %in% data$ZIP,]   # extract only zips in your Sample.csv
map.df   <- fortify(map)        # convert to data frame suitable for plotting
# merge data from Samples.csv into map data frame
map.data <- data.frame(id=rownames(map@data),ZIP=map@data$ZCTA5CE10)
map.data <- merge(map.data,data,by="ZIP")
map.df   <- merge(map.df,map.data,by="id")
# load state boundaries
states <- readOGR(dsn=".","gz_2010_us_040_00_5m")
states <- states[states$NAME %in% c("New York","New Jersey"),] # extract NY and NJ
states.df <- fortify(states)    # convert to data frame suitable for plotting

ggMap <- ggplot(data = map.df, aes(long, lat, group = group)) 
ggMap <- ggMap + geom_polygon(aes(fill = Probability_1))
ggMap <- ggMap + geom_path(data=states.df, aes(x=long,y=lat,group=group))
ggMap <- ggMap + scale_fill_gradientn(name="Probability",colours=brewer.pal(9,"Reds"))
ggMap <- ggMap + coord_equal()
ggMap

说明:

该rgdal软件包有助于从ESRI shapefile创建R Spatial对象.在您的情况下,我们将多边形shapefile导入R中的SpatialPolygonDataFrame对象.后者有两个主要部分:多边形部分,其中包含将在地图上创建多边形而连接的纬度和经度点,以及数据部分其中包含有关多边形的信息(因此,每个多边形一行).例如,如果我们调用Spatial对象map,那么这两个部分可以作为map@polygons和引用map@data.制作等值区域图的基本挑战是将Sample.csv文件中的数据与相关的多边形(邮政编码)相关联.

所以基本的工作流程如下:

1. Load polygon shapefiles into Spatial object ( => zips)
2. Subset if appropriate ( => map).
3. Convert to data frame suitable for plotting ( => map.df).
4. Merge data from Sample.csv into map.df.
5. Draw the map.

第4步是导致所有问题的那个.首先,我们必须将邮政编码与每个多边形相关联.然后我们必须Probability_1与每个邮政编码相关联.这是一个三步过程.

Spatial数据文件中的每个多边形都有唯一的ID,但这些ID 不是邮政编码.多边形ID存储为行名称map@data.邮政编码存储map@data在列中ZCTA5CE10.首先,我们必须创建一个将map@data行名(id)与map@data$ZCTA5CE10(ZIP)相关联的数据框.然后我们Sample.csv使用两个数据框中的ZIP字段将您的结果与结果合并.然后我们将其结果合并到map.df.这可以用3行代码完成.

绘制地图涉及告诉ggplot要使用的数据集(map.df),用于x和y(长和纬度)的列以及如何按多边形(group = group)对数据进行分组.列long,lat以及group在map.df都通过调用创建fortify(...).调用geom_polygon(...)告诉ggplot绘制多边形并使用中的信息填充map.df$Probability_1.调用geom_path(...)告诉ggplot创建一个具有状态边界的图层.调用scale_fill_gradientn(...)告诉ggplot使用基于颜色酿造者"Reds"调色板的配色方案.最后,调用coord_equal(...)告诉ggplot对x和y使用相同的比例,这样地图就不会失真.

NB:状态边界层,使用美国TIGER文件.