Lookup table

[Thardy]

Dear all,
I am trying to read and write data from a lookup table(binary file LUT) based on geographical coordinates from binary files "LON"(longitude) and "LAT"(latitude). The program is running but does not provide real values after the firt reading(0.013231), i.e. the remaining values are 0.00...
I think from line 36 ("get LUT data") of the program below, the reading is problematic. Also in the line "if(ABS(latit(x,y)-(knt-1)*DIF_LATG+90.0", I can not use "latit" in replacement of "latit(x,y)" as there would be a compiler error asking for a scaler.
I am surely missing something. Any help will be appreciated:

program LUT
implicit none
integer*4:: lat, lon,x,y
INTEGER, PARAMETER :: NBYTE_G=4 !byte of albedo data
INTEGER, PARAMETER :: NLONG_G=2880, NLAT_G=1441 !size of albed data
REAL, PARAMETER :: DIF_LONGG=0.125, DIF_LATG=0.125 !size of each pixel of albedo data
REAL*4 :: gal
REAL(KIND=NBYTE_G), ALLOCATABLE :: galraw,
REAL(KIND=4)::latit(1236,1656),longit(1236,1656)
REAL(KIND=4)::longdi(1236,1656)
INTEGER :: iost, knt, knt2
INTEGER :: ilat, ilong
INTEGER :: inlon, inlat
integer, parameter ::irec1=1236*1656*4
C-------
open(10,FILE='LON',FORM='UNFORMATTED',
& ACCESS='DIRECT',RECL=irec1)
read(10, rec=1) longit
open(20,FILE='LAT',FORM='UNFORMATTED',
& ACCESS='DIRECT',RECL=irec1)
read(20, rec=1) latit

allocate(galraw(NLAT_G, NLONG_G))
open(30,FILE='LUT',FORM='UNFORMATTED',STATUS='OLD'
& ,ACCESS='DIRECT',ACTION='READ',RECL=NLONG_G*4)

do knt=1, NLAT_G
read(30,REC=knt) (galraw(knt,knt2)
& , knt2=1, NLONG_G)

C---------open output files
open(40,FILE='out1',FORM='UNFORMATTED',STATUS='UNKNOWN'
& ,ACCESS='DIRECT',ACTION='WRITE',RECL=irec1)
open(50,FILE='out2',FORM='FORMATTED')
C----------- get LUT data
if(ABS(latit(x,y)-(knt-1)*DIF_LATG+90.0)
& .le. DIF_LATG*0.5) then
ilat=knt

exit
endif
enddo

do knt=1, NLONG_G
longdi(x,y)=(knt-1)*DIF_LONGG
if(longit(x,y) .gt. 0) then
if(longdi(x,y) .gt. 180.0) then
longdi(x,y)=longdi(x,y)-360.0
endif
else if(longit(x,y) .lt. 0.0) then
if(longdi(x,y) .ge. 180.0) then
longdi(x,y)=longdi(x,y)-360.0
endif
endif
if(ABS(longit(x,y)-longdi(x,y)) .le. DIF_LATG*0.5) then
ilong=knt
exit
endif
end do
gal=galraw(ilat, ilong)
write(40,rec=1) gal
write(50,'(1236F12.6)') gal

deallocate(galraw)
stop
End program LUT

 

[Thardy]

LON file link:

 

[Thardy]

LAT file link:

 

[gummibaer]

Thardy,

you should rethink the structure of your program:

do knt=1, NLAT_G
read(30,REC=knt) (galraw(knt,knt2)
& , knt2=1, NLONG_G)

C---------open output files
open(40,FILE='out1',FORM='UNFORMATTED',STATUS='UNKNOWN'
& ,ACCESS='DIRECT',ACTION='WRITE',RECL=irec1)
open(50,FILE='out2',FORM='FORMATTED')
C----------- get LUT data
if(ABS(latit(x,y)-(knt-1)*DIF_LATG+90.0)
& .le. DIF_LATG*0.5) then
ilat=knt

exit
endif
enddo

You open two files for each run of this loop, superceding the ones that were opened in the run before without ever using them.

You never assign a value to x and y but use them as indices. So you allways reference latit (0,0) which produces an error.

Norbert



The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.

 

[Thardy]

Thanks for your reaction to my problem.
Even after assigning x and y respectively equal to 1236 and 1656, the issue is not solved.
An example of program structure different from mine would be appreciated.
Thardy

 

[salgerman]

Well, here is what I think you wanted to do...but it is up to you to decide whether I assumed the correct thing or not...

The thing is that you had a single loop for knt and nothing for x,y...so, I speculated that you meant to traverse the entire latit and longit matrices and evaluate every entry versus all your conditionals...so, I inserted a double loop for (x,y) inside the knt loop.

A run of this programs writes out 0.35298297

program LUT
      implicit none 
      integer*4:: lat, lon,x,y 
      INTEGER, PARAMETER :: NBYTE_G=4 !byte of albedo data
      INTEGER, PARAMETER :: NLONG_G=2880, NLAT_G=1441 !size of albed data
      REAL, PARAMETER :: DIF_LONGG=0.125, DIF_LATG=0.125 !size of each pixel of albedo data
      REAL*4 :: gal
      REAL(KIND=NBYTE_G), ALLOCATABLE :: galraw(:, :)
      REAL(KIND=4)::latit(1236,1656),longit(1236,1656)
      REAL(KIND=4)::longdi(1236,1656)
      INTEGER :: iost, knt, knt2, i, j
      INTEGER :: ilat, ilong
      INTEGER :: inlon, inlat
      integer, parameter ::irec1=1236*1656*4

      open(10,FILE='LON',FORM='UNFORMATTED',ACCESS='DIRECT',RECL=irec1)
      read(10, rec=1) longit
      close(10)  
      write(*,*)
      write(*,'("longit, sample elements:")')
      write(*,'("First 5: ",5F10.4)') (longit(1,j),j=1,5)
      write(*,'(" Last 5: ",5F10.4)') (longit(1236,j),j=1652,1656)

      open(20,FILE='LAT',FORM='UNFORMATTED',ACCESS='DIRECT',RECL=irec1) 
      read(20, rec=1) latit
      close(20)      
      write(*,*)
      write(*,'("latit, sample elements:")')
      write(*,'("First 5: ",5F10.4)') (latit(1,j),j=1,5)
      write(*,'(" Last 5: ",5F10.4)') (latit(1236,j),j=1652,1656)
      
      ! get LUT data
      allocate(galraw(NLAT_G, NLONG_G))
      open(30,FILE='LUT',FORM='UNFORMATTED',STATUS='OLD', &
                         ACCESS='DIRECT',ACTION='READ',RECL=NLONG_G*4)    
      do knt=1, NLAT_G
          read(30,REC=knt) (galraw(knt,knt2), knt2=1, NLONG_G)
! sample LUT data
!          if ( 100<knt .and. knt<105 ) then
!            write(*,*) (galraw(knt,j),j=1000,1005)
!          end if
!          if ( 1000<knt .and. knt<1005 ) then
!            write(*,*) (galraw(knt,j),j=NLONG_G-5,NLONG_G)
!          end if
      enddo
      100 continue
      close(30)
      
      write(*,*) "NLAT_G = ", NLAT_G
      write(*,*) "NLONG_G = ", NLONG_G
          
      ! find ilat index
      do knt=1, NLAT_G
        if (mod(knt,100)==0) write(*,*) "LAT, knt = ", knt, " of 1441" 
        do x=1,1236
        do y=1,1656
          if( ABS(latit(x,y)-(knt-1)*DIF_LATG+90.0) .le. DIF_LATG*0.5 ) then
              ilat=knt
              write(*,*) "Found ilat = ", ilat
              go to 200
          endif
        end do
        end do
      enddo
      200 continue
      
      ! find ilong index
      do knt=1, NLONG_G
        if (mod(knt,100)==0) write(*,*) "LON, knt = ", knt, " of 2880" 
        do x=1,1236
        do y=1,1656
      
        longdi(x,y)=(knt-1)*DIF_LONGG
        
        if(longit(x,y) .gt. 0) then
            if(longdi(x,y) .gt. 180.0) then
                longdi(x,y)=longdi(x,y)-360.0
            endif
        else if(longit(x,y) .lt. 0.0) then
            if(longdi(x,y) .ge. 180.0) then
                longdi(x,y)=longdi(x,y)-360.0
            endif
        endif
        
        if(ABS(longit(x,y)-longdi(x,y)) .le. DIF_LONGG*0.5) then
            ilong=knt 
            write(*,*) "Found ilong = ", ilong
            go to 300
        endif
        
        end do
        end do
      end do
      300 continue      
            
      gal=galraw(ilat, ilong)
      write(*,*) gal
      
      deallocate(galraw)
      stop
End program LUT

 

[Thardy]

Salgerman,
Thank you very much for re-arranging this program.I appreciate that. It works well, but I'm still facing the same issue as initially i.e. it processes only one value for latitude and longitude; that produces one value from the lookup table (0.35298297).
An array of 1236*1656 (same size as LON or LAT files) values are expected as result, rather than one value. May be I should explain better.
1)-The lookup table(LUT) is a global map file(array of 2280*1441). There are two sets of values in this file are ratios i.e. in the range 0.0 to 1.00, and -999.00 (replace missing data). All these values are along the longitude direction (x direction) and latitude direction (y direction) and each value is separated from the next value by 0.125 degrees (DIF_LONGG or 360/2880 and DIF_LATG or (90+90)/1441) either in the -x(supposed to cover a total range 0 to 360 degrees) or the -y(supposed to cover a total range of 90 to -90 degrees) direction.

2)-The LON(longitude) and LAT(latitude) are values in degrees (geographical location). These are first read. each LON value is combined to each LAT value at similar position in the arrays, to search for the nearest ratio value in the LUT coreeponding to the geographical location.
Let me mention that the LUT has only values from 0.0 to 1.00 and -999.00, these are not geographical values but are aligned to match global geographical coordinates i.e. from 0 to 360degres(left to right or x direction) and 90 to -90 degrees(up to down or y direction).
I will appreciate any help.
Thardy

 

[Thardy]

Hope to have expressed in detail the issue in the program I sent.
Anybody's help to re-arrange this program or propose a new program that will solve the lookup table problem will be appreciated.
Thanks
Thardy

 

[salgerman]

Thardy....don't give up on yourself. Also, we are not programmers for hire around here...we are here to assist you but you are the one to do the work...most of it, anyway.

Use what I indicated and figure out what else you need to accomplish what you need to accomplish.

Maybe you need to take a step back and with pencil and paper draw bubbles and connect them with arrows, describe in the bubbles what you need to do when, etc, ...this is your algorithm...THEN code it.

Here is more help...but I am not going to code for you, anymore.

I don't quite understand what you are trying to do but you may need a big outer loop for finding you next ilat,ilong pairs and retrieve the corresponding galraw(). Right now, it is being done only once.

Inside this outer loop is where you put the loops above to find ilat and ilong.

I think the trick to not keep finding the same ilat and ilong is to make it so that the starting value for the loop variables should be one more than the last index used.

So, instead of always starting from 1, you should set ilat and ilong to zero outside the outer big loop and then set your loops to start the first one from ilat+1 and the other from ilong+1.

I think that is what you need, but you know better.

good luck

 

[gummibaer]

Drawing bubbles and arrows is one option on how to proceed in your case, might well be what you need to do.

The other is - and you must do this yourself as you are the only living expert on what you expect your program to do - hunting down bugs. It's the same old story over and over again.
> locate the statement that should build the array you want
> find out why it does not do its job, is it executed only once in the run of your prog or is it simply that the index never changes
> find out what is the source of this problem, the parameter, tha variable that is not what it should be
> find out why it is not what you expect it to be
> correct this
> ...
> ...

You will have to decide, if your code does not work at all and you should redesign your algorithm, or if your algorithm seems fine but seems influenced by coding bugs.

This might be a lengthy process but it is the only one to help and you are the only one to do the job. Redesigning your algorithm and coding it , will sure lead to other bugs, so you will find yourself debugging your code anyway nad doing this will practically compute the first set of data manually parallel to what your computer does.

Or find somebody that does it for you on a payment basis.

And yes, if you have problems with Fortran, ('What would be the best way to ... ', 'Why does ... not work ', 'I do not understand how ....')
feel free to come back and ask for help, which you sure will find here. But do not expect us to optimise your algorithm and worm our ways in latitudes and longitudes and geodetical data.

$ 0.02

Norbert

The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.