Sentinel2-L2A数据的下载、读取及基础操作

哨兵2数据处理流程 专栏收录该内容
1 篇文章 1 订阅

写在前面
本文结束于主成分分析,中间因为疫情、网络等等原因一拖再拖,到今天终于告一段落。人在忙的时候祈求假期,闲下来又盼望着做点事情。因为网络原因获取不到影像,陷于停滞之中,几乎不敢打开电脑,生怕罪恶感汹涌而来吞噬自己。好在天灾不会永远持续,生活总会踟蹰着向前。2020年6月3日,山东省部分高校解禁,暌违数月的实验室里,第一次打开电脑的我,点击链接,下载影像,查找资料,光标继续跳动,我的生活又开始向前。
2020年8月10日,我在这篇文章里敲下最后一行字,某片不知名的云翳也就此消散。

本文姑且算是我在学习过程中的笔记,可能一次写不完,也可能会遇到各种各样的错误,希望大家批评指正,排错的过程我也会记录下来。有解决不了的问题,我也会发出来,因为大家都可能遇到,希望看到的人可以参与到讨论中来。

2020年2月1日开始,连接欧空局服务器一直失败。 换了两处网络,还分别挂了VPN,所以应该可以排除网络问题; 旧账号连接失败,注册的新账号也一直连接不上,也可以排除账号问题; 分别使用opera、chrome、IDM、迅雷下载过,都失败了,大概也不是下载软件问题。大家有知道是怎么回事的吗?是被封了吗?

2020年6月3日,搁置数月后再度尝试连接欧空局成功。单幅影像下载能跑满带宽,单账号最多可同时下载3幅影像,在使用IDM时应调整最大通道数为3或2。使用订阅链接进行多幅下载仍有15%的失败概率,同样是在IDM中调整为2或3通道。

下面进入正题

哨兵2号是高分辨率多光谱成像卫星,携带一枚多光谱成像仪(MSI),用于陆地监测,可提供植被、土壤和水覆盖、内陆水路及海岸区域等图像,还可用于紧急救援服务。分为2A和2B两颗卫星。

第一颗卫星哨兵2号A于2015年6月23日01:52 UTC以“织女星”运载火箭发射升空。
第二颗卫星哨兵2号B于2017年3月07日北京时间9时49分 UTC以“织女星”运载火箭发射升空。

欧洲航天局介绍说,“哨兵-2B”卫星与2015年6月发射的“哨兵-2A”卫星为同一组,携带高分辨率多光谱成像装置,主要用于监测土地环境,可提供有关陆地植被生长、土壤覆盖状况、内河和沿海区域环境等信息,不仅对改善农林业种植、预测粮食产量、保证粮食安全具有重要意义,还可用于监测洪水、火山喷发、山体滑坡等自然灾害,为人道主义救援提供帮助。两者同时进入运行状态后,每5天可完成一次对地球赤道地区的完整成像,而对于纬度较高的欧洲地区,这一周期仅需3天。

哨兵-2号卫星携带一枚多光谱成像仪(MSI),高度为786km,可覆盖13个光谱波段,幅宽达290千米。地面分辨率分别为10m、20m和60m、一颗卫星的重访周期为10天,两颗互补,重访周期为5天。从可见光和近红外到短波红外,具有不同的空间分辨率,在光学数据中,哨兵-2号数据是唯一一个在红边范围含有三个波段的数据,这对监测植被健康信息非常有效。
其波段如下:
在这里插入图片描述从2015年12月3日起,哨兵2A(Sentinel2A)数据正式向全球用户提供免费下载。

产品级别:

Level-0: 原始数据。

Level-1A:包含元信息的几何粗校正产品。

Level-1B:辐射率产品,嵌入经GCP优化的几何模型但未进行相应的几何校正。

Level-1C:经正射校正和亚像元级几何精校正后的大气表观反射率产品。

欧空局(ESA)仅发布了哨兵2号(S2)的L1C级多光谱数据(MSI),Sentinel-2 L1C是经过正射校正和几何精校正的大气表观反射率产品,并没有进行大气校正。同时,ESA还对S2 L2A级数据进行了定义,L2A级数据主要包含经过大气校正的大气底层反射率数据(Bottom-of-Atmosphere corrected reflectance),但这个L2A数据需要用户根据需求自行生产,为此,ESA发布了专门生产L2A级数据的插件Sen2cor。目前,SNAP对Sen2cor的支持并不好,很难在SNAP中直接调用Sen2cor,所以,这里介绍下sen2cor的命令行配置步骤。当然,L2A级数据还包含一些别的产品,如气溶胶厚度(Aerosol Optical Thickness, AOT)、大气水蒸气(Water Vapour Map, WVM)等。

注:目前国内用户仅能下载Level-1C数据。2A级产品要用户自己进行处理生产。

一、数据下载

1、https://scihub.copernicus.eu/(官网下载)

open hub > 进入先注册,善用谷歌翻译。
然后按步骤可以进行下载。(建议通过IDM下载器,从https://earthexplorer.usgs.gov/,科学上网下载。)

2、USGS下载:https://blog.csdn.net/qq_41718357/article/details/100092856
在这里插入图片描述

二、Windows下预处理

关于哨兵2号预处理:

主要是对1C级产品进行辐射定标和大气校正,制作为L2A数据。Sen2cor是ESA发布的专用于哨兵预处理的插件,下面以sen2cor为例。主要参考http://blog.sciencenet.cn/home.php?mod=space&uid=3367669&do=blog&id=1085133一文并附上截图。

Sen2cor官方下载地址http://step.esa.int/main/third-party-plugins-2/sen2cor/在这里插入图片描述在这里插入图片描述
在这里插入图片描述进入数据所在文件夹下,键入命令:
L2A_Process J:\s\L1C_T50QNM_A008566_20170211T025103\S2A_MSIL1C_20170211T024831_N0204_R132_T50QNM_20170211T025103.SAFE
在这里插入图片描述在这里插入图片描述批处理可参考:https://blog.csdn.net/qq_41718357/article/details/89023618

Windows下使用SNAP进行其他操作

哨兵系列数据推荐在SNAP中进行操作。SNAP可以在欧空局进行下载,这里放出一个百度云链接: https://pan.baidu.com/s/1uvjPZtZjxydXPej1TgaMmw 提取码: ak2p
我下载的是哨兵L2A的数据,无需在进行预处理。如果是L1数据需要在命令行调用Sen2cor预处理后再导入SNAP中。不再赘述。
导入SNAP后可进行诸如叠加,镶嵌,裁减以及分类等等操作,并且可以方便地使用搜索框进行功能搜索(前提是你能拿的准对应的英文)。但是对于我来说,相比SNAP,还是ENVI更加符合我的使用习惯。故本次只是将原影像进行格式转换并导出。导出操作需要对数据进行重采样,搜索框键入“S2 resampling processor”,打开相应功能。界面及设置如下:
在这里插入图片描述
在这里插入图片描述
重采样可以直接改变保存格式,故无需再进行导出。展示一下导出操作的位置。
在这里插入图片描述

Windows下使用ENVI进行处理

虽然哨兵2只有区区12波段,在进行过辐射矫正后更是只剩11个,但我还是用不上。故而第一步就是对数据进行降维。直接选择数个波段对于高精度分类来说太过粗糙,最常见的还是对影像进行特征转换与提取。常见的方法包括最常见的PCA和ICA,核函数引入和ISOMAP,以及较新的等。这里我将上述的方法依次使用了一遍,核函数也选取了多种,分别成图。下面以PCA为例。

1.数据集制作

SNAP转换的ENVI格式是分波段存储的,对于后续操作较为不便。因此我们要进行转换。打开ENVI的classic界面,如图,选择File一栏,打开要处理的影像。
在这里插入图片描述

而后点击File中的Save File As操作,子菜单选择“ENVI Standard”
在这里插入图片描述
点击Import File导入文件,Recorder File更改文件排序,选择合适的路径及名称,点击OK。关闭Classic。

2.主成分分析

在ENVI中打开上一步导出的影像,找到PCA运行程序,路径如图
在这里插入图片描述
打开后首先导入原数据
在这里插入图片描述
而后调整参数,主成分数量根据具体使用范畴而定,我的选择是5个。
在这里插入图片描述
成图。图示为第三主成分
在这里插入图片描述

图形化界面可以灵活地对数据进行需要的操作,但是对于排错等并不占优。下面我们在Ubuntu环境中使用代码对Sentinel2-L2A数据进行操作。

Ubuntu下处理哨兵2数据

一、GDAL读取

打开终端,切换到配置了GDAL的环境,激活打开jupyter notebook。新建空代码。
在这里插入图片描述
哨兵2数据的结构和HDF格式类似,即一个压缩包,包含jpeg2000格式的栅格影像,一些mosaic和几个元数据文本。

首先,调用相关的模组。

from osgeo import gdal
from osgeo.gdalconst import *

注册所有格式,或单个注册

gdal.AllRegister() #注册所有

driver = gdal.GetDriverByName('HFA')
driver.Register()#注册ERDAS数据格式

打开数据

dt = gdal.Open('*********.***')
print(type(dt))#某些如jupyter中可直接省略print
#读取失败的dt会是Nonetype格式,如果dt的格式为osgeo.database或类似,则证明读取成功。

下面我们可以查看数据的相关信息,比如使用dt.RasterCount来获取波段数目,该命令不是方法而是属性,因此无需加括号。我们以我们要分析的哨兵2为例,使用该命令返回的波段数目是0,这证明该数据以数据集形式存储且以子数据集形式组织,需要进一步读取其子数据集。

sd = dt.GetSubDatasets()

返回如下值:
[(‘SENTINEL2_L2A:/vsizip//home/hao/文档/S2B_MSIL2A_20200720T024549_N0214_R132_T51STA_20200720T062139.zip/S2B_MSIL2A_20200720T024549_N0214_R132_T51STA_20200720T062139.SAFE/MTD_MSIL2A.xml:10m:EPSG_32651’, ‘Bands B2, B3, B4, B8 with 10m resolution, UTM 51N’), (‘SENTINEL2_L2A:/vsizip//home/hao/文档/S2B_MSIL2A_20200720T024549_N0214_R132_T51STA_20200720T062139.zip/S2B_MSIL2A_20200720T024549_N0214_R132_T51STA_20200720T062139.SAFE/MTD_MSIL2A.xml:20m:EPSG_32651’, ‘Bands B5, B6, B7, B8A, B11, B12 with 20m resolution, UTM 51N’), (‘SENTINEL2_L2A:/vsizip//home/hao/文档/S2B_MSIL2A_20200720T024549_N0214_R132_T51STA_20200720T062139.zip/S2B_MSIL2A_20200720T024549_N0214_R132_T51STA_20200720T062139.SAFE/MTD_MSIL2A.xml:60m:EPSG_32651’, ‘Bands B1, B9 with 60m resolution, UTM 51N’), (‘SENTINEL2_L2A:/vsizip//home/hao/文档/S2B_MSIL2A_20200720T024549_N0214_R132_T51STA_20200720T062139.zip/S2B_MSIL2A_20200720T024549_N0214_R132_T51STA_20200720T062139.SAFE/MTD_MSIL2A.xml:TCI:EPSG_32651’, ‘True color image, UTM 51N’)]

这说明哨兵2a数据有4个子数据集,上述说明了子数据集归属、路径、分辨率、波段、名称和投影等。使用type查看发现该返回值是一个4x3的list,每一行第一个值是子数据集路径,第二个值是子数据集注释,第三个是子数据投影。进一步读取子数据集数据:

R10m = gdal.Open(sd[0][0])#按列表读取,00为第一子数据集的路径
print(type(R10m))#输出值是<class 'osgeo.gdal.Dataset'>,证明成功
b10m = R10m.ReadAsArray()#将波段以数组形式进行读取
#b10m = R10m.ReadRaster()#将波段以二进制形式进行读取

数组形式比较直观,在进行处理操作时我更倾向于使用第一个函数。在数据压缩等应用领域,第二个函数可能更加适合。

二、主成分分析

该代码获取自网络,时间久远,早已丢失来源,在这里鸣谢这位作者。我对源代码进行了一些魔改,可能并不适合阅读,故对原代码添加注释,有需要的可自行拿取更改。需要注意的是,如同我在前一章节所说,由于哨兵波段分辨率不同,在进行PCA之前必须进行重采样。GDAL对于超采样是有支持的,但效果不好,且效率较低,故建议在SNAP中进行超采样之后,再进行PCA,这里我与原作者都默认输入影像全波段分辨率一致。

from osgeo import gdal
from osgeo.gdalconst import *
import numpy as np

gdal.AllRegister()

#读取单幅影像
rt = '/home/hao/文档/Win_P/20200720a'
ris = gdal.Open(rt)

sd = ris.GetSubDatasets()

#PCA需要将所有波段导入
def meanX(dataX):
    
    return np.mean(dataX, axis=0) #axis=0表示按照列来求均值,如果输入list,则axis=1

def variance(X):
    m, n = np.shape(X)
    mu = meanX(X)
    muAll = np.tile(mu, (m, 1))
    X1 = X - muAll
    variance = 1. / m * np.diag(X1.T * X1)
    return variance

def normalize(X):
    m, n = np.shape(X)
    mu = meanX(X)
    muAll = np.tile(mu, (m, 1))
    X1 = X - muAll
    X2 = np.tile(np.diag(X.T * X), (m, 1))
    XNorm = X1 / X2
    return XNorm
#定义PCA函数,实现主成分分析
def pca(XMat, k):
    average = meanX(XMat) #计算均值,为进行归一化做准备
    m, n = np.shape(XMat) #行,列(7)
    data_adjust = []
    avgs = np.tile(average, (m, 1))
    covX = np.cov(data_adjust.T) # 计算协方差矩阵
    featValue, featVec = np.linalg.eig(covX) # 求解协方差矩阵的特征值和特征向量
    index = np.argsort(-featValue) # 按照featValue进行从大到小排序
#确定主成分的个数

sum=0
max=0

for i in range(len(index)):
    sum=sum+featValue[index[i]]

if max/sum<0.99:
    max=max+featValue[index[i]]
    k=k+1
    finalData = []
    print ("k = ",k)
    
if k > n:
    print ('k must lower than feature number')
    return
else:
    # 注意特征向量是列向量,而numpy的二维矩阵(数组)a[m][n]中,a[1]表示第1行值
    selectVec = np.matrix(featVec[index[:k]]) # 找出k个特征值对应的特征向量
    reconData = (data_adjust * selectVec.T) + average[:,6]
    #将m*n的数据集乘以k个n维的特征向量的特征向量(m * k),得到最后降维的数据。
    print( np.shape(reconData))

array = []
array = np.array(array) # 列表转数组

for j in range(k):
    array = np.append(array, reconData[:,j]) #将reconData从二维(3列)数组转化从一维数组(1维)
    print (np.shape(array))
    array1 = array.reshape(k, im_height, im_width) #将 一维数组 转成 3维矩阵(k,高(行),宽(列))
    out = ga.SaveArray(array1, os.path.join(path, 'after.img'), format='GTiff', prototype=img)
    return reconData
#打开img或者tif格式的图像进行PCA
path = “… 主成分分析”
input = os.listdir(path)
data = []

for i in input:
    if i == ‘before.img’:
        img = gdal.Open(os.path.join(path, i))
        im_width = img.RasterXSize # 栅格矩阵的列数
        im_height = img.RasterYSize # 栅格矩阵的行数
        im_bands = img.RasterCount # 波段数
        im_data = img.ReadAsArray(0, 0, im_width, im_height) # 获取数据 将数据写成数组,对应栅格矩阵
        #im_data:(波段数,行,列)

for j in range(im_bands):
    picture = im_data[j,:,:].flatten() # 变成一维数组
    data.append(picture) #将各个波段合并在一个数组里,一行为一个波段(n个波段n行)
    print (“data维数”, np.shape(data))
    data = np.mat(data) #数组转换为矩阵
    data1 = data.T #进行转置
    print (“data1维数”, np.shape(data1))
    XMat = data1
  • 5
    点赞
  • 8
    评论
  • 39
    收藏
  • 打赏
    打赏
  • 扫一扫,分享海报

评论 8 您还未登录,请先 登录 后发表或查看评论
©️2022 CSDN 皮肤主题:1024 设计师:白松林 返回首页

打赏作者

大耳朵的猩猩

你的鼓励将是我创作的最大动力

¥2 ¥4 ¥6 ¥10 ¥20
输入1-500的整数
余额支付 (余额:-- )
扫码支付
扫码支付:¥2
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

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

余额充值