克里金插值详细步骤_R语言克里金插值（7）

气象水文科研猫公众号交流邮箱：leolovehydrometeor@hotmail.com欢迎投稿&批评指正如有侵权且本公众号未能正确引用原文，请联系删除，谢谢理解、谢谢配合。

R安装包：https://cran.r-project.orgRStudio安装包：https://rstudio.com/products/rstudio/R-ggplot2-cran：https://cran.r-project.org/web/packages/ggplot2/

Matlab之colorbar：

clc;clear;close alldata_xls=xlsread('testdata.xlsx');x=data_xls(:,1);y=data_xls(:,2);[X,Y] = meshgrid(x,y);Z = sin(X)+cos(Y);contour(X,Y,Z)colorbarcolormap('hsv');colorbar('ticklength',0) box offax2 = axes('Position',get(gca,'Position'),'XAxisLocation','top','YAxisLocation','right','Color','none','XColor','k','YColor','k');set(ax2,'YTick', []);set(ax2,'XTick', []);box on

R语言处理CMIP5：

rm(list = ls())library(stars)      # To process the raster datalibrary(sf)         # To work with vector datalibrary(ggplot2)    # For plottinglibrary(patchwork)  # To combine different ggplot plotsfile_path raster::getData(name = 'worldclim', var = 'bio', res = 10,path = file_path)raster::getData(name = 'CMIP5', var = 'bio', res = 10,rcp = 45, model = 'IP', year = 70,path = file_path)annual_T annual_T annual_T_70 annual_T_70 # The result, temp_colors, is a function with argument n for the number of# colors.temp_colors nbreaks {
      par(mfrow = c(1,2))  plot(annual_T, main = "Annual temperature - 1960-1990",     nbreaks = nbreaks,     col = temp_colors(nbreaks - 1))  plot(annual_T_70, main = "Annual temperature - RCP 4.5 projection for 2061-2080",     nbreaks = nbreaks,     col = temp_colors(nbreaks - 1))}south_america_map # The precision has to be set to a value > 0 to resolve internal boundaries.st_precision(south_america_map) south_america_map {
    par(mar = c(0,0,0,0))plot(south_america_map)}annual_T_SA # CRS for projected T and south america map are the same (EPSG 4326) but the# proj4string includes more details for annual_T_70. Thus, they have to# be made identical before cropping.st_crs(annual_T_70) annual_T_70_SA {
      par(mfrow = c(1, 2))  plot(annual_T_SA, main = "Annual temperature - 1960-1990",       nbreaks = nbreaks,       col = temp_colors(nbreaks - 1))  plot(main = "Annual temperature - RCP 4.5 projection for 2061-2080",       annual_T_70_SA, nbreaks = nbreaks,       col = temp_colors(nbreaks - 1))}annual_T_SAmean(annual_T$bio1.bil, na.rm = T)names(annual_T_SA) annual_T_SA$projected annual_T_SA$change annual_T_SArecent_T_plot   geom_stars(data = annual_T_SA) +  scale_fill_gradientn(name = "Annual T [°C]",                       colors = temp_colors(5),                       limits = c(-7, 32),                       na.value = "white") +  geom_sf(data = south_america_map, fill = NA) +  scale_x_discrete(expand = c(0, 0)) +  scale_y_discrete(expand = c(0, 0)) +  ggtitle("a) 1960-1990") +  theme_void() +  theme(legend.position = "none")projected_T_plot   geom_stars(data = annual_T_SA["projected"]) +  scale_fill_gradientn(name = "Annual T [°C]",                       colors = temp_colors(5),                       limits = c(-7, 32),                       na.value = "white") +  geom_sf(data = south_america_map, fill = NA) +  scale_x_discrete(expand = c(0, 0)) +  scale_y_discrete(expand = c(0, 0)) +  ggtitle("b) 2061-2080 (projected)") +  theme_void() +  theme(legend.position = "bottom")projected_change_T_plot   geom_stars(data = annual_T_SA["change"]) +  scale_fill_gradientn(name = "Change in T [°C]",                       colors = temp_colors(5)[3:5],                       limits = c(1, 5),                       na.value = "white") +  geom_sf(data = south_america_map, fill = NA) +  scale_x_discrete(expand = c(0, 0)) +  scale_y_discrete(expand = c(0, 0)) +  ggtitle("c) Projected change") +  theme_void() +  theme(legend.position = "bottom") (recent_T_plot / projected_T_plot + plot_layout(guides = "keep")) | projected_change_T_plot +  theme(plot.margin = margin(c(0, 0, 0, 0)))

---

R语言处理ERA5—Reanalysis：

#使用欧洲气候中心CDS Toolbox#https://cds.climate.copernicus.eu/#!/home#https://cds.climate.copernicus.eu/cdsapp#!/search?type=dataset#https://cds.climate.copernicus.eu/cdsapp#!/toolbox#https://cds.climate.copernicus.eu/toolbox-editor/60789/rea2##################################################################import cdstoolbox as ct@ct.application()@ct.output.download()def download_application():    data = ct.catalogue.retrieve(        'reanalysis-era5-land-monthly-means',        {
                'product_type': 'monthly_averaged_reanalysis',            'variable': 'total_precipitation',            'year': '2003',            'month': '06',            'time': '00:00',        }    )    return data

library(raster)setwd("F:/Rpeng/30/map_data") dset #dset dsetprint(dset)plot(dset)dset_r dset_rplot(dset_r)library(ggplot2)#ggplot(data = , mapping = aes()) +  ()df str(df)ggplot() +  geom_raster(data = df , aes(x = x, y = y, fill = Total.precipitation.rate)) +   coord_quickmap()ggplot() +  geom_raster(data = df , aes(x = x, y = y, fill = Total.precipitation.rate)) +  scale_fill_viridis_c(limits = c(0.0, 0.02)) +  coord_quickmap()# define jet colormapjet.colors # use the jet colormapggplot() +  geom_raster(data = df, aes(x=x, y=y, fill=Total.precipitation.rate)) +   scale_fill_gradientn(colors = jet.colors(7), limits = c(0.0, 0.000005)) +   coord_quickmap()ggplot() +  geom_raster(data = df, aes(x=x, y=y, fill=Total.precipitation.rate)) +   scale_fill_gradientn(colors = jet.colors(7), limits = c(0.0, 0.000005)) +   borders() +   coord_quickmap()ggplot(df, aes(y=y, x=x, color=Total.precipitation.rate)) +  geom_point(size=2, shape=15) +  borders('world', xlim=range(df$x), ylim=range(df$y), colour='black') +    scale_color_distiller(palette='Spectral') +  coord_map('ortho', orientation = c(40, 20, 0))+theme_void()

library(ggplot2)library(plyr)library("maptools")china_map = readShapePoly("F:/Rpeng/29/data/bou2_4p.shp")x xs china_map1 china_map_data mydata china_data ggplot(china_data, aes(x = long, y = lat, group = group,fill = zhibiao)) +     geom_polygon(colour="grey40") +     scale_fill_gradient(low="white",high="steelblue") +  #指定渐变填充色，可使用RGB     coord_map("polyconic")+theme_bw()province_city ggplot(china_data,aes(long,lat))+geom_polygon(aes(group=group,fill=zhibiao),colour="grey60")+     scale_fill_gradient(low="white",high="steelblue") +coord_map("polyconic") + geom_text(aes(x = jd,y = wd,label = province), data =province_city)+theme_bw()

library(ggplot2)library(plyr)library("maptools")CHN_adm2 CHN_adm2_1 x xs china_map_data anhui mydata=read.csv("F:/Rpeng/29/data/anhui.csv")                  #打开自定义的指标文件。anhui_data ggplot(anhui_data, aes(x = long, y = lat, group = group,fill=zhibiao)) +     geom_polygon(colour="grey40") +     scale_fill_gradient(low="white",high="steelblue") +     coord_map("polyconic") +theme_bw()

library(ggplot2)       ###绘图函数library(plyr)          ###数据合并工具library(maptools)      ###地图素材导入library(sp)#library(Cairo)         #图片高清导出library(RColorBrewer)  ###有一些高质量的地图配色模板可以参考CHN_adm2 CHN_adm2_1 data1 data2 china_map_data dongsansheng aaaaaadongsansheng_permdongsansheng_perm$zhibiao#使用write函数将东三省的行政单位信息导出到数据文件目录下，并取名dongsanshengwrite.table (dongsansheng_perm, file ="F:/Rpeng/29/data/CHN_adm/dongsansheng.csv", sep =",", row.names =FALSE)#使用read.csv函数将刚才整理好的带有业务数据的文件再次导入：mydata#将业务数据与地理信息数据再次合并dongsansheng_map_data ggplot(dongsansheng_map_data, aes(x = long, y = lat, group = group,fill=zhibiao)) +     geom_polygon(colour="black") +scale_fill_gradient(low="blue",high="red") +     coord_map("polyconic") +theme_bw()