python 克里金插值_R语言克里金插值（三）

事实证明： 自拍后， 只要你用心P图3小时，大家就看不出来你是 肥仔 ！！！

‍

setwd("F:/Rpeng/23")#### 首先生成一套模拟数据simdat, 通过Volcano数据随机提取1000个点, 用raster的extract函数提取library(raster)rrr 0, xmx= x11()plot(rrr, col = terrain.colors(100))set.seed(12345)loc x = round(runif(value simdat ### The simulated data#### Krig准备工作开始, 导入程序包library(sp)library(gstat)#### simdat需要转换为sp格式, 需要指定坐标coordinates(simdat) x+### bubble(simdat, "alt")spplot(simdat, "alt") ### 查看每个点的值 #### 创建克里格所用的网格, 20m方格的中心位置xgrid q(ceiling(min(simdat$x)), floor(max(simdat$x)), by = ygrid q(ceiling(min(simdat$y)), floor(max(simdat$y)), by = basexy #### plot(y ~ x, basexy)colnames(basexy) "x", coordinates(basexy) x+gridded(basexy) vgm1 1, simdat)plot(vgm1, plot.numbers = TRUE) m 01,plot(vgm1, model=m)### 设定vgm模型的初始值m 01,### 进行克里格krige_res 1, simdat, basexy, model = ### 查看克里格插值的结果spplot(krige_res, zcol = "var1.pred", main = "Predictions of altitude based on the randomly sampled data", col.regions = terrain.colors(100))

setwd("F:/Rpeng/23")library(maps)      ## 提供绘制世界地图的函数library(mapdata)   ## 提供世界地图library(rgdal)library(ggplot2)   ## 绘图核心函数library(maptools)  ## 提供 readShapePoly 函数， 读取中国多边形数据china0    "F:/Rpeng/23/shp/bou1_4p") china_df  ## 将list转换为ggplot可以使用的dataframe## 读取数据 https://simplemaps.com/static/data/world-cities/basic/simplemaps-worldcities-basic.csvcities "F:/Rpeng/23/shp/simplemaps-worldcities-basic.csv")cities 1:nrow(cities), size = pdf(file = "sampling sites location.pdf", width = 11.5, height = 5.8)cities$tempvar "Kingjames Kingpeng"rrr  borders(colour='darkgrey') +                  ## 添加世界地图 geom_polygon(data = china_df,   aes(x = long, y = lat, group = group), colour = "darkgrey",   fill = "white") +                           ## 添加中国地图 geom_point(data = cities,   aes(x=lng, y=lat, colour=pop), size = 2) +  ##点的颜色代表人口   theme_bw()+scale_colour_gradientn(colours=c("green", "yellow", "orange", "red"))+facet_grid(. ~ tempvar) +        theme(strip.background = element_rect(fill="#999999"),              strip.text = element_text(size=15, colour="white"))rrrdev.off()

setwd("F:/Rpeng/23/sample_map")library(tmap)library(tmaptools)library(sp)library(rgdal)rm(list = ls())## 读取地图country "bou1_4l.shp")province "province_polygon.shp")world "ne_50m_admin_0_countries.shp")# Source: https://www.naturalearthdata.com/downloads/50m-cultural-vectors/ world "longlat")    country "longlat")    province "longlat")    province "longlat")    data(World)# 城市city "simplemaps-worldcities-basic.csv", header = city 1:nrow(city), coordinates(city) proj4string(city) "+proj=longlat +datum=WGS84")World$name[World$name == "Taiwan"] # 转换为spatial dataframe# 绘图tm_shape(world,          xlim = c(60, 140),          ylim = c(0, 60)) +   tm_borders("grey40", lwd = 1.5) +   tm_grid()+ tm_shape(province) +  tm_fill(col = "lightgrey") +   tm_borders("grey60", lwd = 0.8) + tm_shape(country) +   tm_lines(col="grey40", lwd = 1.5) + tm_scale_bar(position=c(0.05,0.0)) + tm_compass(type = "4star",            position=c("left", "top")) +     tm_layout(inner.margins=c(0.12,0.03,0.08,0.03),           legend.stack = "horizontal") + tm_shape(city) +    tm_bubbles("pop",                col = "red",                scale=.5,               border.col = "red") +     tm_text("city", size=0.5) + tm_legend(legend.position = c(0.05, 0.08)) + tm_shape(World) +    tm_text("name", size=1.0)

setwd("F:/Rpeng/23")library(ggplot2)library(rgdal)rm(list = ls())  ## 删除所有对象## 去过的地方provinces_visited "Beijing",                "Liaoning", "Jilin", "Jiangxi", "Hunan",               "Zhejiang", "Hubei", "Chongqing", "Sichuan",               "Henan", "Guangxi", "Guangdong", "Fujian",               "Hainan", "Yunnan", "Guizhou")provinces_visited_df                                    visited_status = rep("Visited", length(provinces_visited)))## 读取地图country      "F:/Rpeng/23/chinamap/bou1_4l.shp")province     "F:/Rpeng/23/chinamap/province_polygon.shp")    province$ID  as.character(province$ID) ## 转换为ggplot2绘图用的data.framecountry_df  province_df ## 提取省级数据中的信息province_dat ## 改正省级数据中关于香港和澳门的错误index_HK "Xianggang")province_dat$X[index_HK] 114.156121province_dat$Y[index_HK] 22.37725province_dat$ID[index_HK] "Hong Kong"index_MC "Aomen")province_dat$X[index_MC] 113.545681province_dat$Y[index_MC] 20.197303province_dat$ID[index_MC] "Macau"## id与province_df数据表中的省出现的顺序相同， 这里id都从0开始， 所以需要减去1，以便匹配province_dat$id as.character(## 增加省拼音province_dat2 "ID",                        by.y = "province_pinyin", all.x = TRUE)## 由于visited status 默认是 factor 类型， 这里需要先转换为字符串， 再处理province_dat2$visited_status as.character(province_dat2$visited_status)## 所有没有匹配上的省份都是没有去过的， 状态为 not yetprovince_dat2$visited_status[is.na(province_dat2$visited_status)] "Not Yet"province_df_merged "id")                   ## 绘图## 用geom_path绘制读取的polyline，因为国界线是polyline## 多圆锥投影 ## http://desktop.arcgis.com/zh-cn/arcmap/10.3/guide-books/map-projections/polyconic.htmcountry_df$tempvar "Kingjames Kingpeng"ggplot() + geom_path(data = country_df,              aes(x = long, y = lat, group=group), colour="darkblue") +  coord_map("polyconic") +                                                geom_polygon(data = province_df_merged,              aes(x = long, y = lat, group=group,                  fill = visited_status), colour="grey") + geom_text(mapping = aes(x = X, y = Y, label = ID),           data = province_dat, colour = 'Black')+xlab("Longitude") + ylab("Latitude") + theme(legend.position = "bottom") + scale_fill_discrete(name = "Province")+facet_grid(. ~ tempvar) +        theme(strip.background = element_rect(fill="#999999"),              strip.text = element_text(size=15, colour="white"))

setwd("F:/Rpeng/23")library(maps)library(ggplot2)thismap = map_data("world")thismap$tempvar "Kingjames Kingpeng"ggplot(thismap, aes(long, lat, group=group)) +  geom_polygon(fill="white", colour="gray") +facet_grid(. ~ tempvar) +        theme(strip.background = element_rect(fill="#999999"),              strip.text = element_text(size=15, colour="white"))ggplot(thismap, aes(long, lat, group=group, fill=region)) +   geom_polygon(show.legend = F) +facet_grid(. ~ tempvar) +        theme(strip.background = element_rect(fill="#999999"),              strip.text = element_text(size=15, colour="white"))library(mapdata) # with china map# map("china")ggplot(map("china", plot=F), aes(long, lat, group=group, fill=region)) +   geom_path(show.legend = F) +  ggtitle("Map of China")thismap "state")thismap$tempvar "Kingjames Kingpeng"ggplot(thismap, aes(long, lat, group=group, fill=region)) +   geom_polygon(show.legend = F)+facet_grid(. ~ tempvar) +        theme(strip.background = element_rect(fill="#999999"),              strip.text = element_text(size=15, colour="white"))