[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
20040129: application for calculating sun and satellite angles (cont.)
- Subject: 20040129: application for calculating sun and satellite angles (cont.)
- Date: Tue, 03 Feb 2004 10:52:33 -0700
>From: "Miguel A. Martinez" <address@hidden>
>Organization: INM
>Keywords: 200401261012.i0QACDp2029502 McIDAS NAVCALC
Hi Miguel,
I apologize for not getting back to you sooner on this...
I am writing to echo Melanie Wetzel's thanks for the help you provided
her in getting access to solar and satellite viewing angles for her
radiative transfer modelling efforts. It is this kind of sharing that
makes the worldwide meteorological community such a joy to be a part
of!
Once again, thank you for the help you provided!
Cheers,
Tom
Date: Thu, 29 Jan 2004 08:18:56 -0800
From: Melanie Wetzel <address@hidden>
To: "Miguel A. Martinez" <address@hidden>
Subject: Re: 20040126: application for calculating sun and satellite angles
(cont.)
Thank you very much, Miguel, for taking the effort to prepare that
information in English and to send the details to me. I really
appreciate your advice and help.
Melanie Wetzel / Desert Research Institute / University of Nevada
Date: Thu, 29 Jan 2004 10:47:52 +0000
From: "Miguel A. Martinez" <address@hidden>
To: Tom Yoksas <address@hidden>
CC: Melanie Wetzel <address@hidden>
Subject: Re: 20040126: application for calculating sun and satellite angles
(cont.)
References: <address@hidden>
Dear Tom Yoksas:
I don't write before because i need to recover all the information and
translate it to english. I send to you this mail describing the zenith
and sun angles calculations. There are several parts:
a) the whole procedure to generate an AREA with zenith angles from the
NAVCALC output.
From McIDAS 7.5, the McIDAS command NAVCALC permit you to obtain an
ASCII output with several geometric parameters of an AREA file. The
problem is that the output of this command is ASCII with head and tail
comments. I tried to directly modify the NAVCALC command adding a simple
write to a file with the derired output only, but the program is in C
and very complicated to do this in few time. Then I made this easy
procedure:
a.1) From the operating system i redirected the ouput to a file, then i
filtered the lines with the time of the image (in this case an GVAR
AREA) and the file was then processed with PV-WAVE (you can made the
same with IDL):
navcalc.k MIDATA/ALL_AREA.nnnn FWD TYPE=I LINELE=3741 5993 LININC=10
189 ELEINC=10 324 > zenital_22_march_2000
cat zenital_22_march_2000|grep "11:46:00"|awk '{print $8}' >
zenital_22_march_2000.txt
(in this command:
MIDATA/ALL_AREA.nnnn is the dataset
TYPE=I this indicates that we work in AREA
coordinates
LINELE=3741 5993 this parameters are the IU coordinates
of the
image,
LININC=10 189 resolution (=10) of the AREA followed
of the
number of lines
ELEINC=10 324 resolution (=10) of the AREA followed
of the
number of elements)
The purpouse with these commands it to sincronize the output with the
source image. If you change in the awk the 8th position for another you
can obtain the latitude, longitude, sun height, etc. (See the NAVCALC
help for more details)
NOTE: One problem of this procedure is that if one pixel is not in
earth the NAVCALC not work.
a.2) In order to convert this ASCII file to AREA I had made in WAVE a
procedure as this (I have tried to locate the original but i don't
found):
The first part here is to read your file:
Read the ASCII file s=dc_read_free('zenital_22_march_2000.t',l)
reform to zen=reform(l, 324,189)
convert to byte zen=byte(zen)
The second part is to made the conversion to AREA:
One operation important is to obtain with IMGCOPY or IMGOPER an AREA in
the system that you are running the WAVE or IDL; the best is ti try
to
change the calibration to BRIT. Once you have the array in memory you
can use a procedure as the described bellow to save as AREA (this was a
version of a mail that I send to the mcidas-users mail list).
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
; first copy the AREA to local disc
; the position will be pointed to file AREAnumber
spawn,' imgcopy.k ADDE_SERVER/ADDE_NAME MIDATA/MI_NAME.position
SIZE=nlines nele ...'
;open the file AREA
; read the 64 long integer of directory block
name_area='$HOME/mcidas/data/AREA+string(number_area,format='(I4.4)')
openr,unit,name_area,/get_lun
header=lonarr(64)
readu,unit,header
header2=header
; change in header for example the number of SS, or coment, or ...
header2(...)=....
; also can change the number of lines, element (remenber to adjust them
with navigation)
; read the navigation codicil
size_nav=header(63-1)-header(35-1)
navegacion= lonarr( size_nav/4)
readu,unit,navegacion
; adjust the parameter with your header if necessary
navegacion(...)=.....
; read the calibration
calibracion=lonarr(512/4)
readu,unit,calibracion
; now open your final AREA
nombre_area=
strcompress('AREA'+string(final_area_number,format='(I4.4)'),/remove_all)
openw,unit,nombre_area,/get_lun
; write your header2
writeu,unit,header2
; write your navegacion2
writeu,unit,navegacion2
; write your calibration
writeu,unit,calibracion
; scale to integer or byte your data (e.g. your_array=byte( your_data)
;;;;;;;;;; here you your_array can be the zenith angle
; write your array
writeu,unit,your_array
;close the file
free_lun,unit
:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
Note 1: (also if you want more precision can use AREA with two bytes
and create or use another calibration)
Note 2: In the file "mail_describing_read_write_of_area.txt" you have
more details. It is a mail that i sent to mcidas-users mail list.
b) Our internal routine to calculate sun angles.
The name of the routine is altura.for. (altura means height in spanish).
It was written by a colleague called Mº Milagros Pertierra, she does not
work with us now. The algorithm is only an aproximation. Some comments
in spaniss are in the code. Any comments or modifications proposed by
you to the routine are wellcome.
c) A McIDAS program, zen.pgm, which calculates the zenith angle for a
METEOSAT image. I did write it to avoid the NAVCALC problem of the
pixels out of earth. The program zen.pgm generates an ASCII file with
the angle for one METEOSAT AREA (number 145). You can adapt to your own.
It use the aproach to use directly the navigation module. Also comments
or modifications are wellcome.
The AREA1146 is an AREA file with the zenith angle for METEOSAT.
Cheers
Miguel Angel
>
> >From: "Miguel A. Martinez" <address@hidden>
> >Organization: INM
> >Keywords: 200401261012.i0QACDp2029502 McIDAS NAVCALC
>
> Hola Miguel,
>
> >In order to calculate zenith angles for radiative transfer studies(I
> >work with the radiative transfer model RTTOV for MSG, MODIS and GOES
> >satellites), I have used the NAVCALC command.
>
> Thanks for pointing us in the right direction.
>
> >There is a mode you can obtain the zenith or view angle for a sector or
> >all image (read the help).
> >The output of NAVCALC command is ASCII; with some UNIX command (grep,
> >..) you can filter the header and tailer.
> >Once you have the data, then I had converted to AREA using a procedure
> >(PV-WAVE or IDL) that creates it with the data.
>
> OK. Melanie will probably end up doing something similar.
>
> >Also, we have internal program (from another people) to calculate the
> >sun angle.
>
> If it is not too much bother, can you send along the internal program
> and supporting modules?
>
> >If you need more information on zenith angle can mail to me.
>
> Your help is greatly appreciated!
>
> >Dr. Miguel A. Martinez Rubio
> > Jefe Unidad Satelites. Servicio de Teledeteccion
> > Instituto Nacional de Meteorologia (I.N.M.). Madrid. SPAIN
> > e-mail: address@hidden Tfno: +34 91 5819 662
> > Fax: +34 91 5819 846
>
> Cheers,
>
> Tom
--
Dr. Miguel A. Martinez Rubio
Jefe Unidad Satelites. Servicio de Teledeteccion
Instituto Nacional de Meteorologia (I.N.M.). Madrid. SPAIN
e-mail: address@hidden Tfno: +34 91 5819 662
Fax: +34 91 5819 846
--------------155D27BFA94494415E6829A6
Content-Type: text/plain; charset=us-ascii;
name="mail_describing_read_write_of_area.txt"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
filename="mail_describing_read_write_of_area.txt"
From - Thu Jan 29 09:29:25 2004
X-Mozilla-Status: 0001
X-Mozilla-Status2: 00000000
Message-ID: <address@hidden>
Date: Thu, 26 Jun 2003 16:16:47 +0000
From: "Miguel A. Martinez" <address@hidden>
Organization: inm.es
X-Mailer: Mozilla 4.8 [en] (X11; U; SunOS 5.8 sun4u)
X-Accept-Language: en
MIME-Version: 1.0
To: "address@hidden" <address@hidden>
Subject: Re: Matlab reader
References: <address@hidden>
Content-Type: text/plain; charset=us-ascii
Content-Transfer-Encoding: 7bit
Anthony James Wimmers wrote:
>
> Hi everyone,
>
> Does anybody have a matlab script that reads AREA, GRID or MD
> files (for instance, using the fread command).
>
> If not, I'll work on writing a matlab AREA file reader. Anyone interested
> in that?
>
> Tony Wimmers
> CIMSS/SSEC
I think that you have two posibilities:
a) Export the AREA to binary file with AXFORM.
This option has the advantage that at the same time you can read the
data calibrated and get also navigation. This example has been made with
PV-WAVE or IDL but I suppose similar could be made on Matlab.
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
; first copy the AREA to local disc
; the position will be pointed to file AREAnumber
spawn,' imgcopy.k ADDE_SERVER/ADDE_NAME MIDATA/MI_NAME.position
SIZE=nlines nele ...'
; call to AXFORM McIDAS program in order to write a binary file with the
data
; the most important is that you can select the calibration units
spawn,' axform.k number name_of_binary_file UNIT=your_unit'
; the binary files will be in your mcidas/data
cd,'$HOME/mcidas/data'
; then you must know the dimension of the array to create and the type
; the type could be (byte, integer or long integer) this depends on
calibration unit and the AREA
; this details are in the name_of_file.HDR
wv=intarr( nlines , nele)
openr,unit,'name_of_binary_file.001',/get_lun
readu,unit,brillo_wv
free_lun,unit
spawn,' lwu.k DEL name_of_binary_file.001'
spawn,' lwu.k DEL name_of_binary_file.HDR'
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
As you can see the type of UNIT and the sizes could be adjusted by you.
Also if you run AXFORM with NAV=YES (see help of AXFORM for details) you
can
read also the latitude and longitude data.
b) Read the file AREA directly.
The pseudocode would be:
:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
;open the file AREA
; read the 64 long integer of directory block
; read the number of lines, element and bytes by pixels from directory
; allocate an array with this dimension
; read or calculated the pointer to data array
; point the file to beginning of the data
; read the array
:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
This has the disadvantage that the data are RAW data and you must to
create your own calibration process.
The best of this is the possibility to use it to write one array on
McIDAS AREA format using an AREA as host.
For example you have used option-a and you have read the data and make
something with them, then you can write the processed data to McIDAS
AREA format.
In order to do this you could make something like this:
:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
;open the file AREA
; read the 64 long integer of directory block
name_area='$HOME/mcidas/data/AREA+string(number_area,format='(I4.4)')
openr,unit,name_area,/get_lun
header=lonarr(64)
readu,unit,header
header2=header
; change in header for example the number of SS, or coment, or ...
header2(...)=....
; also can change the number of lines, element (remenber to adjust them
with navigation)
; read the navigation codicil
size_nav=header(63-1)-header(35-1)
navegacion= lonarr( size_nav/4)
readu,unit,navegacion
; adjust the parameter with your header if necessary
navegacion(...)=.....
; read the calibration
calibracion=lonarr(512/4)
readu,unit,calibracion
; now open your final AREA
nombre_area=
strcompress('AREA'+string(final_area_number,format='(I4.4)'),/remove_all)
openw,unit,nombre_area,/get_lun
; write your header2
writeu,unit,header2
; write your navegacion2
writeu,unit,navegacion2
; write your calibration
writeu,unit,calibracion
; scale to integer or byte your data (e.g. your_array=byte( your_data)
; write your array
writeu,unit,your_array
;close the file
free_lun,unit
:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
--
Dr. Miguel A. Martinez Rubio
Jefe Unidad Satelites. Servicio de Teledeteccion
Instituto Nacional de Meteorologia (I.N.M.). Madrid. SPAIN
e-mail: address@hidden Tfno: +34 91 5819 662
Fax: +34 91 5819 846
--------------155D27BFA94494415E6829A6
Content-Type: text/plain; charset=us-ascii;
name="altura.for"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
filename="altura.for"
From - Thu Jan 29 10:20:42 2004
Return-Path: <address@hidden>
Received: from safdes.inm.es ([127.0.0.1]) by eolosmtp.inm.es
(Netscape Messaging Server 4.15) with ESMTP id HS8XDS00.KHT for
<address@hidden>; Thu, 29 Jan 2004 11:05:52 +0100
Received: (from pif@localhost)
by safdes.inm.es (8.11.6+Sun/8.11.6) id i0TA7cS29841
for address@hidden; Thu, 29 Jan 2004 10:07:38 GMT
Date: Thu, 29 Jan 2004 10:07:38 GMT
From: Pilar Fernandez <address@hidden>
Message-Id: <address@hidden>
Content-Type: text
X-Mozilla-Status: 8001
X-Mozilla-Status2: 00000000
X-UIDL: 11344-1013009125
C *** McIDAS Revision History ***
C* MRY MGA 21/2/90; CALC LA ALTURA SOLAR SOBRE EL HORIZONTE
C* ; ESTA EN EL USER
C 1 altura.for 21-Nov-95,9:00:00,`fpe' Version inicial McIDAS-X
C *** McIDAS Revision History ***
SUBROUTINE ALTURA(IDIA,HORA,LAT,LON,H,Z,COSZ)
C THIS IS INM PROPRIETARY SOFTWARE - ITS USE IS RESTRICTED.
C $ SUBROUTINE ALTURA(IDIA,HORA,LAT,LON,H,Z)
C $ MRY - CALCULA LA ALTURA SOLAR SOBRE EL HORIZONTE
C $ INPUT:
C $ IDIA = (I) DIA DE LA IMAGEN
C $ HORA = (R) HORA ASOCIADA A LA IMAGEN
C $ LAT = (R) LATITUD DEL PIXEL
C $ LON = (R) LONGITUD DEL PIXEL
C $ OUTPUT:
C $ H = (R) ALTURA DEL SOL SOBRE EL HORIZONTE PARA ESE PIXEL
C $ Z = (R) ANGULO CENITAL DEL SOL PARA ESE PIXEL
C EMPIEZA AQUI EL CODIGO NUEVO
C CALCULO DE LA ALTURA SOLAR SOBRE EL HORIZONTE
C***********************************************************************
REAL ALFA,DELTA,HORA,LON,LAT,ET,ANGSOL,SENOH,H,Z
PARAMETER (PI=3.141592,E=0.409280)
IDIAN=MOD(IDIA,1000)
JDIA=IDIAN
C***********************************************************************
C INTRODUCIR DIA DE 1 A 365
C INTRODUCIR HORA COMO HORA Y FRACCION DE HORA
C***********************************************************************
C***********************************************************************
C PRIMERO CALCULO DE LA ASCENSION RECTA
C DEPENDE DEL DIA
C***********************************************************************
ALFA=(JDIA-80.)*24./365.
C WRITE (6,*) 'ASCENSION RECTA ES', ALFA
ALFA=ALFA*15.*PI/180.
C***********************************************************************
C CALCULO DE LA DECLINACION
C DEPENDE DE LA ASCENSION RECTA
C**********************************************************************+
DELTA=ATAN(TAN(E)*SIN(ALFA))
C WRITE(6,*) 'LA DECLINACION ES', DELTA
C***********************************************************************
C CALCULO DE LA ECUACION DEL TIEMPO
C DEPENDE DEL DIA
C***********************************************************************
IF ((JDIA.GE.1).AND.(JDIA.LE.20)) THEN
ET=(15.+2.*JDIA)/5.
ELSEIF ((JDIA.GE.21).AND.(JDIA.LE.61)) THEN
ET=13.
ELSEIF ((JDIA.GE.62).AND.(JDIA.LE.110)) THEN
ET=(1382.-13.*JDIA)/48.
ELSEIF ((JDIA.GE.111).AND.(JDIA.LE.161)) THEN
ET=-2.
ELSEIF ((JDIA.GE.162).AND.(JDIA.LE.190)) THEN
ET=(3.*JDIA-500.)/14.
ELSEIF ((JDIA.GE.191).AND.(JDIA.LE.231)) THEN
ET=5.
ELSEIF ((JDIA.GE.232).AND.(JDIA.LE.290)) THEN
ET=(2204.-9*JDIA)/29.
ELSEIF ((JDIA.GE.291).AND.(JDIA.LE.331)) THEN
ET=-15.
ELSEIF ((JDIA.GE.332).AND.(JDIA.LE.365)) THEN
ET=(16.*JDIA-5741.)/33.
ENDIF
C***********************************************************************
C CALCULO DEL ANGULO HORARIO SOLAR
C DEPENDE DE LA LONGITUD,DE LA HORA TMG Y DEL TIEMPO DE BARRIDO
C EL TIEMPO DE BARRIDO DEPENDE DE LA LATITUD
C DEPENDE DE LA ECUACION DEL TIEMPO
C***********************************************************************
C ANGSOL=(HORA-12.-ET/60.-0.2)*15.+LON
ANGSOL=(HORA-12.-ET/60.-0.2)*15.-LON
C WRITE(6,*) 'EL ANGULO SOLAR ES',ANGSOL
ANGSOL=ANGSOL*PI/180.
C***********************************************************************
C CALCULO DE LA ALTURA SOLAR Y DEL ANGULO CENITAL
C DEPENDE DEL ANGULO HORARIO SOLAR, DE LA DECLINACION,
C DE LA LATITUD Y DE LA LONGITUD
C***********************************************************************
XXLAT=LAT*PI/180.
SENOH=SIN(DELTA)*SIN(XXLAT)+COS(DELTA)*COS(XXLAT)*COS(ANGSOL)
H=ASIN(SENOH)
H=H*180./PI
Z=90.-H
COSZ=COS(PI*Z/180.)
RETURN
END
--------------155D27BFA94494415E6829A6
Content-Type: image/x-portable-graymap;
name="zen.pgm"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
filename="zen.pgm"
subroutine main0
C MXCDSZ : LARGEST ALLOWED CODICIL SIZE
INTEGER MXCDSZ
PARAMETER (MXCDSZ = 5*128)
INTEGER DIR(64),IARR(MXCDSZ)
integer num_area,ival
character*12 cfi
integer lit
character*4 clit
real xlin,xele,xd
real xlon,xlat,xz
real xin(4),xout(3)
num_area=145
CALL READD(num_area,DIR)
NAVSIZ = DIR(63) - DIR(35)
IF (DIR(63).EQ.0) NAVSIZ = DIR(34) - DIR(35)
CALL ARAGET(num_area,DIR(35),NAVSIZ,IARR)
is=NVPREP(1,IARR)
IF(NV1INI(2,LIT('LL ')).NE.0) RETURN
do 1 i=3997,3997+(524-1)*40,4*10
do 2 j=2581,2581+(270-1)*40,4*10
xlin=j
xele=i
xd=0.0
is=nv1sae(xlin,xele,xd,xlat,xlon,xz)
c is=nv1eas(xlat,xlon,xd,xlin,xele,xz)
if (is .ge. 0 ) then
ifunc=lit('ANG ')
xin(1)=dir(4)
xin(2)=dir(5)
xin(3)=xlat
xin(4)=-xlon
is=nv1opt(ifunc,xin,xout)
else
xlat = 9999.0
xlon = 9999.0
xout(1) = 9999.0
endif
write (6,*) xlin,xele,xlat,xlon,xout(1)
2 continue
1 continue
return
end