计算两个经纬度直接的函数
这里用到的是m_map包里带的distance函数，m_map的安装参考这个帖子MATLAB--m_map库的安装 - 简书 (jianshu.com)
lat_wh = 30.60; lon_wh = 114.05;
earth_rad = 6371.004; % 地球半径
dist_rad = distance(lat_1,lon_1,lat_2,lon_2); % 计算两组经纬度之间的弧度
distance = dist_rad/180*pi*earth_rad; % 根据弧长公式，弧度*半径=弧长
想根据经纬度范围选取目标城市的卫星过境数据，这里以L2的气溶胶廓线数据及武汉为例
将100km半径范围内的数据平均，有数据的时候输出为txt
%% 处理原始hdf文件，以武汉为中心，100km为搜索半径
%% daytime
clear;clc;
for iyear=2013:2013
iyear
path_dir = 'E:\CALIOP\profile\';
filenames = ls([path_dir, num2str(iyear),'\*ZD_Subset.hdf']); %文件地址
% filenames = ls([path_dir, num2str(iyear),'\*ZN_Subset.hdf']); %文件地址
for ifile = size(filenames,1):size(filenames,1)
file_name = [path_dir, num2str(iyear),'\',filenames(ifile,:)];
info_L2 = hdfinfo(file_name); % 读取hdf文件信息
dsetsL1 = info_L2.SDS; % 读取SDS信息，包括具体变量的名字
% 读取具体变量
Time= hdfread(file_name, 'Profile_Time'); % TAI时间，用秒记录
lat = hdfread(file_name, 'Latitude');%纬度
lon = hdfread(file_name, 'Longitude');%经度
ext_532 = hdfread(file_name, 'Extinction_Coefficient_532');
CAD = hdfread(file_name, 'CAD_Score');
qc_532 = hdfread(file_name, 'Extinction_QC_Flag_532');
un_ext = hdfread(file_name, 'Extinction_Coefficient_Uncertainty_532');
% quality assurance
ext_532(ext_532<0 | ext_532>1.25) = -9999; % 剔除明显异常值
% CAD 在-100和-20之间，un_ext<10, qc_flag = 0,1,2,16,18
for irow = 1:size(ext_532,1)
for icol = 1:size(ext_532,2)
if CAD(irow,icol,1)>=-100 & CAD(irow,icol,1)<=-20 & un_ext(irow,icol)<=10 ...
& (qc_532(irow,icol) == 0 | qc_532(irow,icol) == 1 | qc_532(irow,icol) == 2 | qc_532(irow,icol) == 16 | qc_532(irow,icol) == 18)
ext_532_final(irow,icol) = ext_532(irow,icol);
else
ext_532_final(irow,icol) = -9999;
end
end
end
% 计算所有格点到武汉的距离
earth_rad = 6371.004;
lat_wh = 30.60; lon_wh = 114.05;
dist = distance(lat_wh,lon_wh,lat(:,2),lon(:,2))/180*pi*earth_rad;
% 挑选100km范围内的数据
index_cir = find(dist<=100);
ext_532_final(ext_532_final==-9999)=nan;
if length(index_cir)>1
ext_532_avg = nanmean(ext_532_final(index_cir,:),1);
ext_532_avg(isnan(ext_532_avg))=-9999;
output_filename = [file_name(end-31:end-28) file_name(end-26:end-25) file_name(end-23:end-22) '_WH_ext532.csv'];
csvwrite(strcat('S:\2.研究方向Ongoing\1.两湖地区\CALIOP卫星\out_daily_wh_150km\daytime_',output_filename), ext_532_avg');
% csvwrite(strcat('S:\2.研究方向Ongoing\1.两湖地区\CALIOP卫星\out_daily_wh_150km\nighttime_',output_filename), ext_532_avg');
end
end
end
