将哨兵5P原始数据进行地理校正
#!usr/bin/env python
# -*- coding: utf-8 -*-
"""
date:2021/6/23
author:甲戌_Tr
email: liu_xxxi@163.com
"""
import os
import random
import string
from osgeo import gdal
import numpy as np
class S5Preprocess:
'''
将S5P原始数据转换为tif格式,默认坐标系为84坐标系
'''
def __init__(self):
# 默认处理波段,也可以自行更改或添加
self.field = {'HCHO':['formaldehyde_tropospheric_vertical_column_precision'],
'NO2': ['nitrogendioxide_tropospheric_column_precision'],
'O3': ['ozone_total_vertical_column_precision']}
def write_s5p_tif(self,in_filepath, variables, output_folder, EPSG_code="4326", spatial_res=[]):
"""
将特定变量的S5P nc文件转换为带坐标系的tif文件
S5P默认栅格大小为3500*7000
in_filepath: str nc文件全路径
variables: list of str 所需转换的变量名称
EPSG_code: str 坐标系
spatial_res: list of float, x 和 y 像元大小 默认为3500*7000
"""
flag = False
kwargs = {}
kwargs['EPSG_code'] = EPSG_code
kwargs['spatial_res'] = spatial_res
# Create vrt files for latitude and longitude variables
geo_params = {}
geo_params['outputSRS'] = f"EPSG:{EPSG_code}"
gdal.Translate("lat.vrt", f'HDF5:"{in_filepath}"://PRODUCT/latitude', **geo_params)
lat_ds = gdal.Open("lat.vrt")
gdal.Translate("lon.vrt", f'HDF5:"{in_filepath}"://PRODUCT/longitude', **geo_params)
random_data = ''.join(random.sample(string.ascii_letters + string.digits, 6)) + '.tif'
random_mask = ''.join(random.sample(string.ascii_letters + string.digits, 6)) + '.tif'
data_var = random_data
mask_var = random_mask
for variable in variables:
output_file = self.generate_out_filepath(in_filepath, output_folder, variable, ".tif")
# Georeference variable datset
self.write_var_to_tif(data_var, in_filepath, variable, "lon.vrt", "lat.vrt", lat_ds,
**kwargs)
# Georeference quality_value datset
self.write_var_to_tif(mask_var, in_filepath, "qa_value", "lon.vrt", "lat.vrt", lat_ds,
**kwargs)
# Apply quality mask to variable dataset
if os.path.isfile(data_var) and os.path.isfile(mask_var):
self.write_masked_data(output_file, data_var, mask_var, self.threshold)
flag = True
# Clean
os.remove(data_var)
os.remove(data_var.split('.')[0] + '_.vrt')
os.remove(mask_var)
os.remove(mask_var.split('.')[0] + '_.vrt')
else:
flag = False
# Clean
os.remove(data_var.split('.')[0] + '_.vrt')
os.remove(mask_var.split('.')[0] + '_.vrt')
# Remove vrt files
os.remove('lon.vrt')
os.remove('lat.vrt')
return flag
def write_var_to_tif(self,out_filepath, in_filepath, variable, lon_file, lat_file, ds, **kwargs):
'''
动态撰写vrt虚拟栅格文件
:param out_filepath: 输出文件路径
:param in_filepath: 输入文件路径
:param variable: 变量名称
:param lon_file: 经度vrt文件
:param lat_file: 纬度vrt文件
:param ds: 包含输出信息的dataset
:param kwargs:不定长参数 包含坐标信息 分辨率等
:return: None
'''
vrt_filepath = out_filepath.split('.')[0] + '_.vrt'
with open(vrt_filepath, "w") as text_file:
text_file.write(f"""
<VRTDataset rasterXSize="{ds.RasterXSize}" rasterYSize="{ds.RasterYSize}">
<metadata domain="GEOLOCATION">
<mdi key="X_DATASET">{lon_file}</mdi>
<mdi key="X_BAND">1</mdi>
<mdi key="Y_DATASET">{lat_file}</mdi>
<mdi key="Y_BAND">1</mdi>
<mdi key="PIXEL_OFFSET">0</mdi>
<mdi key="LINE_OFFSET">0</mdi>
<mdi key="PIXEL_STEP">1</mdi>
<mdi key="LINE_STEP">1</mdi>
</metadata>
<VRTRasterBand band="1" datatype="Float32">
<SimpleSource>
<SourceFilename relativeToVRT="0">HDF5:{in_filepath}://PRODUCT/{variable}</SourceFilename>
<SourceBand>1</SourceBand>
<SourceProperties RasterXSize="{ds.RasterXSize}" RasterYSize="{ds.RasterYSize}" DataType="Float32" BlockXSize="{ds.RasterXSize}" BlockYSize="{ds.RasterYSize}" />
<SrcRect xOff="0" yOff="0" xSize="{ds.RasterXSize}" ySize="{ds.RasterYSize}" />
<DstRect xOff="0" yOff="0" xSize="{ds.RasterXSize}" ySize="{ds.RasterYSize}" />
</SimpleSource>
</VRTRasterBand>
</VRTDataset>
""")
# Add georeferencing to vrt file
self.georef_data(out_filepath, vrt_filepath, vrt=False, **kwargs)
def write_masked_data(self,out_filepath, data_file, mask_file, mask_threshold):
'''
掩膜文件输出
:param out_filepath: 输出文件全路径
:param data_file: 输入处理数据
:param mask_file: 输入掩膜数据
:param mask_threshold: 掩膜阈值
:return: None
'''
data_ds = gdal.Open(data_file)
data = data_ds.ReadAsArray()
mask = gdal.Open(mask_file).ReadAsArray()
data[mask <= mask_threshold] = np.nan
driver = gdal.GetDriverByName('GTiff')
dataset = driver.Create(
out_filepath,
data_ds.RasterXSize,
data_ds.RasterYSize,
1,
gdal.GDT_Float32, )
dataset.SetGeoTransform(data_ds.GetGeoTransform())
dataset.SetProjection(data_ds.GetProjectionRef())
dataset.GetRasterBand(1).WriteArray(data)
dataset.FlushCache() # Write to disk.
def georef_data(self,out_filepath, in_filepath, vrt, EPSG_code, spatial_res):
'''
文件写出
:param out_filepath: 输出文件路径
:param in_filepath: 输入文件
:param vrt: True: 输入文件不是vrt,False: 输入文件为vrt
:param EPSG_code:坐标系编码
:param spatial_res:分辨率
:return:None
'''
params = {}
params['dstNodata'] = -9999
params['dstSRS'] = f"EPSG:{EPSG_code}"
if vrt:
params["format"] = "VRT"
out_filepath = out_filepath.replace('.tif', '.vrt')
else:
params["format"] = "Gtiff"
out_filepath = out_filepath.replace('.vrt', '.tif')
if not spatial_res:
if EPSG_code == "4326":
params['xRes'] = 0.06288 # equivalent to 7000 meters
params['yRes'] = 0.06288 # equivalent to 7000 meters
else:
params['xRes'] = 7000 # meters
params['xRes'] = 7000 # meters
else:
params['xRes'] = spatial_res[0]
params['xRes'] = spatial_res[1]
try:
gdal.Warp(out_filepath, in_filepath, **params)
except Exception as e:
print(e)
def generate_out_filepath(self,in_filepath, output_folder, v_name,ext='.tif'):
'''
输出文件命名
:param in_filepath:输入文件名称全路径 str
:param output_folder:输出文件夹 str
:param v_name:变量名称 str
:param ext:输出格式
:return:输出文件名
'''
f_path, f_name = os.path.split(in_filepath)
filename, txt = os.path.splitext(f_name)
var_name_short = filename[13:20].replace('_','')
timestamp = filename[20:35]
return output_folder + os.path.sep + var_name_short + '_' + timestamp + '_' + v_name + ext
def get_product_name(self,in_filepath):
f_path,f_name = os.path.split(in_filepath)
basename,txt = os.path.splitext(f_name)
var_name_short = basename[13:20].replace('_', '')
return var_name_short
def main(self,files,outpath,value=75):
self.threshold = value
crashed_file = {'crashed_image':[]}
x = []
msg_outpath = []
for f in files:
name = S5Preprocess().get_product_name(f)
var = self.field[name]
outpath_ = outpath + os.sep + name
if not os.path.exists(outpath_):
os.mkdir(outpath_)
msg_outpath.append(outpath_)
flag = self.write_s5p_tif(f, var, outpath_)
if not flag:
x.append(f)
crashed_file['crashed_image'] = x
if x ==[]:
msg = None
else:
msg = crashed_file
msg_out = {'outpath': msg_outpath}
return msg,msg_out
if __name__ == '__main__':
files= ['待处理影像1全路径','待处理影像2全路径',.......] # List
outpath = '输出路径' # str
value = 75 # 质量控制值, int
result = S5Preprocess().main(files,outpath,value)
处理代码参考stackoverflow上内容,
但由于原始链接找不见,后续查到后会进行补充。
欢迎指正~