首先，我想展示如何使用R制作简单的地图。此外，我想展示包装中的一个重要功能“gdalUtils”，其中包含地理空间数据抽象库（GDAL）实用程序的包装。我将使用gdal_rasterize函数对向量进行栅格化，因为对于大型向量数据集来说它非常快。

library(ggplot2)
library(maptools)
library(rgdal)
library(raster)
library(plyr)
library(gdalUtils)
library(RColorBrewer)
library(rasterVis)

加载世界数据（包装maptools）并对其进行绘制以获得全景。

数据(wrld_simpl)
plot(wrld_simpl)

子集数据。选择奥地利国家（aut_poly）&奥地利，其邻居在西部（aut_neigh）。

aut<-wrld_simpl[[email protected]$NAME=="Austria",]
aut_neigh<-wrld_simpl[[email protected]$SUBREGION==155,]

Ggplot仅支持data.frames，这就是我们将空间多边形类转换为数据框的原因。

aut_gg<-fortify(aut)
aut_neigh_gg<-fortify(aut_neigh)
p1<-ggplot()+geom_map(data=aut_neigh_gg,map=aut_neigh_gg,aes(x=long,y=lat,map_id=id),fill="white",color="black", size=0.25)+geom_map(data=aut_gg,map=aut_gg,aes(map_id=id,fill=id))
p1

现在，我们更改地图的主题并添加标签。标签需要一个向量来定义
位置。因此，我们计算平均值X坐标和Y坐标，并将它们保存在data.frame中“labs.out_neigh”.
p2 向我们显示了一个新地图，其中包含每个国家的标题和标签。如果您想为标签找到更好的位置，可以扩展功能 geom_text 带有更多参数：一些参数是角度，垂直，尺寸….

labs.aut_neigh<-ddply(aut_neigh_gg,.(id),summarise, meanx = mean(range(long)),meany=mean(range(lat)))
p2<-p1+ggtitle("Part of Europe")+geom_text(data=labs.aut_neigh,aes(meanx,meany,label=id))
p3<-p2+theme_bw()+theme(legend.position = "right")+
scale_fill_discrete(name="Country-Id",l=50)
p3

我经常面临将矢量数据（形状文件）转换为栅格文件的任务。有几种解决方案的方法。我准备了两个例子。

解决方法一：

extent(aut_neigh) # Delivers us the extent of our area of interest.
raster_aut_neigh_1<-raster(extent(aut_neigh)) # Creates an 栅格 based on the extent
raster_aut_neigh_1 # Shows use the number of rows and columns. Used in the following row
raster_aut_neigh_1[] = runif(10*10) # Fills the 栅格s with randomly selected values. 
plot(raster_aut_neigh_1) # Plot 栅格
plot(aut_neigh,add=TRUE) # Add polygon to the plot

res(raster_aut_neigh_1)=c(0.5,0.5) # We increase the spatial resolution
raster_aut_neigh_1[] = runif(28*44) 
plot(raster_aut_neigh_1)
plot(aut_neigh,add=TRUE)

栅格_aut_neigh_final<-rasterize(aut_neigh,raster_aut_neigh_1,field="ISO3")
proj4string(raster_aut_neigh_final)[email protected]
plot(raster_aut_neigh_final)

解决方案二：

让’假设我们已经有一个形状文件。所以我们用 盖达尔Utils软件包 并将形状文件转换为栅格文件。首先，我们创建感兴趣区域的形状文件。

diroutput<-"C:/Users/YourName/Documents/Foo/World_map" #define your path
# Writes a shape file into your directory
writeOGR(aut_neigh,dsn=diroutput,layer="aut_neigh",driver="ESRI Shapefile") 
# Load shape file in R
shape_file<-readOGR(dsn=diroutput,"aut_neigh")

现在我们变换形状文件shape_file）放入栅格文件。

#find EPSG code for the Raster
EPSG <- make_EPSG() # Functions to create a 数据.frame which contains all EPSG-Codes
[email protected]@projargs # obtains the projection
# Looking after the EPSG-Codes: you will get one row with EPSG: 4326
EPSG[grepl("+proj=longlat", EPSG$prj4, fixed=TRUE) & grepl("+datum=WGS84", EPSG$prj4, fixed=TRUE), ]
#> 249 4326 # WGS 84 +proj=longlat +datum=WGS84 +no_defs

最后，我们为最后的计算定义了两条路径。请在我的系统属性中修改我的代码。

inputdir<-"C:/Users/Foo/Documents/Ba/World_map"
outputdir<-"C:/Users/Foo/Documents/Ba/World_map/gdal_shp_2_r.tif"

#tr= target layer resolution; a_srs=EPSG-Code; l= shape file layer name; 
#a = shape file attribute field to copy
r_id<-gdal_rasterize(inputdir, outputdir,tr=c(0.5,0.5),a_nodata=-9999,a_srs="epsg:4326",l="aut_neigh",a="UN",output_Raster = FALSE)
#load 栅格化d polygon
r_id2<-raster(outputdir)
# define 栅格 field attribute name
names(r_id2)<-"ID"

现在我们绘制最终结果。

# We take a look to the relationship between country and the 栅格 cell value "UN"
sort_aut<[email protected][order([email protected]$UN),]

使用以下命令，我们将数字栅格调用值重新分类。
图例显示了与数值属性相关的类别，而不是数字单元格值“UN”

r_test <- ratify(r_id2) 
rat <- levels(r_test)[[1]]
rat$country <- c(as.character(sort_aut[sort_aut$UN%in%rat[,1],"ISO3"]))
levels(r_test)<-rat
my.color<-brewer.pal(7,"Dark2") # define your own color palette
levelplot(r_test,col.regions=my.color)