前言
本机环境:R4.0 RStudio 4.0
师兄用S5Processor包提取哨兵5p某些波段的数据,但是下载的数据里面有些打不开,程序就中断了。在看了R语法后加了一个异常捕获,遇到异常暂时不处理,直接跳过。程序利用shape裁剪得到图像,方便后期上色,同时保存了一些默认配色方案的png,方便查看。
用途和方法
利用R中的S5Processor,提取sentinel5p中一些波段的shape部分区域的数据。
方法如下:
不制作shape会提取波段中的全部数据,跳过第一步即可。
- 准备shape文件,可以在envi中制做。在envi中打开数据,选择roi,勾选区域,最后点击roi窗口中的file-export-export to shapefile可得到需要的shape文件。
- 下载哨兵5p数据,参考教程:Sentinel-5P数据介绍与预处理。
- 安装S5Processor包。参考教程。先安装devtools,在R Studio的终端中键入:
install.package("devtools")
,在键入devtools::install_github("MBalthasar/S5Processor")
。网络问题可能会安装失败,可以下载到本地安装。 - 在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() 条形图