MATTS Tutorials netCDF R LOAD RNCDF #a becomes name of data a=open.nc("D:\\BEP02.cdf") #print ncdf information print.nc(a) #get variable tke var.get.nc(a,'tke') #plot tke b=var.get.nc(a,'tke') plot(var.get.nc(a,'tke')) #get first 100 values of tke var.get.nc(a,'tke',count=100) #get 100 values of tke from 99 var.get.nc(a,'tke',start=99,count=100) ReadFile R FOR FILE WITH HEADER and DATA: dectime MODEL11 MODEL12 MODEL13 MODEL14 MODEL16 MODEL17 MODEL21 MODEL22 MODEL25 MODEL28 MODEL30 MODEL31 MODEL32 MODEL33 MODEL35 MODEL36 MODEL37 MODEL39 MODEL42 MODEL43 MODEL45 MODEL46 MODEL47 MODEL49 MODEL50 ALPHA Quality TimePeriod Quality 0 = GOOD FOR ALL FLUXES day = 2 night = 1 transition = 0 ======================PROGRAM=================================== #read in data from directory dataSet<-read.csv(file="V:\\2009\\ALPHA OUTPUT\\Stage1_QE.csv",head=TRUE,sep=",") #print decimal time - by file and column name print(dataSet$dectime) #print model 16 data print(dataSet$MODEL16) #model 16 data = mod16 mod16=dataSet$MODEL16 #decimal time data = decTime decTime=dataSet$dectime #plot mod16 plot(mod16) #plot mod16 against decTime - careful here! plot(decTime,mod16) #only good ALPHA data - i.e. array of line numbers that meet this condition goodData=which(dataSet$Quality == 0) #plot only good values plot(dataSet$ALPHA[goodData]) #plot between fixed time period ie. 100-300 plot(dataSet$MODEL12[100:300]) #mean of good data for model 13 print(mean(dataSet$MODEL13[goodData])) #mean of day data dayTime=which(dataSet$TimePeriod == 2) print(mean(dataSet$MODEL13[dayTime])) Ensemble MEAN Derive ENSAMBLE MEAN for each month With 25 models of data, we want to end up with 25 plots per month Cut data to 225/year 1 – 333/year 2 Need to consider number of days in month and how hours fit in: e.g. Jan = 31 days = 744 hrs = 1488 observations. Need to extract these observations array[x:x+1488] Break into days of 48 values, eg. Day1 Put into matrix Jan=matrix(NA,31,48) Matrix of NA values of 31 rows and 48 columns To put day 1 info into matrix: Jan[1,]=day1 First row becomes first day Do for all days using for statement: F (in 1:31) { Jan[i,]=dayi } Get mean for each half hour using: arrayMeanJan=mean(Jan[,1...48]) Plotting Get 5 by 5 plots for each month par(mfcol=c(5,5), mar=c(1.5,1.5,1,1.1), oma=c(5,5,5,2.5)) Then for each model: Plot(hours,arrayMeanJan) lines(hours,meanObs)