計(jì)算兩個(gè)經(jīng)緯度直接的函數(shù)
這里用到的是m_map包里帶的distance函數(shù)，m_map的安裝參考這個(gè)帖子MATLAB--m_map庫(kù)的安裝 - 簡(jiǎn)書(shū) (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);? % 計(jì)算兩組經(jīng)緯度之間的弧度
distance = dist_rad/180*pi*earth_rad;? ? ? ? ? ? ?% 根據(jù)弧長(zhǎng)公式誓军，弧度*半徑=弧長(zhǎng)
想根據(jù)經(jīng)緯度范圍選取目標(biāo)城市的衛(wèi)星過(guò)境數(shù)據(jù)拗馒，這里以L2的氣溶膠廓線數(shù)據(jù)及武漢為例
將100km半徑范圍內(nèi)的數(shù)據(jù)平均愉老，有數(shù)據(jù)的時(shí)候輸出為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時(shí)間踪少，用秒記錄
? ? ? ? lat = hdfread(file_name, 'Latitude');%緯度
? ? ? ? lon = hdfread(file_name, 'Longitude');%經(jīng)度
? ? ? ? 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
? ? ? ? % 計(jì)算所有格點(diǎn)到武漢的距離
? ? ? ? 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范圍內(nèi)的數(shù)據(jù)
? ? ? ? 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.兩湖地區(qū)\CALIOP衛(wèi)星\out_daily_wh_150km\daytime_',output_filename), ext_532_avg');
%? ? ? ? ? ? csvwrite(strcat('S:\2.研究方向Ongoing\1.兩湖地區(qū)\CALIOP衛(wèi)星\out_daily_wh_150km\nighttime_',output_filename), ext_532_avg');
? ? ? ? end
? ? end
end
