Skip to content

reading into R IPCC data files

10 messages · Mario Bellinato, Thomas Adams, Barry Rowlingson +2 more

Original
Thomas Adams
# 
Marco,

I just took a look at the data; unfortunately, it's not in a very usable format without reformatting it for
importing the data into either R or a GIS, such as GRASS GIS. I'm sure an experienced R programmer could
write a clever R script to do what you want. Personally, I would write a Perl script to reformat the data to
import directly into GRASS GIS. I looked at the HADCM3_A2a_PREC_1980.tar.gz data and it is cell-centered.
GRASS grids are referenced by the lower-left grid corner. I noticed that there are many different grid spacings
and not all the grids are square, they are rectangular. Bottom line is that a lot of care is going to be needed
to manipulate these data sources. I have to think that somewhere there must be this same data available in grib
format, which imports into GRASS very easily.

Regards,
Tom
On 1/15/11 6:54 PM, Mario Bellinato wrote:

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Barry Rowlingson
    

    

      
      #
    

  


  

      
On Sun, Jan 16, 2011 at 3:28 PM, Thomas Adams <Thomas.Adams at noaa.gov> wrote:

            

      
      
      
    

        
That's a horrendous data file format, typical of the sort of thing you
get from Fortran programmers :)

There's a very complex-looking "web 0.1" style data browser thing here
which lets you eventually download .csv files which might be easier to
read into R. I think you can only get a month at a time though, and it
might not be the same data as in your original link (but the Hadley
centre stuff seems to be there...):

http://www.ipcc-data.org/cgi-bin/ddc_nav/dataset=ar4_gcm

 You'd almost think they don't want us to know about climate change....

Barry

  
  


  


      

    

    

      
      

        


  

        

      Thomas Adams
    

    

      
      #
    

  


  

      
Barry,

WRT your last line ? I could not agree more; the same thought came to me 
as well!

Tom

        
On 1/16/11 11:04 AM, Barry Rowlingson wrote:

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Barry Rowlingson
    

    

      
      #
    

  


  

      
On Sun, Jan 16, 2011 at 5:35 PM, Thomas Adams <Thomas.Adams at noaa.gov> wrote:

            

      
      
      
    

        
The counter-argument of course is that these guys are hardcore
programmers who don't have time to put a wussy presentation layer on
top of the data! *Their* Fortran code tells them climate change is
real, why should they have to convert everything to ESRI Raster Grids
when it's so obvious? :)

Barry

  
  


  


      

    

    
1 day later

    

      
      

        


  

        

      Andy Bunn
    

    

      
      #
    

  


  

          

      
      
      
    

        
Two things:
1. Casual perusal of the AR4 data can be done via the Visualisation Tools. I think it's well done overall in terms of accessing and plotting a ludicrous amount of GCM data. Look here: 
http://www.ipcc-data.org/ddc_visualisation.html 

2. If you want the actual zeros and ones, model output from the AR4 is available from the Coupled Model Intercomparison Project (CMIP3 - specifically). Unlike the IPCC DDC link given originally in this thread, the CMIP3 data are easier to parse for non F90 folks and available as netCDF. Look here: http://www-pcmdi.llnl.gov/ipcc/about_ipcc.php 

There is a nice overview of the model simulation archive here:
http://www.realclimate.org/index.php/archives/2008/02/ipcc-archive/ 

HTH, -AB

  
  


  


      

    

    

      
      

        


  

        

      Robert J. Hijmans
    

    

      
      #
    

  


  

      
I agree with Andy, although dealing with ncdf files has its own problems
((although you can open most with raster: brick(filename))

So, for the record, reading these ASCII files is not that hard if you have
R. See below for an example that works on the one file I looked at. You
could easily extend this to automate downloading and processing all the
files with functions download.file() and uzip() 

Robert


library(raster)

filename <- HADCM3_A1F_DSWF_1980.mea

x <- readLines(filename)
x <- apply(matrix(x),1, function(x) gsub("  ", " ", x))
x <- apply(matrix(x),1, function(x) trim (gsub("  ", " ", x)))

start = which(substr(x,1,5)=="7008 ")+1
end = c(which(substr(x,1,4)=="IPCC")[-1] - 1, length(x))

nc = as.integer(substr(x[2], 8, 11))
nr = as.integer(substr(x[2], 13, 16))

# "Grid is 96 * 73  Month is Jan"         

if (length(start) != 12) {
	stop('something wrong')
} else {
	d = matrix(nrow=nc *nr, ncol=12)
	month = NULL
	for (m in 1:12) {
		d[,m] <- as.numeric(unlist(strsplit(x[start[m]:end[m]], " ")))
		mm <- x[start[m]-4]
		mm <- trim(substr(mm, nchar(mm)-3, nchar(mm)))
		month <- c(month, mm)
	}
}

b <- brick(nrow=nr, ncol=nc)
b <- setValues(b, d)
b <- rotate(b)
ln <- paste(month, ";", x[1], ";", x[3], ";", x[4])
layerNames(b) <- ln

plot(b, 1)
plot(b)

b <- writeRaster(b, 'filename.tif')

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Barry Rowlingson
    

    

      
      #
    

  


  

      
On Mon, Jan 17, 2011 at 11:04 PM, Robert Hijmans <r.hijmans at gmail.com> wrote:

        

            

      
      
      
    

        
The important point here is that you say it works on the _one_ file
you looked at. Probably works on others. You don't know until you've
tried. It might read it in without giving any errors but get all the
wrong numbers and then polar bears die.

 There's no excuse for publishing ASCII data sets like this when there
are open standards for data in this format, including all the metadata
hacked on to the top of the ASCII. This makes automated, repeatable
analysis hard to do.

 Probably preaching to the converted here anyway. I have done some
work with the environmental modelling community and they are in their
own little world sometimes with regard to code development (complete
lack of any concept of a revision control system) and standards usage
(ascii in, ascii out). Made me appreciate R even more!

Barry

  
  


  


      

    

    

      
      

        


  

        

      Mario Bellinato
    

    

      
      #
    

  


  

      
Thank you for all for your useful suggestions and the interesting
discussion! I am currently trying use Professor Hijmans's code to read
in the same file. I have already encountered one problem though. The
code runs fine until:

b <- rotate(b)

when I get the following error message:

Error in intersectExtent(x, y) : Invalid extent

This is the the RasterBrick object obtained from
b <- setValues(b, d)

b

class       : RasterBrick
filename    :
nlayers     : 12
nrow        : 73
ncol        : 96
ncell       : 7008
projection  : +proj=longlat +datum=WGS84
min value   :  0  0 10  0  0  0  0  0  6  0 ...
max value   : 427 336 324 352 360 396 375 355 323 354 ...
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
resolution  : 3.75, 2.465753  (x, y)


Any help would be greatly appreciated! I am running R 2.12.1 on Ubuntu 10.04

Best wishes,

Mario

  
  


  


      

    

    

      
      

        


  

        

      Robert J. Hijmans
    

    

      
      #
    

  


  

          

      
      
      
    

        
I think this goes away if you update the raster package from CRAN. 

I noticed one other thing that needs to be done. After

b <- brick(nrow=nr, ncol=nc) 

You need to do:

xmin(b) = 0
xmax(b) = 360

And then continue. 

Barry is right in saying that:

            

      
      
      
    

        
Although it should not be too hard to fix the script as you muddle along,
you obviously need to very careful. 

R