vlambda博客
学习文章列表

library(rgdal)# 读取文件v = read.csv('data.csv', 1, encoding = 'UTF-8')vlibrary(gstat)library(sp)data(meuse)head(meuse)coordinates(meuse) <- c("x","y")#定义坐标 class(meuse)spplot(meuse,"cadmium")plot(v,plot.number = T)v <- variogram(log(meuse$cadmium) ~ 1,meuse, width = 100)plot(v, plot.number = T)m1 <- vgm(psill = 1.5,model = "Sph",#Sph球形曲线 range = 1400,nugget = 0.5)m2 <- fit.variogram(v,m1)#拟合半变异函数"plot(v,model = m2)m2data(meuse.grid)data(meuse)head(meuse.grid)kr <- krige(log(meuse$cadmium) ~ 1, loc = ~ x+y, #位置信息 data = meuse,  newdata = meuse.grid,  model = m2)head(kr)library(lattice)levelplot(var1.pred ~ x + y,kr )

#R语言进行降水站点的克里金插值setwd("F:/Rpeng/21")## 1---加载程序包library(raster)library(sp)library(rgdal)library(gstat)library(raster)library(maptools)## 2---设置工作空间setwd("D:/CBR/2016Krige")## 3---读取站点csv数据#导入进行插值的点, 位置在工作空间内，#格式暂设为code，lon，lat,ST值,year,month,day(1~31),day1(1~366).data<-read.table("2016_ST.csv",header=TRUE,sep=",")#看一下表头head(data)#建立新列time，格式为2016_1_1data$time<-paste(data$year,data$month,data$day,sep='_')head(data)## 4---导入研究区矢量地图bound<-readOGR("CBH.shp")#plot(bound,col="grey") 用于显示矢量图的代码## 5---读取遥感tiff make插值模板#读取一个tifraster0 <- raster("D:/CBR/2016CCS_CUT/clip_CCS_1d20160621.tif")summary(raster0)#plot(raster0) 用于显示tif的代码blank_raster<-raster0 #新建空栅格values(blank_raster)<-0bound_raster<-rasterize(bound,blank_raster) #空栅格裁剪bound_raster[!(is.na(bound_raster))] <- 0#plot(bound_raster,main=paste("CBH"))#获取栅格的格网数据Grid<-as(bound_raster,"SpatialGridDataFrame")grd_LL<-rasterToPoints(raster0,spatial=F)grd_LL<-grd_LL[,-3] #第三列是数值0，进行删除#普通克里金插值#克里金法是通过一组具有 z 值的分散点生成估计表面的高级地统计过程。#与插值工具集中的其他插值方法不同，选择用于生成输出表面的最佳估算方法之前，#有效使用克里金法工具涉及 z 值表示的现象的空间行为的交互研究。## 6---进行克里金插值for(i in 1:366){# i<-1#获取某日数据dataone=NULLindex=which(data$day1==i)dataone=data[index,]time_i<-unique(dataone$time)st_judgement<-length(unique(dataone$ST))#当插值的点的非零值大于5时进行插值if (st_judgement>=5){##选取空间数据[,2:3]表示二三列的全部数据dsp <- SpatialPoints(dataone[,2:3], proj4string=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))dsp <- SpatialPointsDataFrame(dsp,dataone)#计算半方差#v<- variogram(log(ST) ~ 1, data =dsp)v<- variogram(ST ~ 1, data =dsp)plot(v)#不同的半方差模拟函数#v.fit<-fit.variogram(v,model=vgm(1,"Lin",0))#v.fit<-fit.variogram(v,model=vgm(model = "Exp", psill=20000, range = 1000, kappa = 10))v.fit<-fit.variogram(v,model=vgm(model = "Sph", psill=1.5, range = 6, kappa = 0.5),fit.sills = 4500)plot(v,v.fit)v.fit#location为已知点的坐标；#newdata为需要插值的点的位置；#nmax和nmin分别代表最多和最少搜索点的个数kri<-krige(formula=ST~1,model=v.fit,locations=dsp,newdata=Grid,nmax=12, nmin=8)#kr <- krige(formula = ST ~ 1, dataone,basexy,model =v.fit)spplot(kri["var1.pred"])writeRaster(raster(kri["var1.pred"]),paste("./","ST_K_",time_i,".tif",sep=""))write.csv(data.frame(grd_LL,kri$var1.pred),paste("./","ST_pred_",time_i,".csv",sep=""))write.csv(data.frame(grd_LL,kri$var1.var),paste("./","ST_var_",time_i,".csv",sep=""))} else#否则输出0文件{writeRaster(raster(blank_raster),paste("./","ST_0_",time_i,".tif",sep=""))write.csv(data.frame(grd_LL,blank_raster@data@values),paste("./","ST_0pred_",time_i,".csv",sep=""))write.csv(data.frame(grd_LL,blank_raster@data@values),paste("./","ST_0var_",time_i,".csv",sep=""))}}

#作者：北师大申泽西setwd("F:/Rpeng/21")library(rgdal)library(tmap)library(raster)library(ggplot2)library(maptools)library(maps)library(gstat)library(RColorBrewer)ex = extent(105,112,42,47)shp = as(ex,'SpatialPolygons')crs(shp) = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"long = runif(20,105,112)lat = runif(20,42,47)pr = runif(20,0,15)df = data.frame(long = long,lat = lat,pr = pr)head(df) long_shp = extent(shp)[1:2] lat_shp = extent(shp)[3:4]  long_p = seq(min(long_shp),max(long_shp),by =resolution )  lat_p = seq(min(lat_shp),max(lat_shp),by = resolution) s_p = expand.grid(long_p,lat_p) inusa = map.where(x = s_p[,1],y = s_p[,2]) inusa = !is.na(inusa) s_p = s_p[inusa,]  colnames(s_p) = c('long','lat') coordinates(s_p) = ~ long + lat gridded(s_p) = T fullgrid(s_p) = F proj4string(s_p) = crs(shp)pr = as.numeric(df[,3])  na_index = which(is.na(pr)) df = data.frame(long=df[,1],lat=df[,2],pr=as.numeric(df[,3])) df2 = df if(length(na_index)>0){ df = df[-na_index,]  } coordinates(df) = ~ long+lat proj4string(df) = crs(shp) print("Pass Phase 2")  f1 = as.formula(pr ~ long+lat) var.sampl = variogram(f1,df,cloud =F, cutoff = 100000, width = 89900) dat.fit = fit.variogram(var.sampl,fit.ranges = F, fit.sills = F, vgm(psill = 14,model = 'Sph', range = 590000,                              nugget = 0)) dat.krig = krige(f1,df,s_p,dat.fit) dat.r = raster(dat.krig) writeRaster(dat.r,'raster.tif',format = 'GTiff',overwrite = T)krg.res = as.data.frame(dat.r,xy = T)  colnames(krg.res) = c('long','lat','pr')  mycolor = colorRampPalette(brewer.pal(11,'Spectral'))(12) p = ggplot()+ geom_tile(data = krg.res,aes(x = long,y = lat,fill = pr), color = 'transparent')+ geom_point(data = df2,aes(x = long,y = lat),color = 'black', shape =1, size = 5)+ scale_fill_gradientn(colours = mycolor)+ theme_bw()+ theme( axis.text = element_text(face = 'bold',color = 'black',size = 12), axis.title = element_text(face = 'bold',color = 'black',size = 14, hjust = .5), legend.text = element_text(face = 'bold',color = 'black',size = 12), legend.title = element_blank(), legend.position = 'bottom', legend.direction = 'horizontal', legend.key.width = unit(2,'cm') )+ xlab('Longitude')+ ylab('Latitude') p  png('krige.png', height = 20, width = 20, units = 'cm', res = 800) print(p) dev.off()

library(ggplot2)data("USArrests") # 美国犯罪率数据集crimes <- data.frame(state = tolower(rownames(USArrests)), USArrests)crimes$tempvar <- "Kingjames Kingpeng"crimesm <- reshape2::melt(crimes, id = 1) # 宽转长# 绘制地图states_map <- map_data("state") # ggplot2内置州级地图数据集USA_map <-ggplot(crimes, aes(map_id = state)) +geom_map(aes(fill=Murder), map=states_map,colour="cyan",size = 0.5) + # 指定map参数scale_fill_distiller(palette = "RdPu",direction= 1) + # 修改颜色标度expand_limits(x = states_map$long, y = states_map$lat) # 扩大显示，显示所有在经纬度内的USA_mapUSA_map+ facet_grid(. ~ tempvar) + theme(strip.background = element_rect(fill="#999999"), strip.text = element_text(size=15, colour="white"))

library(ggplot2)library(sf)library(magrittr)library(rgdal)library(dplyr)#读取各省中心点坐标数据path1<-"F:/Rpeng/21/prov_centroids.csv"Prov_centers<-read.csv(path1,stringsAsFactors=FALSE)str(Prov_centers)#读取地图坐标数据path2<-"F:/Rpeng/21/china6.json"China_1<-st_read(path2,stringsAsFactors=FALSE)#编造业务数据mydata_1<-data.frame(id=1:34,name=Prov_centers$name,scale=runif(34,100,200)%>%round(),Scope=rep(LETTERS[1:5],length=34))#合并数据集mydata_2<-left_join(Prov_centers,mydata_1,by="name")mydata_3<-left_join(China_1,mydata_2,by="name")mydata_3$tempvar <- "Kingjames Kingpeng"#连续变量ggplot(mydata_3)+geom_sf(aes(fill=scale),color="cyan")+coord_sf()+scale_fill_continuous(type="viridis")#调整颜色标度###ggplot(mydata_3)+geom_sf(aes(fill=scale),color="cyan")+coord_sf()+scale_fill_continuous(type="viridis")+ facet_grid(. ~ tempvar) + theme(strip.background = element_rect(fill="#999999"), strip.text = element_text(size=15, colour="white"))

setwd("F:/Rpeng/21")library(ggplot2)#创建数据集df <- data.frame(treatment = factor(c(1, 1, 1, 2, 2, 2, 3, 3, 3)),response = c(2, 5, 4, 6, 9, 7, 3, 5, 8),group = factor(c(1, 2, 3, 1, 2, 3, 1, 2, 3)),se = c(0.4, 0.2, 0.4, 0.5, 0.3, 0.2, 0.4, 0.6, 0.7))head(df) #查看数据集# 使用geom_errorbar()绘制带有误差棒的条形图# 这里一定要注意position要与`geom_bar()`保持一致，由于系统默认dodge是0.9，# 因此geom_errorbar()里面position需要设置0.9，width设置误差棒的大小df$tempvar <- "Kingjames Kingpeng"ggplot(data = df, aes(x = treatment, y = response, fill = group)) +geom_bar(stat = "identity", position = "dodge") +geom_errorbar(aes(ymax = response + se, ymin = response - se),position = position_dodge(0.9), width = 0.15) +scale_fill_brewer(palette = "Set1")+ facet_grid(. ~ tempvar) + theme(strip.background = element_rect(fill="#999999"), strip.text = element_text(size=15, colour="white"))