[R] problem with loop to put data into array with missing data forsome files
Jenny Barnes
jmb at mssl.ucl.ac.uk
Mon Nov 20 15:39:02 CET 2006
Dear R-help,
Thank you to all who answered my question on how to solve my problem (original
message is shown below). I think I have solved it now with your help. For those
who are interested:
I inserted another if loop within the code at the line before I called upon one
of my assigned variables (ny) that would have caused a problem as it would not
have been there as the record did not exist:
#This bit was in the code before:
param <- status[2]
dims <- unlist(strsplit(status[3], " "))
nx <- as.numeric(dims[13])
ny <- as.numeric(dims[15])
test <- dims[15]
#Here's the new bit:
if(is.na(test)){
data.out$date[k] <- paste(i, j, mon2, sep="")
NAindata <- matrix(NA, 94*192)
data.out$data[k,] <- NAindata
}else{
#then I continued with the code as before but with an extra } at the end
to close off the if loop.
Thanks once again for everyone's time and email inbox space,
Jenny
>>-----Original Message-----
>>From: r-help-bounces at stat.math.ethz.ch
>>[mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Jenny Barnes
>>Sent: 20 November 2006 13:03
>>To: r-help at stat.math.ethz.ch
>>Subject: [R] problem with loop to put data into array with missing data
>>forsome files
>>
>>
>>Dear R-help community,
>>
>>My main goal of this message is to find a way of skipping a file of a
>>month/year
>>in a loop that does not exist (and making it's output into an data.out
>>array
>>would be NA) and moving onto the next year/month in the loop to carry
>on
>>filling
>>data.out with real precipitation data.
>>
>>The situation so far:
>>I downloaded 50 years worth of GRIB data files from the NCEP data site
>>http://nomad3.ncep.noaa.gov/pub/reanalysis-1/month/grb2d.gau/
>>
>>I then created a loop in R to read each month of each of the 50 years
>>worth of
>>files and only extract the precipitation records using wgrib and grep
>as
>>show in
>>the code at the end of this message. I had to use grep to extract the
>>PRATE
>>records each time, as unfortunately with these grib files precipitation
>>is not
>>always the same record number each time.
>>
>>This worked fine until Feb 1972 (which is file
>>/home/jmb/sst_precip/gribdata/197202.grib). This file did not download
>>completely and as a result the file that saved to my computer did not
>>include
>>the record for precipitation - only the first 13 records. At this point
>>in the
>>loop it just stopped with error message:
>>Error in seq.default(lat.beg, lat.end, length = ny) :
>> length must be non-negative number
>>Which is obviously because ny does not actually exist for this file as
>>it cannot
>>be extracted from source as this record does not exist for this file!
>>
>>My question is:
>>Does anybody have any suggestions as to how I could code it so that my
>>loop will
>>allow this file to just output NA into the data.out array in
>>data.out$data and
>>maybe keep FALSE in data.out$date for the month&years where this has
>>occured and
>>then continue to loop around all the other months and years in
>>succession?
>>
>>Things that don't work that I have tried:
>>1)I cannot just restart the loop after any anomalous years as this
>>"resets"
>>data.out and changes its dimensions depending on the number of years
>>looped
>>(nyrs) etc.
>>2)The first point of trouble is when I call upon ny in the line
>>>lats <- seq(lat.beg, lat.end, length=ny)
>>as ny is NA in this record - so I tried putting the line
>>>test <- dims[15]
>>>if(test != "NA"){
>>with the closing } at the end of the rest of the loop
>>...this gives the error message
>>>Error in if (test != "NA") { : missing value where TRUE/FALSE needed
>>3) My thought was that maybe I would need something like
>>>if (parameter == "PRATE" ){
>>then continue with loop; but this will obviously not work by itself and
>>gives
>>the error message
>>>Error in if (test != "NA") { : missing value where TRUE/FALSE needed
>>- would need to specify what to do if parameter is not =="PRATE". I
>>encourage
>>anybody with another idea to contact me with their suggestions!
>>
>>Thank you for taking the time to read my problem,
>>
>>Sincerely,
>>
>>Jennifer Barnes
>>
>>r-help at stat.math.ethz.ch
>>Here is my code:
>>
>>#takes precip data for every year and month needed
>>#from file saved on computer already
>>#and sorts it into a list (array)
>>
>>#set some control parameters
>>
>>startyr <-1954 #note needs to be year before data wanted
>>nyrs <- 50
>>months <- nyrs*12
>>k <- 0
>>library(chron)
>>library(utils)
>>data.out <- list(lats=seq(88.542, -88.542, length=94),
>> lons=seq(0, 360-1.875, length=192),
>> date=vector(length=months),
>> data=array(NA, c(months, 94*192))
>>)
>>#the loop
>>for(yr in 1:nyrs){
>> for(mon in 1:12){
>> i<- startyr+yr
>>
>> if(mon < 10) {j <- 0}
>> if(mon < 10) {mon2 <- mon}
>> if(mon == 10) {(j <- 1)}
>> if(mon == 10) {mon2 <- 0}
>> if(mon == 11) {(j <- 1)}
>> if(mon == 11) {mon2 <- 1}
>> if(mon == 12) {(j <- 1)}
>> if(mon == 12) {mon2 <- 2}
>>
>> outfile <- paste("/home/jmb/sst_precip/gribdata/",i, j, mon2,
>>".grib",
>>sep="")
>> #decode the grib
>> cmd <- paste("/usr/local/bin/wgrib ", outfile, "| grep
>>\":PRATE:\" |
>>/usr/local/bin/wgrib", outfile, "-i -V -text -o /tmp/wgrib.output")
>> status <- system(cmd, intern=T)
>>
>> #extract metadata needed for working with julian dates:
>> dt <- unlist(strsplit(status[1], " "))[3]
>> dt.yr <- as.numeric(substr(dt, 1,4))
>> dt.mon <- as.numeric(substr(dt, 5,6))
>> dt.day <- as.numeric(substr(dt, 7,8))
>> dt.hr <- as.numeric(substr(dt, 9,10))
>> dt.jul <- julian(dt.mon, dt.day+(dt.hr/24), dt.yr)
>>
>> parameter <- unlist(strsplit(status[1], " "))[4]
>>
>> param <- status[2]
>> dims <- unlist(strsplit(status[3], " "))
>> nx <- as.numeric(dims[13])
>> ny <- as.numeric(dims[15])
>> temp <- unlist(strsplit(status[5]," "))
>> lat.beg <- as.numeric(temp[6])
>> lat.end <- as.numeric(temp[8])
>> lats <- seq(lat.beg, lat.end, length=ny)
>> temp <- unlist(strsplit(status[6]," "))
>> lon.beg <- as.numeric(temp[14])
>> lon.end <- as.numeric(temp[16])
>> if(lon.beg > lon.end){lon.beg <- (lon.beg -360)} #since sp
>>prefers -180
>>to 180 encoding
>> lons <- seq(lon.beg, lon.end, length=nx)
>> lons[lons > 180] <- lons[lons > 180]-360 #since sp
>>prefers -180
>>to 180 encoding
>>
>> #read flat file into R (header first)
>> dims <- scan("/tmp/wgrib.output", nlines=1, quiet=T)
>> nrec <- dims[1]*dims[2]
>> indata <- paste("indata.", i, j, mon2, sep="")
>> indata <- scan("/tmp/wgrib.output", 'numeric', nlines=nrec,
>>skip=1,
>>quiet=T)
>> indata <- matrix(as.numeric(indata), nrow=dims[1], ncol=dims[2],
>>
>>byrow=F)
>>
>> #to get precip into mm/month, not mm/sec:
>> if((dt.mon==02)&(dt.yrendnumber%%4==0)){days.mon <- 29} else
>>{days.mon
>><-28} #to take into account leap-years!!!
>> if(dt.mon==01){days.mon <- 31}
>> if(dt.mon==03){days.mon <- 31}
>> if(dt.mon==05){days.mon <- 31}
>> if(dt.mon==07){days.mon <- 31}
>> if(dt.mon==08){days.mon <- 31}
>> if(dt.mon==10){days.mon <- 31}
>> if(dt.mon==12){days.mon <- 31}
>> if(dt.mon==04){days.mon <- 30}
>> if(dt.mon==06){days.mon <- 30}
>> if(dt.mon==09){days.mon <- 30}
>> if(dt.mon==11){days.mon <- 30}
>> sec.mon <- days.mon*24*60*60
>>
>> newindata <- paste("newindata.",i, j, mon2, sep="")
>> newindata <- sec.mon*indata
>>
>> k<-k+1
>> data.out$date[k] <- paste(i, j, mon2, sep="")
>> data.out$data[k,] <- newindata
>>}
>>}
>>
>>
>>~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
>>Jennifer Barnes
>>PhD student - long range drought prediction
>>Climate Extremes
>>Department of Space and Climate Physics
>>University College London
>>Holmbury St Mary, Dorking
>>Surrey
>>RH5 6NT
>>Web: http://climate.mssl.ucl.ac.uk
>>
>>______________________________________________
>>R-help at stat.math.ethz.ch mailing list
>>https://stat.ethz.ch/mailman/listinfo/r-help
>>PLEASE do read the posting guide
>>http://www.R-project.org/posting-guide.html
>>and provide commented, minimal, self-contained, reproducible code.
>
>~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
>Jennifer Barnes
>PhD student - long range drought prediction
>Climate Extremes
>Department of Space and Climate Physics
>University College London
>Holmbury St Mary, Dorking
>Surrey
>RH5 6NT
>01483 204149
>07916 139187
>Web: http://climate.mssl.ucl.ac.uk
>
>
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Jennifer Barnes
PhD student - long range drought prediction
Climate Extremes
Department of Space and Climate Physics
University College London
Holmbury St Mary, Dorking
Surrey
RH5 6NT
01483 204149
07916 139187
Web: http://climate.mssl.ucl.ac.uk
More information about the R-help
mailing list