NOTICE: This version of the NSF Unidata web site (archive.unidata.ucar.edu) is no longer being updated.
Current content can be found at unidata.ucar.edu.
To learn about what's going on, see About the Archive Site.
>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