我想根据一个县内疾病的频率绘制一张有 n 个正方形的地理地图。就像这里的这个:
但我不知道如何用 R 或 qGIS 做到这一点。谢谢大家的帮助。
首先编写一个名为“stackbox”的小函数,它在正确的位置绘制(使用“rect”)堆栈的正方形。
这是该函数的第一行:
stackbox <- function(x,y,n,size,maxheight=5){
然后为每个有病例的县调用该函数。
你到底在哪里卡住了这个过程?
这是完整的功能:
stackbox <- function(x,y,n,size,maxheight=5,...){
stackheight = seq(0,n,by=maxheight)
stackheight=diff(unique(c(stackheight,n)))
for(col in 1:length(stackheight)){
xl=rep(x+(col-1)*size,stackheight[col]) - (length(stackheight)/2)*size
yb=y+size*((1:stackheight[col])-1) - (max(stackheight)/2)*size
xr=xl+size
yt=yb+size
rect(xl,yb,xr,yt,...)
}
}
例子:
plot(1:10)
for(i in 1:10){
stackbox(i,i,i,3,size=.1,col="red",border="white")
}
要在地图上执行此操作,您需要 sp 和 maptools 包,以及其中包含您的数据的 shapefile 或其他地理空间数据源:
africa=readShapeSpatial(file.path(mapLib,"africa.shp"))
plot(africa,border="gray")
coords=coordinates(africa)
for(i in 1:nrow(africa)){
if(cases[i]>0){
stackbox(coords[i,1],coords[i,2],africa$cases[i],1,border="#606060",col="#0083FE")
}
}
我选择了看起来有点像你原来的颜色。请注意,这些框位于绘图坐标中,因此我必须将它们设为 0.1 度。您可能希望将您的地图转换为欧几里得投影(使用 package:gdal 中的 spTransform),然后它们将采用这些单位并适当地平方。
不过,像原来那样做漂亮的文本标签会很棘手......
作为替代方案,我建议使用气泡图,使用 ggplot2 使用点几何很容易制作。帮助页面的示例(http://had.co.nz/ggplot2/geom_point.html):
p <- ggplot(mtcars, aes(wt, mpg))
p + geom_point(aes(size = qsec))
这导致:
气泡的大小代表测量的量。该示例不使用地理坐标,但可以轻松使用它们。ggplot2 包还支持添加空间多边形作为附加层,例如使用 geom_path。coord_map的帮助文件显示了一些使用地图数据的示例。另请参见将 SpatialPolygons 对象转换为 data.frame 的 fortify 函数(ggplot2 需要该函数)。