今天介紹一下Sentinel衛星以及一些預處理的方法。
1.基本信息（成像儀/重訪周期/波段數/分辨率）
哨兵2號是高分辨率多光譜成像衛星，攜帶一枚多光譜成像儀(MSI)，用於陸地監測，可提供植被、土壤和水覆蓋、內陸水路及海岸區域等圖像，分為2A和2B兩顆衛星,哨兵，2B與2015年6月發射的哨兵-2A衛星為同一組。哨兵-2號衛星高度為786km，覆蓋13個光譜波段，幅寬達290千米。地面分辨率分別為10m、20m和60m、一顆衛星的重訪周期為10天，兩顆互補，重訪周期為5天。從可見光和近紅外到短波紅外，具有不同的空間分辨率，在光學數據中，哨兵-2號數據是唯一一個在紅邊范圍含有三個波段的數據，這對監測植被健康信息非常有效。
下載網站：https://scihub.copernicus.eu/dhus/#/home
Ps:夜里1:00下載會很快
成像儀：MSI 重訪周期：5-10天 波段數：13 分辨率：10m/20m/60m
2.產品等級及插件介紹
歐空局僅發布了哨兵2號的L1C級多光譜數據（MSI），Sentinel-2 L1C是經過正射校正和幾何精校正的大氣表觀反射率產品，並沒有進行大氣校正。同時，ESA還對S2 L2A級數據進行了定義，L2A級數據主要包含經過大氣校正的大氣底層反射率數據（Bottom-of-Atmosphere corrected reflectance），這個數據可以通過Sen2cor插件自行生產。
插件說明：1.Sen2cor 分為2.05/2.08版本，前者用於處理290km寬幅的16年老數據，后者用於處理新數據。L2A Process + 文件名完成輻射定標+大氣校正。
2.Sen2Res，提供超分辨率合成功能。這是一個nbm文件，由於Sentinel2有10/20/60m三個分辨率的遙感數據，所以在對圖像進行其他處理之前，需要先統一到一個分辨率，可以用重采樣或者是超分辨率合成。具體操作步驟是在Snap中的plugins里添加已下載的插件並安裝，此外該nbm文件還可以通過python調用實現批處理功能。
3.文件夾介紹
介紹一下數據的內容，這里分為16年和19年的數據，由於16年的影像是290km幅寬，所以存儲方式和19年有所不同。
首先是16年的數據。其中xml文件是索引文件，索引到存儲的數據，具體存放波段數據的文件是GRANULE文件夾，用10個文件夾把大圖幅切片分別存儲，打開某個圖幅的IMG_DATA文件夾存儲了13個波段的影像數據，jp2格式。通過最外層索引文件控制內層的文件，直接打開所有數據。也可以通過內層文件夾的xml文件打開單個圖幅的數據。AUX_DATA為輔助文件，可能存儲 太陽高度角、日地距離等信息。
19年的數據存儲大小發生改變，減小了單幅圖像的幅寬。GRANULE文件夾中只存放了一個圖幅，存儲所有波段數據（圖4）。
接着對上述L1C級產品加工，做大氣校正，得到L2A級產品，數據格式有所不同。校正后IMG_DATA把數據按照分辨率不同分成了三個文件夾進行存儲。Sentinel2A數據沒有全色波段，因此他提供了10m分辨率的四個波段數據，用於完成圖像的超分辨率。此外10波段作為卷雲波段，不進行大氣校正，沒有提供，需要從L1C產品獲取取。
完成圖像大氣校正后，由於不同波段的分辨率各不相同，所以需要重采樣到同一個分辨率，或者用Sen2Res完成超分辨率合成，速度較慢，不過質量很好。這樣才能進行下一步處理所有波段分辨率統一了才能進行下一步處理，可以Snap保存為ENVI格式或者GTiff等等再做運算，這樣計算NDVI等指數也就方便了。
L2A產品的文件里還有相關的分類文件，分類了植物、雲、雪、水體等等，還有對應的index，但是精度太低。
4.數據的讀取。
1）Snap軟件：打開xml索引文件從而打開所有波段數據，或者打開單個文件xml來打開單個波段。
2）Python+GDAL
JP2格式的圖像數據可以用GDAL讀取，並且可以讀取到文件的投影信息和地理參數，有需要可以用GDAL讀取並且轉換成Tiff文件。
這里還發現數據集ds可以直接調用ds.ReadAsArray()讀取數組陣列，真是有點奇葩。。。我還一直以為ds只能ds.ReadRaster呢。
from osgeo import gdal import numpy as np filename = r'S2A.jp2' class sentinel: def read_img(self, filename): ds = gdal.Open(filename) Xsize = ds.RasterXSize Ysize = ds.RasterYSize bnum = ds.RasterCount data = ds.ReadAsArray(0, 0, Xsize, Ysize) #太奇葩了，ds都能ReadAsArray() proj = ds.GetProjection() geotrans = ds.GetGeoTransform() return proj, geotrans, data, Xsize, Ysize, bnum def write_img(self, filename, outfile): proj, geotrans, data, width, height, bands = self.read_img(filename) driver = gdal.GetDriverByName('GTiff') ds = driver.Create(outfile, width, height, bands, gdal.GDT_UInt16) ds.SetGeoTransform(geotrans) ds.SetProjection(proj) if bands==1: ds.GetRasterBand(1).WriteArray(data) else: for i in range(1, bands+1): ds.GetRasterBand(i).WriteArray(data) del ds run = sentinel() proj, geotrans, data = run.read_img(filename) print(np.shape(data))
另外貼一下16年數據的讀取代碼，有所不同。
from osgeo import gdal import os #這里打開了一個哨兵2A L1C zip文件 filename = (' S2A_MSIL1C_***.zip') root_ds = gdal.Open(filename) # 返回結果是一個list，list中的每個元素是一個tuple，每個tuple中包含了對數據集的路徑，元數據等的描述信息 #讀取所有子數據集 ds_list = root_ds.GetSubDatasets() visual_ds = gdal.Open(ds_list[0][0]) # 取出第一個圖幅的第一個波段數據 visual_arr = visual_ds.ReadAsArray() # 將數據集中的數據轉為ndarray del visual_arr # 獲得柵格數據的一些重要信息 print(f'打開數據為：{ds_list[0][1]}') print(f'投影信息：{visual_ds.GetProjection()}') print(f'柵格波段數：{visual_ds.RasterCount}') print(f'柵格列數（寬度）：{visual_ds.RasterXSize}') print(f'柵格行數（高度）：{visual_ds.RasterYSize}')
以上就是本次學習Sentinel2A的一些心得，其實仔細看看發現和Landsat8都是有共通之處的。下次會介紹一下Sentinel3A衛星，這個數據的存儲方式好像都是.nc文件，有點意思。