SEBAL_manual.rar_matlab水文_sebal_蒸散发_蒸散发 matlab

%function [] =SEBAL(rootpath,year,day,latitudeorfile,demfile,Land_usefile) %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~% % Surface Energy Balance Algorithm for Land (SEBAL) % % (Matlab code) % % by % % Xuejuan Chen % % @ SNNU % % 2013.04.07 % %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~% clc clear all; %------------------------------------------------------------------- %input parameter DOY=161; %---------------- Esun=1367; %solar constant (W/m2) PI=3.1415926; [X,R] = geotiffread('F:\遥感蒸散发\MODISnew\DEM\DEM.tif'); %read image with spatially referenced X为影像矩阵，相当于imread,R为地理信息。 info=geotiffinfo('F:\遥感蒸散发\MODISnew\DEM\DEM.tif'); %%记录地理信息 %----------------- %%%%%%%%%%%%%%%%%计算卫星过境瞬时温度%%%%%%%%%%%%%%%%%%% latitudeor=imread('F:\遥感蒸散发\纬度\latitude-1.tif'); gaoqiwen=imread('气象数据new\MaxTemTif\2002\2002_06_26__MaxTem1.tif'); diqiwen=imread('气象数据new\MinTemTif\2002\2002_06_26__MinTem1.tif'); RH=imread('气象数据new\AvgRhuTif\2002\177\data.tiff'); %平均相对湿度 AP=imread('气象数据new\AvgPrsTif\2002\177\data.tiff'); U2=imread('气象数据new\AvgWinTif\2002\177\data.tiff'); idx=find(U2<0.5); U2(idx)=0.5; % idx=find(U2>4.5); % U2(idx)=4.5; P=1.8; %日最高气温出现时刻与太阳时正午时刻的时间差 （h); time difference of solar noon to highest air temperature in a day PZ=2; %P取整数 TC=4; %夜温变化时间常数，一般取4小时； constant of night temperature variation Tk=15; %计算瞬时温度时的参数 [M,N]=size(gaoqiwen); T_hour1=zeros(M,N); T_hour2=zeros(M,N); Ta=zeros(M,N); latitude=zeros(M,N); omigaS=zeros(M,N); DL=zeros(M,N); NL=zeros(M,N); TimeS=zeros(M,N); TimeR=zeros(M,N); es=zeros(M,N); de=zeros(M,N); delta=zeros(M,N); lamda=zeros(M,N); PsyCon=zeros(M,N); disp(66); for m=1:M for n=1:N solardec=0.409*sin(2*PI*DOY/365-1.39); %太阳磁偏角(赤纬）(solar declination)(rad); latitude(m,n)=PI/180*latitudeor(m,n); omigaS(m,n)=acos(tan(solardec).*-tan(latitude(m,n))); %日落时角度;(rad); sunset angle DL(m,n)=24/PI*omigaS(m,n); %最大日照时数（hour) the maximum sun light hour NL(m,n)=24-DL(m,n); %夜长 night length TimeS(m,n)=12+DL(m,n)/2; %日落时刻； time of sunset TimeR(m,n)=12-DL(m,n)/2; %日出时刻； time of sun rise T_hour1(m,n)=diqiwen(m,n)+(gaoqiwen(m,n)-diqiwen(m,n))*sin(PI*(10-TimeR(m,n))/(DL(m,n)+2*P)); T_hour2(m,n)=diqiwen(m,n)+(gaoqiwen(m,n)-diqiwen(m,n))*sin(PI*(11-TimeR(m,n))/(DL(m,n)+2*P)); Ta(m,n)=0.5*(T_hour1(m,n)+T_hour2(m,n)); %卫星过境时的温度 es(m,n)=0.6108*exp(17.27*Ta(m,n)/(Ta(m,n)+237.3)); %kpa 此处所用温度单位为摄氏度 饱和水汽压 de(m,n)=es(m,n)-es(m,n)*RH(m,n); %水汽压差 delta(m,n)=4098*es(m,n)/(Ta(m,n)+237.3)^2; %Kpa/k 饱和水汽压斜率 lamda(m,n)=2.5e6-2.32e3*Ta(m,n); %J/Kg 蒸发系数 PsyCon(m,n)=1013*AP(m,n)/0.622/lamda(m,n); %Kpa/k 干湿表常数 end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %------------------------------------------------------------------- %Read MODIS images Land=imread('MODISnew\land_use\Land_use.tif'); %land use map Land=single(Land); DEM_0=imread('MODISnew\DEM\DEM.tif'); %DEM map DEM=single(DEM_0); NDVI_0=imread('MODISnew\NDVI-LST\NDVI\MOSAIC_TMP_2001161.hdfout.1_km_16_days_NDV.tif'); NDVI_0=sinqgle(NDVI_0); NDVI_0=NDVI_0.*0.0001; idx=find(NDVI_0==-0.3); NDVI_0(idx)=0; SR=(1+NDVI_0)./(1-NDVI_0); %%%新增新增新增 albedo_0=0.28-0.14.*exp(-6.08./SR.^2); albedo_0=single(albedo_0); geotiffwrite('MODISnew\Albedo-modify\2001161Albedo.tif',albedo_0,R, 'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag); LST_0=imread('MODISnew\NDVI-LST\LST_modify\MOSAIC_TMP_2001161.hdfout.LST_Day_1k_modify.tif'); %land surface temperature %%%新增新增新增 LST_0=single(LST_0); LST_0=LST_0.*0.02; Surface_emissivity_1=imread('MODISnew\MOSAIC_TMP_2001161.hdfout.Emis_31.tif'); %surface emissivity (MOD11A1) Surface_emissivity_1=double(Surface_emissivity_1); Surface_emissivity_2=imread('MODISnew\MOSAIC_TMP_2001161.hdfout.Emis_3.tif'); %surface emissivity (MOD11A1) Surface_emissivity_2=double(Surface_emissivity_2); Sur_Emiss_0=(Surface_emissivity_1.*0.002+Surface_emissivity_2.*0.002)./2+0.49; Sur_Emiss_0=single(Sur_Emiss_0); LAI_0=imread('MODISnew\MOSAIC_TMP_2001161.hdfout.Lai_1k.tif'); LAI_0=single(LAI_0); LAI_0=LAI_0.*0.1; Solar_Zenith_0=imread('MODISnew\MOSAIC_TMP_2001161.hdfout.1_km_16_days_sun_zenith_angl.tif'); %%%太阳天顶角 Solar_Zenith_0=single(Solar_Zenith_0); Solar_Zenith_0=Solar_Zenith_0.*0.01; %----------------------------------------------------------------------- idx=find(LAI_0>=24.9&Land>=11&Land<=19); LAI_0(idx)=0.7271*exp(3.0236*NDVI_0(idx)); idx=find(LAI_0>=24.9&Land==21); LAI_0(idx)=0.5628*SR(idx)+0.3817; idx=find(LAI_0>=24.9&Land>=22&Land<=24); LAI_0(idx)=1.1273*SR(idx)-0.3468; idx=find(LAI_0>=24.9&Land>=31&Land<=33); LAI_0(idx)=0.8253*exp(0.3309*SR(idx)); idx=find(LAI_0>=24.9&Land>33); LAI_0(idx)=0; %----------------------------------------------------------------- [m,n]=size(LAI_0); %------------------------------------------------------------------------------------------------------------------------------- %calculate net radiation and ground_heat_flux 计算净辐射通量和土壤热通量 %------------------------------------------------------------------------------------------------------------------------------- %1.Out going longwave radiation (W/m2) %地表发射长波辐射,与地表温度有关！ Out_going_long=zeros(m,n); Out_going_long=5.67e-8*Sur_Emiss_0.*LST_0.^4; geotiffwrite('ET-best\txt\2001161Out_going_long.tif',Out_going_long,R, 'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag); %----------------------------------------------------------------------------------------------------------------------- %2.income short wave radiation (W/m2) %入射短波辐射 Tsw=zeros(m,n); %atmospheric transmissivity 大气透过率与地形有关！ Rs=zeros(m,n); Tsw=0.75+2e-5.*DEM; Rs=Esun.*Tsw.*(1+0.033.*cos(DOY.*2.*PI./365)).*cos(Solar_Zenith_0.*PI./180); geotiffwrite('ET-best\txt\2001161Rs.tif',Rs,R, 'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag); %----------------------------------------------------------------------------------------------------------------------- %3.net income long wave radiation %入射长波辐射 Net_Income_Long=zeros(m,n); Net_Income_Long=Sur_Emiss_0.*(0.85.*(-1.*log(Tsw)).^0.09).*5.67e-8.*(Ta+273.16).^4; geotiffwrite('ET-best\txt\2001161Net_Income_Long.tif',Net_Income_Long,R, 'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag); %----------------------------------------------------------------------------------------------------------------------- %4.Net radiation (Rn) %净辐射通量 Rn=zeros(m,n); Rn=(1-albedo_0).*Rs+Net_Income_Long-Out_going_