#####Sample R script to open "mann2008original.nc", from NOAA/NCDC Paleoclimatology World Data Center.  Author: Eugene R. Wahl, 9-8-2010.
#####This script can also be used with minimal modification to open any file of the form "mann2008XXXXXXX.nc".

#####NOTE:  The primary data grid variable is a site-by-time matrix.  It is similarly brought into the R workspace as either site-by-time
#####       or time-by-site matrices.  It is not in a more general array format with additional latitude and longitude coordinates.
#####       This was done because to add lat and lon coordinates would have greatly expanded the size of the array, by a factor of the 
#####       square of the number of sites (1209).

#####       Example code to use for sub-partitioning the data array according to lat and lon values, and also according to the 
#####       Mann et al. 2008 proxy data type code system, is provided at the end of the script.



#####This script was developed for R 2.8.1.  It likely is back-compatible with earlier versions of R, but this has not been systematically checked.



##source("PUT APPROPRIATE FILEPATH AND FILENAME HERE")  #USE THIS LINE TO SOURCE THIS FILE DIRECTLY IN R, ONCE FILEPATH AND FILENAME ARE SPECIFIED FOR LOCAL CIRCUMSTANCE.


# rm(list=ls()) #  Commented out, Remove pound sign to activate. ACTIVATION REMOVES ALL OBJECTS IN R WORKSPACE -- BE CAREFUL TO INSURE YOU WANT TO DO THIS.

library (ncdf)


##Set-up Information  


filepath<-("//Dinobast/paleo_fileshare/2000-reconstruction/Mann2008Proxy/")  #REPLACE WITH APPROPRIATE FILEPATH FOR LOCAL CIRCUMSTANCE.

filename<-("mann2008original.nc")


##Open netCDF File


nc<-open.ncdf(paste(filepath,filename,sep=""))
print (nc)

time<-get.var.ncdf(nc,"time")
dim(time)
time

site<-get.var.ncdf(nc,"site")
dim(site)
site

site_name<-get.var.ncdf(nc,"site_name")
dim(site_name)
site_name

network_info<-get.var.ncdf(nc,"proxy_network_info")  #Provides metadata of latitude, longitude, and proxy data type code for each site in one variable
dimnames(network_info)<-list(c("lat","lon","data_type"),site)
dim(network_info)
network_info

lat<-get.var.ncdf(nc,"lat")
dim(lat)
lat

lon<-get.var.ncdf(nc,"lon")
dim(lon)
lon

data_type<-get.var.ncdf(nc,"data_type")
dim(data_type)
data_type


proxy_data<-get.var.ncdf(nc,"proxy_data")
dimnames(proxy_data)<-list(site,time)   
dim(proxy_data)

proxy_data_trans<-t(proxy_data)
dim(proxy_data_trans)


sample_depth<-get.var.ncdf(nc,"sample_depth")
dimnames(sample_depth)<-list(site,time)
dim(sample_depth)

sample_depth_trans<-t(sample_depth)
dim(sample_depth_trans)


##rm(proxy_data)                   #Removes specified variable if activated:  used to conserve memory space if needed
##rm(sample_depth)



##View R array of sample of proxy data


proxy_data_trans[time<=2003 & time>=1955,site<=1209 & site>=1200]   



##An example way to sub-select proxies by latitude and longitude


lat_partition<-(network_info[1,])
lat_partition<-cbind(lat_partition,site)
lat2<-round(lat,digits=2)
dimnames(lat_partition)<-list(lat2,c("lat","site"))

lon_partition<-network_info[2,]
lon_partition<-cbind(lon_partition,site)
lon2<-round(lon,digits=2)
dimnames(lon_partition)<-list(lon2,c("lon","site"))


lat_partition2<-lat_partition[lat2>=-50 & lat2<=10,]     #Example values input into first dimension's logic statement to define lat partition

lon_partition2<-lon_partition[lon2>=-90 & lon2<=-30,]    #Example values input into first dimension's logic statement to define lon partition


lat_lon_partition<-match(lat_partition2[,2],lon_partition2[,2], nomatch = 0)
   for (i in 1:length(lat_lon_partition)) {
     if (lat_lon_partition[i]!=0) {
       lat_lon_partition2<-lat_lon_partition[i]
          break }}
   for (j in (i+1):length(lat_lon_partition)) {
     if (lat_lon_partition[j]!=0) {
       lat_lon_partition2<-cbind(lat_lon_partition2,lat_lon_partition[j]) }}
lat_lon_partition3<-lon_partition2[lat_lon_partition2,]


proxy_data_trans2<-proxy_data_trans[,lat_lon_partition3[,2]]   

proxy_data_trans2[1:20,1:40]



##An example way to sub-select proxies by proxy data type


data_type_partition<-(network_info[3,])
data_type_partition<-cbind(data_type_partition,site)
dimnames(data_type_partition)<-list(data_type,c("data_type","site"))


data_type_partition2<-data_type_partition[data_type==2000,]    #Example value input into first dimension's logic statement to define data type partition


proxy_data_trans3<-proxy_data_trans[,data_type_partition2[,2]]

proxy_data_trans3[1:10,1:40]



#####END