brightness_temp.pro

从NOAA卫星传感器AVHRR探测数据计算获得大气亮温的程序源码

;+
; NAME:
;	brightness_temp
;
; PURPOSE:
;	calculate the brightness temperatures for a NOAA AVHRR image of an
;	IR channel (3,4 or 5)
; 	Brightness temperatures can be calculated from Planck's radiation equation
; 	using the slopes and intercept calibration values stored with each scan
; 	(cf. Polar Orbiting Users Guide Section 3.3 or ask Dominik)
;	
; CATEGORY:
;	NOAA satellite image processing
;
; CALLING SEQUENCE:
;	brightness_temp,image,channel,cal_coeff,satnr=satnr,temp=temp,eightbit=eightbit
;
; EXAMPLE:
;
; INPUTS:
;	image	: the GAC or LAC image produced with gac_channel.pro or lac_channel.pro
;		  Array(scansize,no_of_scans) of Pixel_Struct
;	channel : the channel number (must be 3,4, or 5)
;	cal_coeff : LonArr(2,5,no_of_scans)
;
; OPTIONAL INPUT PARAMETERS: 
;	satnr:	(BYTE) the number of the NOAA-satellite (11,12, or 14)
;
; KEYWORD INPUT PARAMETERS:
;	eightbit:	set this keyword if image is an 8-Bit resolution image only
;			(usually when it is a LAC image)
;
; OUTPUTS:
;	temp	: FltArr(scansize,no_of_scans): the brightness temperatures for each pixel
;		  of the image
;
; COMMON BLOCKS:
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;	
; PROCEDURE:
;	
; MODIFICATION HISTORY:
;	All information describing the format
;	of GAC data can be found on WWW in
;	<A HREF="http://www2.ncdc.noaa.gov/docs/podug/">Polar Orbiter Guide</A>
;	First implementation Oct, 25, 1997 by Dominik Brunner, ETH Zurich
;-

pro brightness_temp,image,channel,cal_coeff,satnr=satnr,temp=temp,eightbit=eightbit

print,'Calculate brightness temperature for an IF channel'

; check for input parameters
IF (n_elements(image) EQ 0) OR (n_elements(channel) EQ 0) OR $
   (n_elements(cal_coeff) EQ 0) THEN return

IF n_elements(image(0,*)) NE n_elements(cal_coeff(0,0,*)) THEN BEGIN
   print,'parameter cal_coeff has wrong size!'
   return
ENDIF

IF channel LT 3 THEN BEGIN
   print,'no brightness temperature can be calculated for this channel'
   return
ENDIF

; enter the satellite number (because individual calibration factors aply)
IF n_elements(satnr) EQ 0 THEN BEGIN
   print,''
   read,'Enter Satellite number (11,12, or 14): ',satnr
   IF NOT ((satnr EQ 11) OR (satnr EQ 12) OR (satnr EQ 14)) THEN BEGIN
   	print,'no correction factors for this satellite number entered in the table yet'
  	print,'therefore the corrections for satellite number 14 will be applied'
   	satnr=14
   ENDIF
ENDIF
CASE satnr OF
     11: satno=0
     12: satno=1
     14: satno=2
ENDCASE
print,''

; enter the temperature range which should be fit best
REPEAT BEGIN
  print,'select the approx. temperature range to fit the best'
  print,'range 1: 190-230 K (CB clouds up to tropopause)'
  print,'range 2: 230-270 K (mid tropospheric clouds and cold earth)
  print,'range 3: 270-310 K (sea surface temperatures)
  print,'range 4: 290-330 K (warm land surface)
  read,'Enter your selection (1-4): ',range
  print,''
END UNTIL ((range GT 0) AND (range LT 5))
range=range-1

; define the satellite specific central wavenumbers (in cm-1) for each channel
; the central wavelengths are different for the four temperature ranges
; range1=190-230K, range2=230-270K, range3=270-310K (sea surface), range4=290-330K
      
cwnr=FltArr(3*4,3)	; (3 channels x 4 temp. ranges per channel) x 3 satellites
cwnr(*,0)=[2663.50,2668.15,2671.40,2670.96,$ 		; channel 3, satnr 11
	   926.81,927.36,927.83,927.75,$		; channel 4, satnr 11
	   841.40,841.81,842.20,842.14]			; channel 5, satnr 11
cwnr(*,1)=[2632.713,2636.669,2639.61,2640.817,$ 	; channel 3, satnr 12
	   920.0158,920.5504,921.0291,921.2741,$	; channel 4, satnr 12
	   836.6847,837.0251,837.3641,837.5612]		; channel 5, satnr 12
cwnr(*,2)=[2638.652,2642.807,2645.899,2647.169,$ 	; channel 3, satnr 14
	   928.2603,928.8284,929.3323,929.5878,$	; channel 4, satnr 14
	   834.4496,834.8066,835.1647,835.374]		; channel 5, satnr 14

; define central wavenumber for current selections
cwnr=cwnr((channel-3)*4+range,satno)
print,'central wavenumber for channel'+STRING(channel,format='(I2)')+$
      ' on NOAA-'+STRING(satnr,format='(I2)')+' is'+strcompress(cwnr)

; define constants for planck's law
c1=1.1910659E-05	; in mW m-2 steradian cm-4
c2=1.438833D		; cm K-1

; define scaling factors for calibration coefficients

; Channel counts must be multiplied by a factor of 4 in eightbit version
; Therefore the scale factor is a factor of 4 lower than usual

IF keyword_set(eightbit) THEN scale1=2D^28 ELSE scale1=2D^30
			; for slope
scale2=2D^22		; for intercept

temp=FltArr(n_elements(image(*,0)),n_elements(image(0,*)))
FOR i=0,n_elements(image(0,*))-1 DO BEGIN
    radiance=Double(image(*,i).value)*cal_coeff(0,channel-1,i)/scale1+$
	     cal_coeff(1,channel-1,i)/scale2
    temp(*,i)=c2*cwnr/alog(1+c1*cwnr^3/radiance)
ENDFOR

index=WHERE(temp EQ 0,count)
IF count GT 0 THEN temp(index)=!values.f_nan

end