[Gdal-dev] IPCC raster format

[Ivan Lucena]

Hi Tim,

I was about to write to Frank Warmerdam and Christopher Baker to say thanks for their suggestions.

Mine implementation works like that:

$ ipcc2gdal -of <writable driver> <input file> <output file "prefix">

Example:

$ ipcc2gdal -of gtiff cpre0130.dat cpre0130.tif
cpre0130_1.tif
cpre0130_2.tif
cpre0130_3.tif
cpre0130_4.tif
cpre0130_5.tif
cpre0130_6.tif
cpre0130_7.tif
cpre0130_8.tif
cpre0130_9.tif
cpre0130_10.tif
cpre0130_11.tif
cpre0130_12.tif

So thanks everybody,

Ivan 

________________________________________
From: Tim Sutton [mailto:timlinuxgis at gmail.com] On Behalf Of Tim Sutton
Sent: 18 June 2006 19:27
To: Ivan Lucena; Peter Brewer
Cc: Tim Sutton; gdal-dev at lists.maptools.org
Subject: Re: [Gdal-dev] IPCC raster format

Hi

I wrote a tool in c++ once to read this and if I recall correctly spit out arc/info ascii grid files. Its been a while since I looked at it so Ill need to dig it up and recheck exactly what it did. Pete can you remember if this is in the climate data processing wizard? I have no network access as I write this so I will only be able to confirm for you tomorrow what I actually have. Ivan if you still  have not found a suitable solution and Ill dig out what I have.

Regards

Tim

On 16/06/2006, at 12:08, Ivan Lucena wrote:


Hi there,
 
I would rather been watching the World Cup, but I have to get this data in a usable format and there is no better place to get it than from GDAL's experts.
 
The "Intergovernmental Panel on Climate Change" (IPCC) offers data on the web in a format described by this Fortran source code: (The file extension is just a generic ".dat")
 
http://ipcc-ddc.cru.uea.ac.uk/obs/observed_fileformat.html
 
> more cpre0130.dat
    grd_sz      xmin      ymin      xmax      ymax    n_cols    n_rows  n_months
      0.50      0.25    -89.75    359.75     89.75       720       360        12
...
-9999-9999-9999-9999-9999-9999-9999-9999-9999   52   88  103   91   65   45   27   32 ...
-9999-9999-9999-9999-9999   11   20   22   25   27   26   27   28   31   29   28   30 ...
...
 
It's not hard to analyze the pattern on the data or the code, so I figured that:
 
- It's an ASCII format;
- There are two lines of reader (easy to read)
- Each cell is 5 characters wide;
- The no-data is "-9999";
- There is no separator between cells (what makes my task difficult);
- It's a multi-band raster dataset with 12 bands; (720x360x12 cells)
- I am not sure if positive values can go beyond 9999 (that reminds me COBOL)
 
I am heading in the direction of writing a Python+GDAL script to do the task, but I am stuck in how to separate the cells. The Fortran code is doing it in just one command as you can see in the code. I am thing that maybe by using the NumPY stuff embedded in GDAL I could do the same.
 
Any Idea? Any experience with this data format?
 
Thanks in advance,
 
Ivan
 
-------------------------------------------------------
This is a Python command line interaction with the data:
 
>>> import gdal
>>> f = open('C:/Data/IPCC/cpre0130.dat', 'r');
>>> f.readline()
'    grd_sz      xmin      ymin      xmax      ymax    n_cols    n_rows  n_months\n'
>>> f.readline()
'      0.50      0.25    -89.75    359.75     89.75       720       360        12\n'
>>> f.readline()
'-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9
 
-------------------------------------------------------
Fortran source code:
 
program rd_ascii       
!
! f90 program to read in an integer ascii grid with
! variable dimensions into a global grid (720x360x12)
!
  parameter :: n_cols=720
  parameter :: n_rows=360
  integer, dimension (n_cols,n_rows,12) :: grid
  character(len=72) :: infl
  character(len=20) :: fmt
!
  call getarg(1,infl)
!       
  if(infl.eq.' ')then
    write(*,*) 'Enter ascii grid file name'
    read(*,'(a72)')infl
  end if
!       
  open(1,file=infl,status='old')
  read(1,*)
  read(1,*)xmin,ymin,xmax,ymax,ncols,nrows,nmonths,missing
  grid=missing
  fmt='( i5)'
  write(fmt(2:4),'(i3)')n_cols
  do im=1,nmonths
    do lat=1,n_rows
      read(1,fmt)(grid(lon,lat,im),lon=1,n_cols)
    enddo
  enddo
!       
end program rd_ascii
 
 
 
_______________________________________________
Gdal-dev mailing list
Gdal-dev at lists.maptools.org
http://lists.maptools.org/mailman/listinfo/gdal-dev

Tim Sutton
tim at linfiniti.com