Subtracting Images

[dl08]

I was testing out the method for requesting data from the server with the desired lat/lon coordinates. I wrote this script that requests the image data at a specific lat/lon location, (29.7276,-85.0274) for Ch2 and Ch4 of each imagetime. (My dataset consists of one imagetime per day so this script loops through the first 6 days of the dataset) :

Code: Select all

datapath1=('C:\\Users\\Daphne\\collocation\\flgoes_resampled\\CH2')
datapath2=('C:\\Users\\Daphne\\collocation\\flgoes_resampled\\CH4')

##
##Define a local ADDE dataset
##
makeLocalADDEEntry(dataset='FL', imageType='CH2', mask=datapath1, format='McIDAS Area', save=True)
makeLocalADDEEntry(dataset='FL', imageType='CH4', mask=datapath2, format='McIDAS Area', save=True)

ch2data=getLocalADDEEntry(dataset='FL', imageType='CH2')
ch4data=getLocalADDEEntry(dataset='FL', imageType='CH4')


beginday='2012001'
endday='2012006'


##
##Define an ADDE request to create a list of images for a range of days
##
ADDE_listRequest = dict(
    localEntry=ch4data,
    descriptor=ch4data,
    dataset='FL',
    day=(beginday,endday),
    position='ALL'
)
dateTimeList=listADDEImageTimes(**ADDE_listRequest)
#
##getADDEImage for each of the directories returned
##
for dateTime in dateTimeList:

   imageDate=dateTime['day']
   
   ADDE_T2_getRequest = dict(
      localEntry=ch2data,
      descriptor=ch2data,
      day=(imageDate,imageDate),
      coordinateSystem=LATLON,
        location=(29.7276,-85.0274),
        place=ULEFT,
      band=2,
      unit='TEMP',
        size=(1,1),
   )
   T2MetaData,T2Data=getADDEImage(**ADDE_T2_getRequest)
   
   ADDE_T4_getRequest = dict(
        localEntry=ch4data,
        descriptor=ch4data,
        day=(imageDate,imageDate),
        coordinateSystem=LATLON,
        location=(29.7276,-85.0274),
        place=ULEFT,
        band=4,
        unit='TEMP',
        size=(1,1),
   )
   T4MetaData,T4Data=getADDEImage(**ADDE_T4_getRequest)
   print 'finished get calls for date %s' % imageDate

##
##getADDEImage returned data, so apply algorithm 
##

   ##use function below for fogproduct in TD values
   IRFog=sub(T4Data,T2Data)   
   
print 'finished data calls'

I then get this returned:

Code: Select all

adde://localhost/imagedirectory?&PORT=8113&COMPRESS=gzip&USER=idv&PROJ=0&VERSION=1&DEBUG=false&TRACE=0&GROUP=FL&DESCRIPTOR=8A7C7E34&BAND=ALL&DAY=2012001 2012006&POS=ALL
adde://localhost/imagedata?&PORT=8113&COMPRESS=gzip&USER=idv&PROJ=0&VERSION=1&DEBUG=false&TRACE=0&GROUP=FL&DESCRIPTOR=8712C574&BAND=2&LATLON=29.7276 -85.0274 E&PLACE=ULEFT&SIZE=1 1&UNIT=TEMP&MAG=1 1&SPAC=4&NAV=X&AUX=YES&DOC=X&DAY=2012001 2012001&TIME=&POS=0&TRACK=0
adde://localhost/imagedata?&PORT=8113&COMPRESS=gzip&USER=idv&PROJ=0&VERSION=1&DEBUG=false&TRACE=0&GROUP=FL&DESCRIPTOR=8A7C7E34&BAND=4&LATLON=29.7276 -85.0274 E&PLACE=ULEFT&SIZE=1 1&UNIT=TEMP&MAG=1 1&SPAC=4&NAV=X&AUX=YES&DOC=X&DAY=2012001 2012001&TIME=&POS=0&TRACK=0
finished get calls for date 2012001
adde://localhost/imagedata?&PORT=8113&COMPRESS=gzip&USER=idv&PROJ=0&VERSION=1&DEBUG=false&TRACE=0&GROUP=FL&DESCRIPTOR=8712C574&BAND=2&LATLON=29.7276 -85.0274 E&PLACE=ULEFT&SIZE=1 1&UNIT=TEMP&MAG=1 1&SPAC=4&NAV=X&AUX=YES&DOC=X&DAY=2012002 2012002&TIME=&POS=0&TRACK=0
adde://localhost/imagedata?&PORT=8113&COMPRESS=gzip&USER=idv&PROJ=0&VERSION=1&DEBUG=false&TRACE=0&GROUP=FL&DESCRIPTOR=8A7C7E34&BAND=4&LATLON=29.7276 -85.0274 E&PLACE=ULEFT&SIZE=1 1&UNIT=TEMP&MAG=1 1&SPAC=4&NAV=X&AUX=YES&DOC=X&DAY=2012002 2012002&TIME=&POS=0&TRACK=0
finished get calls for date 2012002
adde://localhost/imagedata?&PORT=8113&COMPRESS=gzip&USER=idv&PROJ=0&VERSION=1&DEBUG=false&TRACE=0&GROUP=FL&DESCRIPTOR=8712C574&BAND=2&LATLON=29.7276 -85.0274 E&PLACE=ULEFT&SIZE=1 1&UNIT=TEMP&MAG=1 1&SPAC=4&NAV=X&AUX=YES&DOC=X&DAY=2012003 2012003&TIME=&POS=0&TRACK=0
adde://localhost/imagedata?&PORT=8113&COMPRESS=gzip&USER=idv&PROJ=0&VERSION=1&DEBUG=false&TRACE=0&GROUP=FL&DESCRIPTOR=8A7C7E34&BAND=4&LATLON=29.7276 -85.0274 E&PLACE=ULEFT&SIZE=1 1&UNIT=TEMP&MAG=1 1&SPAC=4&NAV=X&AUX=YES&DOC=X&DAY=2012003 2012003&TIME=&POS=0&TRACK=0
finished get calls for date 2012003
adde://localhost/imagedata?&PORT=8113&COMPRESS=gzip&USER=idv&PROJ=0&VERSION=1&DEBUG=false&TRACE=0&GROUP=FL&DESCRIPTOR=8712C574&BAND=2&LATLON=29.7276 -85.0274 E&PLACE=ULEFT&SIZE=1 1&UNIT=TEMP&MAG=1 1&SPAC=4&NAV=X&AUX=YES&DOC=X&DAY=2012004 2012004&TIME=&POS=0&TRACK=0
adde://localhost/imagedata?&PORT=8113&COMPRESS=gzip&USER=idv&PROJ=0&VERSION=1&DEBUG=false&TRACE=0&GROUP=FL&DESCRIPTOR=8A7C7E34&BAND=4&LATLON=29.7276 -85.0274 E&PLACE=ULEFT&SIZE=1 1&UNIT=TEMP&MAG=1 1&SPAC=4&NAV=X&AUX=YES&DOC=X&DAY=2012004 2012004&TIME=&POS=0&TRACK=0
finished get calls for date 2012004
adde://localhost/imagedata?&PORT=8113&COMPRESS=gzip&USER=idv&PROJ=0&VERSION=1&DEBUG=false&TRACE=0&GROUP=FL&DESCRIPTOR=8712C574&BAND=2&LATLON=29.7276 -85.0274 E&PLACE=ULEFT&SIZE=1 1&UNIT=TEMP&MAG=1 1&SPAC=4&NAV=X&AUX=YES&DOC=X&DAY=2012005 2012005&TIME=&POS=0&TRACK=0
adde://localhost/imagedata?&PORT=8113&COMPRESS=gzip&USER=idv&PROJ=0&VERSION=1&DEBUG=false&TRACE=0&GROUP=FL&DESCRIPTOR=8A7C7E34&BAND=4&LATLON=29.7276 -85.0274 E&PLACE=ULEFT&SIZE=1 1&UNIT=TEMP&MAG=1 1&SPAC=4&NAV=X&AUX=YES&DOC=X&DAY=2012005 2012005&TIME=&POS=0&TRACK=0
finished get calls for date 2012005
adde://localhost/imagedata?&PORT=8113&COMPRESS=gzip&USER=idv&PROJ=0&VERSION=1&DEBUG=false&TRACE=0&GROUP=FL&DESCRIPTOR=8712C574&BAND=2&LATLON=29.7276 -85.0274 E&PLACE=ULEFT&SIZE=1 1&UNIT=TEMP&MAG=1 1&SPAC=4&NAV=X&AUX=YES&DOC=X&DAY=2012006 2012006&TIME=&POS=0&TRACK=0
adde://localhost/imagedata?&PORT=8113&COMPRESS=gzip&USER=idv&PROJ=0&VERSION=1&DEBUG=false&TRACE=0&GROUP=FL&DESCRIPTOR=8A7C7E34&BAND=4&LATLON=29.7276 -85.0274 E&PLACE=ULEFT&SIZE=1 1&UNIT=TEMP&MAG=1 1&SPAC=4&NAV=X&AUX=YES&DOC=X&DAY=2012006 2012006&TIME=&POS=0&TRACK=0
finished get calls for date 2012006
finished data calls

I then use the getData() method to return the pixel values at this location: LATLON=29.7276 -85.0274
For example, for CH2 of day 2012002 I use this request string:

Code: Select all

from edu.wisc.ssec.mcidas import AreaFile as af
url="adde://localhost/imagedata?&PORT=8113&COMPRESS=gzip&USER=idv&PROJ=0&VERSION=1&DEBUG=false&TRACE=0&GROUP=FL&DESCRIPTOR=8712C574&BAND=2&LATLON=29.7276 -85.0274 E&PLACE=ULEFT&SIZE=1 1&UNIT=TEMP&MAG=1 1&SPAC=4&NAV=X&AUX=YES&DOC=X&DAY=2012002 2012002&TIME=&POS=0&TRACK=0"
a=af(url);
value=a.getData();
print value

and get this returned:

Code: Select all

array([[I, [array([I, [array('i', [2808, 2819, 2830, 2840])])])

and for CH4 I use this request string:

Code: Select all

from edu.wisc.ssec.mcidas import AreaFile as af
url="adde://localhost/imagedata?&PORT=8113&COMPRESS=gzip&USER=idv&PROJ=0&VERSION=1&DEBUG=false&TRACE=0&GROUP=FL&DESCRIPTOR=8A7C7E34&BAND=4&LATLON=29.7276 -85.0274 E&PLACE=ULEFT&SIZE=1 1&UNIT=TEMP&MAG=1 1&SPAC=4&NAV=X&AUX=YES&DOC=X&DAY=2012002 2012002&TIME=&POS=0&TRACK=0"
a=af(url);
value=a.getData();
print value

and get this returned:

Code: Select all

array([[I, [array([I, [array('i', [2826, 2833, 2841, 2853])])])

I have a couple of questions, I apologize if they are rather trivial but I am new to Jython and working with this type of data.
Why am I getting an array of 4 values for a lat/lon location? For a given lat/lon coordinate, shouldn't I be getting one corresponding pixel value?

Also, I want to subtract Ch4-Ch2 temperature values and get the pixel difference value for a corresponding lat/lon coordinate. Where does the "IRFog=sub(T4Data,T2Data)" come into play in the script? Or would I need to somehow subtract the two arrays I am getting returned?

[tomw]

Why am I getting an array of 4 values for a lat/lon location? For a given lat/lon coordinate, shouldn't I be getting one corresponding pixel value?

ADDE servers always return the number of pixels across a line as a multiple of 4. I cannot tell you why, just that this is the case. The first data value should correspond to the requested latitude-longitude.

Also, I want to subtract Ch4-Ch2 temperature values and get the pixel difference value for a corresponding lat/lon coordinate. Where does the "IRFog=sub(T4Data,T2Data)" come into play in the script? Or would I need to somehow subtract the two arrays I am getting returned?

The AreaFile.getData() method returns a 3-dimensional, integer array (bands, lines, elements) of the un-scaled data values. It appears that your data are temperatures (K) scaled by 10.

To get the difference, you would just subtract the elements...something like this:

Code: Select all

diff = value4[0][0][0] - value2[0][0][0];

assuming that "value4" is the "getData()" from channel 4, and likewise for channel 2.

[abhiram]

Tom,

Thanks for the clarification. We were able to figure out the same using the data-probe tool present in McIDAS-V where its readout was matching with the first value present in the array. To compensate for the scaling, I divide the value by 10.

Thank you for helping with my questions, I now have a script that can compute the temperature difference for each day of the year for a set of weather stations. Will let you know if I have any more questions.

Happy Thanksgiving

[abhiram]

Hello,

I am running into this error, when I try to compute the temperature difference for all the 366 days

Code: Select all

22:32:33.397 [MainThread] INFO  jython - print: java.lang.OutOfMemoryError: java.lang.OutOfMemoryError: Java heap space

The script is running on a i5 laptop with 4 gigs of ram. I am not sure if that's sufficient or not. One approach that I have tried is calling the collectGarbage routine at the beginning of the loop. Is it as simple as calling the routine?

I have attached my script and error log. Tomorrow, I am planning on running the script on a better machine. Will let you know how that goes.

Thanks!

[tomw]

The only thing I can suggest is that when you run this, you set the maximum heap size on the command line. For a 64-bit machine, something like:

java - Xmx3000m ......

for 3GB. If you're on a 32-bit machine, you cannot specify a value larger than about 1500m, regardless of how much more memory you might have.

As for the garbage collection or other considerations, I must defer to the programmers who wrote the routines you are using...sorry...