[
Date Prev][
Date Next][
Thread Prev][
Thread Next][
Date Index][
Thread Index]
20030120: McIDAS-X user command at Santa Clara University (cont.)
- Subject: 20030120: McIDAS-X user command at Santa Clara University (cont.)
- Date: Tue, 21 Jan 2003 11:25:49 -0700
>From: "Mahalingam Haritharan" <address@hidden>
>Organization: Santa Clara University
>Keywords: 200301211713.h0LHD0x09375 McIDAS
Hari,
>Here is my code in fortran to take location.
OK, I took a quick look at the code and see that it is always
requesting a chunk of data that falls within the first 480x640 pixels
from the input image. The code that does this is:
>c It wants a certen number of sleects else it wont do a mcalin
>c Looks like mcaget wants at lest these two slectts
>c Set selection conditions
> selects(1)='SIZE 480 640'
> selects(2)='BAND 4'
So, without looking any harder at your code, I would say that you can
modify the code to include the entire image by changing the SIZE select
condition. For instance, if the input image has 600 lines and 800
elements, change selects(1) to:
selects(1)='SIZE 600 800'
After making the change, of course, you will need to rebuild your
application.
You could even make your application more flexible by allowing the user
to input the size of the LINxELE box that is to be requested. You can
follow the example in the code that gets command line input for the
AREA file numbers. For instance:
C..... Expanded command line syntax:
C
C JBLACK area1 area2 outfile nline nelemet
C
...
integer ib, ie, nch
integer nlin, nele
integer nchars
...
nlin = ipp( 4, 480 ) ! default is 480
nele = ipp( 5, 640 ) ! default is 650
write( cstr, '("SIZE ",I4,1X,I4)' ) nlin, nele
nch = nchars( cstr, ib, ie)
...
selects(1) = cstr(ib:ie)
I hope that this helps.
Tom
> subroutine main0
>C
>C PROGRAM JBLACK
>C
>C ? JBLACK -- Simple feature tracking algorithm for support of lightning
>C ? detection by satellite.
>C ? JBLACK area1 area2 outfile
>C ? Parameters:
>C ? area1 | First of two McIDAS area files
>C ? area2 | Second of two McIDAS area files
>C ? outfile | Output file
>C ?
>C ? Remarks:
>C ? Two area files will be opened and data extracted. Information
>C ? will be written to outfile.
>C ?
>C ? Goes in /home/user/mcideas/src. To compile, "fx jblack l"
>C ? Add this comand(file name of this porgram) to IR404.CMD so that the
>McIDEAS sceduler
>C ? will preform this program as soon as it updateds the IR4 images.
>C ? ---------
>
>c Variable declaration
> character*12 dataset ! Data set name
> character*80 selects(2) ! Array holding values for area selecti
> on
> character*2 format ! Format of the returned data
> character*4 unit ! Units of the returned data
> character*12 cpp ! McIDAS function
> character*12 outfile ! Output file
> character*4 cnvarea ! Character equivalent of nvarea
> character*4 cinarea ! Character equivalent of iadir(1)
> character*120 cstr ! Character string for usin in skeyin
> character*24 etime
> character*24 epochtime
>
> integer nselects ! Number of elements in selects
> integer iadir(64) ! Array for the area directory
> integer maxbyte ! Maximum byte size per line
> integer msgflag ! Message flag
> integer handle ! Data stream indicator
> integer status ! Success/failure to call to mcaget/mca
> lin
> integer data_buff(640) ! Data array
> integer nvarea ! Area number for use in nvset
> integer lbeg ! Upper image line in area
> integer ebeg ! Left image element in area
> integer lres ! Line resolution
> integer eres ! Element resolution
> integer ilin ! Image line
> integer iele ! Image element
> integer work_array(480,640) ! Testing array
> integer threshold ! tempeture below wich ligting
> will occrue, in dk
> integer rows, cols ! start row, start colum varbils
> integer temp_array(480,640) ! Holds the tempture
> integer ccount ! cool count
> integer z
> integer altidue ! Arbtray number of meters asuming the
> strom is
>above the gound
>
> real xlin ! Real equivalent of ilin
> real xele ! Real equivalent of iele
> real xlat ! Latitude
> real xlon ! Longitude
> real xdummy ! Dummy variable
>
> threshold = 2000
altidue = 2000
>
>c=======================================================================
>
>c Turn on scudler to get updated images. Makes sure is no suspned on SKL
>list, asusming is #5 on list.
> Call skeyin("SF 2")
> Call skeyin("SKU ON")
> Call skeyin("SKU REL 5")
>
>c Example. Could suspend other lines in sculer
>c call skeyin("SKU SUSPEND 1")
>
>c Displays the picuter the program is curnching. Comint out if it gets
>annoying
> Call skeyin("IMGDISP LIDIR404 2")
> Call skeyin("MAP VH")
> Call skeyin("EU REST SVGAIR")
>
>c Get positional parameter
> area1=ipp(1,9999)
> area2=ipp(2,9999)
> outfile=cpp(3,'jblack4.out')
>
>c Open output file
> open(unit=11,file=outfile,status='replace')
>
>c Set the dataset
> dataset='LIDIR404'
>
>c It wants a certen number of sleects else it wont do a mcalin
>c Looks like mcaget wants at lest these two slectts
>c Set selection conditions
> selects(1)='SIZE 480 640'
> selects(2)='BAND 4'
>c selects(3)='DAY 02254'
>c selects(4)='TIME 20:59 21:01'
>c selects(3)='POS 1'
>
>
>c Set number of selection conditions
> nselects=2
>
>c Set the units of the returned data
> unit='TEMP'
>
>c Set the format of the returned data buffer
> format='I4'
>
>c Set the maximum number of bytes per line
> maxbyte=9000000
>
>c Set the message flag
> msgflag=1
>
>c Open a connection to area file
> status=mcaget(dataset,nselects,selects,unit,format,maxbyte,
> & msgflag,iadir,handle)
>
> if(status.lt.0) then
> call edest('mcaget failed',0)
> return
> endif
> do 5 i=1,64
>c write(11,*) iadir(i)
>5 continue
>
>
>
>
>c Retrieve navigation information from file. Here an imgcopy must be
>performed
>c in order to designate an area number. This step is necessary as
>apparently the
>c subroutine nvset has not been updated to adde form.
>
>c Note: is a reltve postinl number, could be biger than 9. Put IF to chec
> k
> IF(iadir(1) > 9) then
> write(cinarea, '(i2)') iadir(1)
> else
> write(cinarea,'(i1)') iadir(1)
> endif
>
> nvarea=9999
> write(cnvarea,'(i4)') nvarea
>
> cstr='IMGCOPY LIDIR404.'//cinarea//' LID.'//cnvarea//' SIZE=ALL'
>
> call skeyin(cstr)
>
>c Nvset gets navigation information
> status=nvset('AREA',nvarea)
> if(status.lt.0) then
> call edest('nvset failed',0)
> return
> endif
>
>c Navigation is performed from image coordinates via nv1sae. So we need
>c to convert from area coordinates to image coordinates
>
> lbeg=iadir(6)
> lres=iadir(12)
> ebeg=iadir(7)
> eres=iadir(13)
>
> DO i=0,479
>
> status=mcalin(handle,data_buff)
> if(status.lt.0) then
> call edest('mcalin failed',0)
> return
> endif
>
> ilin=i*lres+lbeg
>
> do j=0,639
> iele=j*eres+ebeg
> xlin=real(ilin)
> xele=real(iele)
> status=nv1sae(xlin,xele,0.,xlat,xlon,xdummy)
> if(status.ne.0) then
>c write(11,*) ilin,iele,data_buff(j+1)
>
> else
>c write(11,*) ilin,iele,data_buff(j+1),xlat,xlon
> endif
>
> temp_array(i+1,j+1) = data_buff(j+1)
>
> If(data_buff(j+1) < threshold) then
> work_array(i+1,j+1)=1
> else
> work_array(i+1,j+1)=0
> endif
>
>
> Enddo
>
>
> EndDO
>
>
>
>
>c Dived up the image into 60x64 grids, wich crates a 16/20 grid. Any grid
>
>with 75% or more cold pixles is a targit for the satlite.
>c The satlite targit point will be the middel of the grid. This is crude
> and
>fast
>c First Linse of the Ephermeris data file, Numberof EphemrisPoints is
>assuming one point per minit over 2 hours
>c EphemerisLLATimePos means put the data in this format <time form epoc
>xxxx.xxx>\t <Latiude> \t<Longidue> \t<altidue>
>c I have arbitraly set alditued to 2000 meaters. Assuming an avrage hight
>
>for stomes. If its not good, change it.
>
>
> write(11,*) "stk.v.4.3"
> write(11,*) "\n BEGIN Ephemeris"
> write(11,*) "ScenarioEpoch ",epochtime(iadir(4),iadir(5
> ))
> write(11,*) "NumberofEphemerisPoints 3000"
> write(11,*) "CentralBody Earth"
> write(11,*) "EphemerisLLATimePos",'\n'
> Do rows=0,479,200
> Do cols=0,639,150
> ccount=0
>
> Do i=0,199
> do j=0,149
> If( work_array(i+rows+1, j+cols+1) > 0) th
> en
> ccount=ccount+1
> endif
> enddo
> enddo
>
>
> if(ccount > 200 ) then
> ilin=rows+(16/2)+lbeg
> iele=cols+(32/2)+ebeg
> xlin=real(ilin)
> xele=real(iele)
> DO z=0,100
> status=nv1sae(xlin,xele,0.,xlat,xlon,xd
> ummy)
> if (status.ne.0) then
> write(11,*) 'T', ccount, temp_a
> rray(rows+1, cols+1), ilin, elin
> else
> write(11,*) z,"0.00",'\t', xlat
> ,'\t',xlon,'\t',altidue
>
> endif
>
> endDO
>
>
> endif
> enddo
>
> enddo
> write(11,*) "END Ephemeris",'\r'
>
> close(unit=11)
> call edest('DONE',0)
> return
> end
>
>
>
>c To convert the time of image from macidas to the format stk whants. int to
>sting
>c form julin year and day: ccyyddd, and nomil time hhmmss to /tday month
>year time = 1 Jan 1999 00:00:00.00
>c since this is writn in 2002, all dates are asume do be 20yy year format.
>
> Function epochtime(oldjulen, oldtime)
>
> character*24 epochtime
> integer oldjulen
> integer oldtime
> integer ascii_off_set
> integer year
> integer day
> integer month
> integer hour
> integer min
> integer second
> character*3 smonth
> integer status
> character*2 sday
> character*2 ssecond
> character*2 smin
> character*2 shour
> character*4 syear
>
> ascii_off_set = 48
>
>!Mccydtodmy = McIdas ccyydd(julin) to day month year
>
> status = Mccydtodmy (oldjulen, day, month, year)
> if(status.lt.0) then
> call edest('Mccydtodmy failed',0)
> endif
>
> hour = oldtime/10000
>
> min = mod(oldtime,10000) / 100
> second = mod(oldtime,100)
>
>! I dont think ther is a swich operateor in fortran 77, so
>
> if(month == 1) then
> smonth='Jan'
> else if(month==2) then
> smonth='Feb'
> else if(month==3) then
> smonth = 'Mar'
> else if(month==4) then
> smonth='Apr'
> else if(month==5) then
> smonth='May'
> else if(month==6) then
> smonth='Jun'
> else if(month==7) then
> smonth='Jul'
> else if(month==8) then
> smonth='Aug'
> else if(month==9) then
> smonth='Sep'
> else if(month==10) then
> smonth='Oct'
> else if(month==11) then
> smonth='Nov'
> else if(month==12) then
> smonth='Dec'
> endif
>
> if(day > 9) then
> sday(1:1)='1'
> else
> sday(1:1) = ' '
> endif
>
> sday(2:2)= CHAR(mod(day,10)+ascii_off_set)
>
> syear(1:1)=CHAR(year/1000+ascii_off_set)
> syear(2:2)=CHAR( MOD(year,1000)/100+ascii_off_set)
> syear(3:3)=CHAR( MOD(year,100)/10 +ascii_off_set)
> syear(4:4)=CHAR( MOD(year,10)+ascii_off_set)
>
> shour(1:1)=CHAR( hour/10+ascii_off_set)
> shour(2:2)=CHAR( mod(hour,10)+ascii_off_set)
>
> smin(1:1)=CHAR(min/10+ascii_off_set)
> smin(2:2)=CHAR( mod(min,10)+ascii_off_set)
>
> ssecond(1:1)=Char(second/10+ascii_off_set)
> ssecond(2:2)=CHAR( mod(second,10)+ascii_off_set)
>
> epochtime= '\t'//sday//' '//smonth//' '//syear//'
>'//shour//':'//smin//':'//ssecond//'.00'
>
> return
> END
> END
>
>------=_NextPart_000_615e_3f3_780--
>
