Sentinel5p提取特定波段数据

目录

前言

本机环境：R4.0 RStudio 4.0
师兄用S5Processor包提取哨兵5p某些波段的数据，但是下载的数据里面有些打不开，程序就中断了。在看了R语法后加了一个异常捕获，遇到异常暂时不处理，直接跳过。程序利用shape裁剪得到图像，方便后期上色，同时保存了一些默认配色方案的png，方便查看。

用途和方法

利用R中的S5Processor，提取sentinel5p中一些波段的shape部分区域的数据。
方法如下：
不制作shape会提取波段中的全部数据，跳过第一步即可。

1. 准备shape文件，可以在envi中制做。在envi中打开数据，选择roi，勾选区域，最后点击roi窗口中的file-export-export to shapefile可得到需要的shape文件。
2. 下载哨兵5p数据，参考教程：Sentinel-5P数据介绍与预处理。
3. 安装S5Processor包。参考教程。先安装devtools，在R Studio的终端中键入：install.package("devtools")，在键入devtools::install_github("MBalthasar/S5Processor")。网络问题可能会安装失败，可以下载到本地安装。
4. 在R studio中运行即可。新建一个script，键入程序，点击Source。

完整程序

#  导入需要的各种包，运行时没有的包按照提示安装即可
library(S5Processor)
library(ncdf4)
library(ggplot2)
library(dismo)
library(maptools)
library(raster)
library(geosphere)
library(rgdal)
library(rgeos)
library(sp)
library(conflicted)

#  部分路径设置。
#  filepath：文件保存的路径，也是工作路径
#  patho：sentinel5p数据所在的顶层文件夹
#  sites：多个地点的子文件夹名字
#  shapename：每个地点对应的shape文件

filepath <- 'F:\\'
patho <- "F:\\"
imagepath <- "F:\\Img\\"
sites <- c("sites1", "sites2", "sites3")
shapename <- c("sites1.shp", "sites2.shp", "sites3.shp")
sitesindex <- length(sites)

for (j in 1:sitesindex) {
  #  setwd：设置工作路径，paste0()方法用于字符串拼接，不带拼接符号。
  #  path：为哨兵数据的目录位置，顶层和地点文件夹拼接得到
  setwd(paste0(filepath, sites[j]))
  path <- paste0(patho, sites[j],"/Sentinel5p")
  shapepath <- paste0(patho,"/Shape/",sites[j],"/",shapename[j])
  vnm_shp <- raster::shapefile(shapepath)
  
  #  file_list：得到所有后缀为nc的文件。
  file_list <- list.files(path, pattern = ".nc")
  #  tmp <- strsplit(file_list, split = ".", fixed = TRUE)  #  暂时未用
  #  file_name <- unlist(lapply(tmp, head, 1))
  
  index <- length(file_list)
  for (i in 1:index) {
    print(file_list[i])
    #  生产保存的png文件绝对路径。substring()字符串切片。
    imagename <- paste0(
      imagepath, sites[j], "/",
      substring(file_list[i], 1, nchar(file_list[i]) - 58), ".png"
    )
    #  将生成的图像保存为png文件
    png(filename = imagename)
    #  R中的异常捕获:try()。silnet=FALSE在终端输出错误详情。
    x.inv <- try(
      #  调用S5Processor()提取vnm_shp区域中的数据，不带my_aoi参数提取全部数据。
      band2 <- S5P_process(
        input = paste(path, file_list[i], sep = "/"),
        product = 6, my_aoi = vnm_shp
      ),
      silent = FALSE
    )
    #  捕获错误，关闭图像，删除生成的空文件，跳过本次循环。
    if ("try-error" %in% class(x.inv)) {
      dev.off()
      file.remove(imagename)
      next
    }
    dev.off()
    tifname <- substring(file_list[i], 1, nchar(file_list[i]) - 58)
    #  writeRaster()用于保存提取到的波段数据，格式为tif。
    writeRaster(band2, tifname, format='GTiff', overwrite = TRUE)
  }
  file_list <- NULL
  tmp <- NULL
  file_name <- NULL
}

附件

部分R语法，参考教程

# 配置参考教程：
#   https://blog.csdn.net/weixin_39287540/article/details/106459170
#   https://sspai.com/post/47386
# 菜鸟教程：R: https://www.runoob.com/r

library(xlsx)
library(XML)
library(rjson)

myString.0 <- "Hello, R"
"Hello,R!" -> myString.1
myString.2 = "hello, R!"
print(myString.0)

# cat(s,s1...)拼接结果
print(cat(myString.0, myString.1, myString.2, sep = "/"))
# cat(context, file="toStorageFile", append=TRUE) 拼接

# sink(toStorageFile) 控制台输出的到文本，不显示
# sink(tofile, split=true) 控制台显示。
# sink() 重置

# readLines(file) 文件读入
lines <- readLines("./readtest.txt")
print(lines[1])

# 工作路径:getwd() setwd()

# 用if做注释
# if(FALSE) {
    # "
    # 多行注释
    # 多行注释
    # "
# }

# 基础运算 +-*/ 逻辑运算&|~！ 整除%/% %% %*% %in% : 求逆矩阵solve()
# 向量、矩阵、列表、数组

# 条件 if...else if...else
# 循环 repeat{break} while() {} for(i in value) {}
# 函数 functinoName <- function(args) {return()}

# 字符串：paste(...,sep='/',collapse=NULL) 指定符号拼接
#        paste(...,collapse="=") 指定符号对应拼接 paste(letterss[1:6], 1:6, collapse="=")
# 格式化：format(x, digits, nsmall, scientific, width, justify = c("left", "right", "centre", "none")) 
# tolower(), toupper(), substring(str, first, last)

# 列表：list
# list(args) names(list) <- c(args) 生产列表，为列表元素命名
# list[index], list$listName <- NULL/NewVaulse 删除/添加元素
# mergedlist <- c(list1, list2, ...) list合并
# unlist() list-->vector

# csv文件：read.csv(file,encoding="UTF-8")
# subset(csvdata, expression'index==1') 查找，返回一行的值
# write.csv(csvvalue, csvfile) 保存

# read.xlsx(xlsxfile,sheetIndex = 1) excel
# write.xlsx() write

# xmlParse(file = "test.xml") read xml file.
# xmlRoot(), xmlSize(), xmlToList(), xmlToDataFrame()

# fromJSON(file = "jsontest.json") json
# as.data.frame() To data frame


# pie(x,labels = label, col = cols) 默认生成pdf，可以通过png(),jpeg(),bmp()生成图像。
# png(filename = "pie.png", width = 100, height = 100) png()
# legend(site,names,cex,fill) 添加颜色样本标注

# barplot() 条形图

本文地址：https://blog.csdn.net/Mr_kuilei/article/details/108252533