R код для чтения данных в формате ncdf

Feb 10, 2014 15:35

1. установить и загрузить библиотеку (ncdf)
2. скачать ncdf файлы по ссылкам на базы данных. Эти ссылки можно найти в моем верхнем посте. В этом примере использованы данные SST (Sea Surface Temperature, http://www.esrl.noaa.gov/psd/data/gridded/data.noaa.oisst.v2.highres.html), несколько этих файлов можно хранить в одной папке, код найдет их и откроет годы один за другим, если нужна выборка за несколько лет.
3. вручную задать год начала и конца периода наблюдений.
4. вручную задать координаты интересующей точки.

В этом коде ввод координат и времени происходит вручную, но в моих рабочих кодах обычно эта информация вытаскивается автоматически из исходных биологических данных.

#########################################################
### Reading climatological data from .ncdf files
### In this example it is Sea Surface temperature (SST).
### To read other files with other parameters simply add pathways to those data in "open.ncdf" function.

library(ncdf) # needed to read this kind of files

### Get data from ncdf files automatically

### INPUT first and last year
year_start <- 2006 # enter the start year manually

year_end <- 2013 # enter the end year manually

list_years<- seq(year_start, year_end, by=1)

### Here years could be selected manually
#list_years<- seq(2005, 2006, by=1)

list_years
l<-length(list_years)
l
#list_years[1]
### end of Find out the time range of data and read

### INPUT one geographic location of the station (or any point of interest), coordinates are in decimal degrees
# code will later select the closest point of the grid to the station (or to the point of interest)
lat<- 56 # input latitude manually, here it is Mejlflak
lon<- 10.5 # input longitude manually

### Prepare the empty dataframe to collect annual data into single multi-annual dataset
data_SST_total = data.frame(matrix(vector(), 0, 4, dimnames=list(c(), c("date_num", "date_date", "date_date_posix", "node_SST"))), stringsAsFactors=F)
str(data_SST_total)

### loop to read large ncdf data, SST, from many years
for(i in 1:l) {
year_global.nc = open.ncdf(paste("X:/24_konstante Daten/hydrographic_datasets/SST_daily_mean_global_grid_0_25x0_25degree/sst.day.mean.",list_years[i],".v2.nc", sep=""))
# this line opens the connection to large file, it does not read it into the memory, so it does not take time
# this folder contains global 1/4 degree grid with daily SST data for the period between 1980 till 2013
# year_global.nc

### create a date format for correct representation of time information from ncdf
date_num<-year_global.nc$dim$time$vals
str(date_num)
date_num_vector<-as.vector(date_num)
date_num_vector
data_SST <- data.frame(cbind(date_num=date_num)) #, dir_list_degrees_factor=dir_list_degrees, dir_list_deg=0, proc=0))
data_SST
data_SST$date_date<-as.Date(data_SST$date_num, origin="1800-01-01") # origin="1800-01-01" is specific for this ncdf files
data_SST$date_date_posix<-as.POSIXct(paste(data_SST$date_date, "00:00:00", "%H:%M:%S"))
str(data_SST)
head(data_SST)
#date_SST$date_date<-as.POSIXct(paste("1899-12-30", format(date_SST$date_num, "%dd:%dd:%yyyy")))

dist_lon = abs(year_global.nc$dim$lon$vals-lon)
dist_lon
index_lon = which.min(dist_lon)
index_lon
#return( index )
dist_lat = abs(year_global.nc$dim$lat$vals-lat)
dist_lat
index_lat = which.min(dist_lat)
index_lat
node_SST<-get.var.ncdf(year_global.nc, start=c(index_lon, index_lat, 1), count=c(1,1,length(year_global.nc$dim$time$vals))) # start= c() tells from where to start reading data, count=c() tells how many cells to read along each axes. Both vectors should be in format c(x,y,z,t)
#year_global.nc$dim$time$vals
#length(year_global.nc$dim$time$vals)
data_SST$node_SST <- node_SST
str(data_SST)
data_SST_total<-rbind(data_SST_total, data_SST)
#data_SST_total

close.ncdf(year_global.nc) # closes connection to currect ncdf file
}
str(data_SST_total)
data_SST_total

### code to plot SST, inflag to take a quick look on data

plot(data_SST_total$date_date, data_SST_total$node_SST, xlab="time", ylab="SST, degrees C")

### end of the code to plot SST

### write table with data to the file
write.csv(data_SST_total, file = paste(getwd(),"/", year_start, "_", year_end, "_lat", lat, "_lon", lon, "_SST.csv", sep=""))

### End of Get data from ncdf files automatically
###################################################

программирование, климат, программы, работа, r, базы данных

Previous post Next post
Up