Last active
October 26, 2016 11:54
-
-
Save aolinto/7d3f282c9fde96133daa4cdd4f1bbcab to your computer and use it in GitHub Desktop.
R script to extract Aqua MODIS Chlorophyll a data from a nc file
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# Antonio Olinto Avila-da-Silva, Instituto de Pesca, Brasil | |
# ver 2016-05-17 | |
# https://gist.github.com/aolinto/7d3f282c9fde96133daa4cdd4f1bbcab | |
# script to process Aqua MODIS Sea Surface Temperature | |
# files downloaded from http://oceancolor.gsfc.nasa.gov/cgi/l3 | |
# Aqua MODIS Clorophyll Concentration, OCI Algorithm Monthly 9 km SMI images | |
# all .L3m_MO_CHL_chlor_a_9km.nc files must be in the working directory | |
# the script will open each nc file, read date, lon, lat and chl data, | |
# then select data from specified area and write them into | |
# a single csv file named MODISA_chl.csv | |
# Some reference pages | |
# http://geog.uoregon.edu/GeogR/topics/netCDF-read-ncdf4.html | |
# https://scottishsnow.wordpress.com/2014/08/24/many-rastered-beast/ | |
# load libraries | |
# ncdf4 needs libnetcdf-dev netcdf-bin in Linux | |
# install.packages(c("ncdf4","reshape2")) | |
library("ncdf4") | |
library("reshape2") | |
# set working directory | |
setwd("/mnt/Dados02/MODISA/CHL_NC") # indicate the path to the files | |
file.exists("MODISA_chl.csv") # caution new data will be appended to this file if it already exists | |
# file.rename("MODISA_chl.csv","MODISA_chl.old") | |
# file.remove("MODISA_chl.csv") | |
# list and remove objects | |
ls() | |
rm(list = ls()) | |
# set the study area in decimal degrees | |
lonmax<- -45.92 | |
lonmin<- -46.45 | |
latmax<- -23.80 | |
latmin<- -24.22 | |
# create a list of files and indicate its length | |
(f <- list.files(".", pattern="*.L3m_MO_CHL_chlor_a_9km.nc",full.names=F)) | |
(lf<-length(f)) | |
# variable | |
var<-"chlor_a" | |
for (i in 1:lf) { | |
# progress indicator | |
print(paste("Processing file",i,"from",length(f),f[i],sep=" ")) | |
# open netCDF file | |
data<-nc_open(f[i]) | |
# extract data | |
lon<-ncvar_get(data,"lon") | |
lat<-ncvar_get(data,"lat") | |
value<-ncvar_get(data,var) | |
unit<-ncatt_get(data,var,"units")$value | |
# matrix to data.frame | |
dimnames(value)<-list(lon=lon,lat=lat) | |
dat.var<-melt(value,id="lon") | |
# select data from the study area taking out missing data | |
# dat.varSAtmp<-na.omit(subset(dat.var,lon<=lonmax & lon>=lonmin & lat<=latmax & lat>=latmin)) | |
dat.varSAtmp<-subset(dat.var,lon<=lonmax & lon>=lonmin & lat<=latmax & lat>=latmin) | |
# extract date information | |
dateini<-ncatt_get(data,0,"time_coverage_start")$value | |
dateend<-ncatt_get(data,0,"time_coverage_end")$value | |
datemean<-mean(c(as.Date(dateend,"%Y-%m-%dT%H:%M:%OSZ"),as.Date(dateini,"%Y-%m-%dT%H:%M:%OSZ"))) | |
year<-substring(datemean,0,4) | |
month<-substring(datemean,6,7) | |
# prepare final data set | |
dat.varSA<-data.frame(rep(as.integer(year,nrow(dat.varSAtmp))),rep(as.integer(month,nrow(dat.varSAtmp))), | |
dat.varSAtmp,rep(unit,nrow(dat.varSAtmp)),rep(var,nrow(dat.varSAtmp))) | |
names(dat.varSA)<-c("year","month","lon","lat","value","unit","var") | |
# save csv file | |
fe<-file.exists("MODISA_chl.csv") | |
write.table(dat.varSA,"MODISA_chl.csv",row.names=FALSE,col.names=!fe,sep=";",dec=",",append=fe) | |
# close netCDF | |
nc_close(data) | |
# clean workspace | |
rm(data,lon,lat,value,unit,dat.var,dat.varSAtmp,dateini,dateend,datemean,year,month,dat.varSA,fe) | |
} | |
rm(var,f,i,latmax,latmin,lf,lonmax,lonmin) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment