- 利用修改后的griddata1.m来插值，避免了拉丝
%function：将S5P的值读取出来，处理成网格存入csv
%author：HeQin
%ncread的原型：ncread(source,varname,start,count,stride)
clc; %清屏
clear; %清空
datadir='D:\Shared Data\newEvaluation\S5P_NO2\RPRO\'; %指定批量数据所在的文件夹
filelist=dir([datadir,'*.nc']); %指定批量数据的类型
for s=1:length(filelist)
filename=[datadir,filelist(s).name];
ncid=netcdf.open(filename,'NC_NOWRITE');
ncdisp(filename); %在命令窗中显示nc文件的变量
LON = ncread(filename,'/PRODUCT/longitude'); %经度
LAT = ncread(filename,'/PRODUCT/latitude'); %纬度
NO = ncread(filename,'/PRODUCT/nitrogendioxide_tropospheric_column');
NO(isnan(NO))= nan;
NO(NO<0) = nan;
%size_1=size(LAT);
%lat = reshape(LAT,1,size_1(1)*size_1(2));
%lon = reshape(LON,1,size_1(1)*size_1(2));
%no = reshape(NO,1,size_1(1)*size_1(2),1);
%locate = find(no<0); %a是存储数据的数组名，find是找到 的数的位置
%lat(locate) = []; %删除数组a中 的元素
%lon(locate) = []; %删除数组a中 的元素
%no(locate) = []; %删除数组a中 的元素
latlat=3:0.125:54;
lonlon=73:0.125:136;
[X,Y] = meshgrid(lonlon,latlat);
vq= griddata1(LON,LAT,NO,X,Y,'linear');
vq(isnan(vq))= -1; %标记出异常值
vq = vq*1000000;
csvwrite([filename,'.csv'],vq);
netcdf.close(ncid); % 关闭文件
end
- 后来改为全球时出现边缘拉丝，就用python写了个最邻近的插值
# -*- coding:utf-8 -*-
import os
import glob
import h5py
import numpy as np
convert_to_E17moleccm2 = 602.214
package = glob.glob('F:\\S5P\\6.28之后\\*.nc')
print("file found {}".format(package))
for i in range(len(package)):
data = package[i]
print("-----------",i,data)
with h5py.File(data,'r') as f:
co = np.zeros((1440,2880))
count = np.zeros((1440,2880),dtype=int)
#layer = f['PRODUCT']['layer'][:]
lat = f['PRODUCT']['latitude'][:,:][0] #[0]是固定加上的，可能它涉及到了不同的layer,但一般只会有这一层的
lon = f['PRODUCT']['longitude'][:,:][0]
no2_trop_origin = f['PRODUCT']['carbonmonoxide_total_column'][:,:][0]
no2_trop_origin = np.array(no2_trop_origin, dtype=np.float64) #否则会溢出
no2_trop = np.where(no2_trop_origin<1,no2_trop_origin*convert_to_E17moleccm2,no2_trop_origin*(-1)) #正常值乘上系数做单位转换
no2_trop = np.where(no2_trop>0,no2_trop,0)
for y in range(lon.shape[0]):
for x in range(lon.shape[1]):
lonG = int((lon[y][x]+180)*8)
latG = 1440 - int((lat[y][x]+90)*8)
if lonG == 2880:
lonG = 2779
if latG == 1440:
latG = 1339
temp = co[latG][lonG] * count[latG][lonG]+no2_trop[y][x]
count[latG][lonG] = count[latG][lonG]+1
co[latG][lonG] = temp /count[latG][lonG]
np.save(data,co)
其他人写的一些也没有很好的解决
https://github.com/braunwiediefarbe/He5ToGeotiff之后处理成月均值的程序，偷了个懒，没有对一天内overlap的平均后再取均值，而是直接取了均值
# -*- coding:utf-8 -*-
import os
import glob
import h5py
import numpy as np
package = glob.glob('H:\\S5P\\npy\\3\\*.npy')
print("file found {}".format(package))
avg = np.zeros((1440,2880))
count = np.zeros((1440,2880),dtype=int)
for i in range(len(package)):
data = package[i]
print("-----------",i,data)
one = np.load(data)
#print(one[1210][2612])
temp = avg * count+ one
#print(temp[1210][2612])
count = np.where(one>0,count+1,count)
#print(count[1210][2612])
avg = np.where(count>0,temp /count,0) #注意避免count为0的情况
#print(avg[1210][2612])
np.savetxt('3.txt',avg)
补充：
- 简单的显示出矩阵对应的图案
import matplotlib.pyplot as plt
plt.imshow(avg)
- 保存为txt时设置自己想要的格式
np.savetxt('77.txt',avg,fmt='%.6f')
