在不规则的网格上绘制轮廓
我已经阅读了R中的等高线图(包括很多关于stackoverflow的提示)的页面和页面,但都没有成功。 这里是我的数据轮廓,包括添加一个卢旺达的地图(数据包括经度,纬度和雨的十四个值,如x,y和z):
Lon Lat Rain 28.92 -2.47 83.4 29.02 -2.68 144 29.25 -1.67 134.7 29.42 -2.07 174.9 29.55 -1.58 151.5 29.57 -2.48 224.1 29.6 -1.5 254.3 29.72 -2.18 173.9 30.03 -1.95 154.8 30.05 -1.6 152.2 30.13 -1.97 126.2 30.33 -1.3 98.5 30.45 -1.81 145.5 30.5 -2.15 151.3
这里是我从stackoverflow尝试的代码:
datr <- read.table("Apr0130precip.txt",header=TRUE,sep=",") x <- datr$x y <- datr$y z <- datr$z require(akima) fld <- interp(x,y,z) par(mar=c(5,5,1,1)) filled.contour(fld)
插值fails.help将不胜感激。
这里有一些不同的使用base
Rgraphics和ggplot
可能性。 这两个简单的轮廓图,生成地图顶部的地块。
插值
library(akima) fld <- with(df, interp(x = Lon, y = Lat, z = Rain))
使用filled.contour
绘制base
R图
filled.contour(x = fld$x, y = fld$y, z = fld$z, color.palette = colorRampPalette(c("white", "blue")), xlab = "Longitude", ylab = "Latitude", main = "Rwandan rainfall", key.title = title(main = "Rain (mm)", cex.main = 1))
基本的ggplot
替代使用geom_tile
和stat_contour
library(ggplot2) library(reshape2) # prepare data in long format df <- melt(fld$z, na.rm = TRUE) names(df) <- c("x", "y", "Rain") df$Lon <- fld$x[df$x] df$Lat <- fld$y[df$y] ggplot(data = df, aes(x = Lon, y = Lat, z = Rain)) + geom_tile(aes(fill = Rain)) + stat_contour() + ggtitle("Rwandan rainfall") + xlab("Longitude") + ylab("Latitude") + scale_fill_continuous(name = "Rain (mm)", low = "white", high = "blue") + theme(plot.title = element_text(size = 25, face = "bold"), legend.title = element_text(size = 15), axis.text = element_text(size = 15), axis.title.x = element_text(size = 20, vjust = -0.5), axis.title.y = element_text(size = 20, vjust = 0.2), legend.text = element_text(size = 10))
由ggmap
创build的Google地图上的ggmap
# grab a map. get_map creates a raster object library(ggmap) rwanda1 <- get_map(location = c(lon = 29.75, lat = -2), zoom = 9, maptype = "toner", source = "stamen") # alternative map # rwanda2 <- get_map(location = c(lon = 29.75, lat = -2), # zoom = 9, # maptype = "terrain") # plot the raster map g1 <- ggmap(rwanda1) g1 # plot map and rain data # use coord_map with default mercator projection g1 + geom_tile(data = df, aes(x = Lon, y = Lat, z = Rain, fill = Rain), alpha = 0.8) + stat_contour(data = df, aes(x = Lon, y = Lat, z = Rain)) + ggtitle("Rwandan rainfall") + xlab("Longitude") + ylab("Latitude") + scale_fill_continuous(name = "Rain (mm)", low = "white", high = "blue") + theme(plot.title = element_text(size = 25, face = "bold"), legend.title = element_text(size = 15), axis.text = element_text(size = 15), axis.title.x = element_text(size = 20, vjust = -0.5), axis.title.y = element_text(size = 20, vjust = 0.2), legend.text = element_text(size = 10)) + coord_map()
ggplot
从shapefile创build的地图上
# Since I don't have your map object, I do like this instead: # get map data from # http://biogeo.ucdavis.edu/data/diva/adm/RWA_adm.zip # unzip files to folder named "rwanda" # read shapefile with rgdal::readOGR # just try the first out of three shapefiles, which seemed to work. # 'dsn' (data source name) is the folder where the shapefile is located # 'layer' is the name of the shapefile without the .shp extension. library(rgdal) rwa <- readOGR(dsn = "rwanda", layer = "RWA_adm0") class(rwa) # [1] "SpatialPolygonsDataFrame" # convert SpatialPolygonsDataFrame object to data.frame rwa2 <- fortify(rwa) class(rwa2) # [1] "data.frame" # plot map and raindata ggplot() + geom_polygon(data = rwa2, aes(x = long, y = lat, group = group), colour = "black", size = 0.5, fill = "white") + geom_tile(data = df, aes(x = Lon, y = Lat, z = Rain, fill = Rain), alpha = 0.8) + stat_contour(data = df, aes(x = Lon, y = Lat, z = Rain)) + ggtitle("Rwandan rainfall") + xlab("Longitude") + ylab("Latitude") + scale_fill_continuous(name = "Rain (mm)", low = "white", high = "blue") + theme_bw() + theme(plot.title = element_text(size = 25, face = "bold"), legend.title = element_text(size = 15), axis.text = element_text(size = 15), axis.title.x = element_text(size = 20, vjust = -0.5), axis.title.y = element_text(size = 20, vjust = 0.2), legend.text = element_text(size = 10)) + coord_map()
降雨量数据的插值和绘图当然可以用复杂得多的方法完成,使用R中好的空间数据工具 。 考虑我的答案是一个相当快速和简单的开始。