r 语言 处理栅格数据

r 语言 处理栅格数据

作为新手,记录自己学到的一些东西,也希望能对需要的人有稍许帮助

常用的包

  • raster #处理栅格数据
  • rasterVis #栅格数据可视化
  • RColorBrewer #颜色设置
  • ncdf4 #nc数据的读取

常用的数据类型

  • tiff
  • nc

数据导入和简单的可视化

  • tiff数据读取

    • 读取一个tiff文件 raster() 函数

      ## 导入使用的包
      library(raster)
      library(rasterVis)
      library(RColorBrewer)
      library(ncdf4)
      
      n_data_path <- dir('数据存放地址',full.names = TRUE)
      n_deposition <- raster(n_data_path[19])
      plot(n_deposition) # 使用基础的plot函数
      levelplot(n_deposition, margin = list(FUN = 'median'),
          par.settings=RdBuTheme) # 使用levelplot函数
      

      使用raster::plot()
      使用levelplot()

    • 读取多个tiff文件 stack() 函数

      n_depo_all <- stack(n_data_path[c(1,6,11,16)])
      plot(n_depo_all) # 使用基础的plot函数
      levelplot(n_depo_all)  # 使用levelplot函数
      

    在这里插入图片描述
    在这里插入图片描述

  • nc数据读取
    同样使用 raster()stack() 函数

    sif_nc <- nc_open(dir('数据存放地址',full.names = TRUE)[2]) # 读取nc文件
    names(sif_nc$var)  # 获取nc文件中的变量名
    
    ## 读取一个时间点的数据(假设z轴为时间)
    sif_raster <-  raster(dir('数据存放地址',full.names = TRUE)[2], 
          varname = 'all_daily_sif',band=45)
    plot(sif_raster)
    
    ## 读取所有时间点的数据
    sif_raster_5 <-  stack(dir('数据存放地址',full.names = TRUE)[2], 
          varname = 'all_daily_sif',bands=40:43)
    
    

栅格数据的简单计算

  • 计算所有年份的氮沉降均值 (还可以使用 sum,min,max等函数)
    n_depo_ave <- mean(n_depo_all,na.rm=TRUE)
    plot(n_depo_ave)
    

在这里插入图片描述

可视化进阶

  • 颜色的设置

    plot(n_deposition,xlim=c(-130,-60),ylim=c(25,50),      
          col=colorRampPalette(c("blue","white","red"))(255))
    
    

在这里插入图片描述

  • 色卡的设置
    设置色卡的间断点
    建立色卡
    # 由于原始值分布不均,因此先将大于20的值设置为22
    n_depo_ave_reclassify <- n_depo_ave
    n_depo_ave_reclassify[n_depo_ave_reclassify>20] <- 22
    cutpts <- c(0,2,4,6,8,10,12,14,16,18,20,22)
    my_labels <- paste(cutpts[1:10],cutpts[2:11],sep = '-')
    myColorkey <- list(at=cutpts, ## where the colors change
                     labels=list(
                         labels=c(my_labels,'>20',''), ## labels
                         at=cutpts+1 ## where to print labels
                     ))
    levelplot(n_depo_ave_reclassify,margin=FALSE,colorkey=myColorkey,
            col.regions=(colorRampPalette(c('#f3e5f5','#9c27b0',"#6a1b9a"))))
    

在这里插入图片描述

  • 在栅格数据上叠加边界 使用 rgdal::readOGR() 函数读取shp
    us_border <- rgdal::readOGR('数据存放地址')
    p <- levelplot(n_depo_all, layers = 1, margin = list(FUN = 'median'))
    p+layer(sp.lines(us_border,lwd=0.8,col='black'))
    

在这里插入图片描述

栅格数据的剪切和重采样

读取的物候数据与氮沉降数据的分辨率和范围都不一样

  1. 使用 resample() 函数以氮沉降数据为标准对物候数据进行重采样
  2. 使用 mask() 截取氮沉降数据范围内的物候数据
pheno_data_path <- dir('数据存放地址',
    full.names = TRUE,pattern = 'EOS')
eos <- raster(pheno_data_path[19]) # 读取物候数据

## 使用eos的栅格分辨率,使用n_depo_ave的范围
std_raster <- resample(n_depo_ave,eos[[1]])  # 以eos对氮沉降进行重采样
std_raster <- crop(std_raster,n_depo_ave)
eos_with_ndepo <- crop(eos[[1]],std_raster)
eos_with_ndepo <- mask(eos_with_ndepo,std_raster)  # 用氮沉降数据范围截取eos
plot(eos_with_ndepo)

栅格重分类

  1. 建立重分类矩阵
  2. 使用reclassify() 函数
# 建立重分类矩阵
reclassmatrix <- matrix(c(0,2,1,
                          2,4,2,
                          4,6,3,
                          6,8,4,
                          8,10,5,
                          10,12,6,
                          12,14,7,
                          14,16,8,
                          16,18,9,
                          18,20,10,
                          20,70,11),11,3,byrow = TRUE)
n_depo_ave_reclassify2 <- reclassify(n_depo_ave,rcl = reclassmatrix)
plot(n_depo_ave_reclassify2)

在这里插入图片描述

  • 31
    点赞
  • 168
    收藏
    觉得还不错? 一键收藏
  • 41
    评论
评论 41
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值