Last active
January 4, 2023 16:57
-
-
Save aolinto/79e184f6c156c6ab21b3 to your computer and use it in GitHub Desktop.
R script to extract Aqua MODIS Sea Surface Temperature 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/79e184f6c156c6ab21b3 | |
# script to process Aqua MODIS Sea Surface Temperature | |
# files downloaded from http://oceancolor.gsfc.nasa.gov/cgi/l3 | |
# Aqua MODIS Sea Surface temperature 11 u daytime Monthly 9 km SMI images | |
# all .L3m_MO_SST_sst_9km.nc files must be in the working directory | |
# the script will open each nc file, read date, lon, lat and sst data, | |
# then select data from specified area and write them into | |
# a single csv file named MODISA_sst.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/MODIS/SST") # indicate the path to the files | |
file.exists("MODISA_sst.csv") # caution new data will be appended to this file if it already exists | |
# file.rename("MODISA_sst.csv","MODISA_sst.old") | |
# file.remove("MODISA_sst.csv") | |
# list and remove objects | |
ls() | |
rm(list = ls()) | |
# set the study area in decimal degrees | |
lonmax<--40 | |
lonmin<--50 | |
latmax<--22 | |
latmin<--29 | |
# create a list of files and indicate its length | |
(f <- list.files(".", pattern="*.L3m_MO_SST_sst_9km.nc",full.names=F)) | |
(lf<-length(f)) | |
# variable | |
var<-"sst" | |
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<-subset(dat.var,lon<=lonmax & lon>=lonmin & lat<=latmax & lat>=latmin & value<45) | |
# 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_sst.csv") | |
write.table(dat.varSA,"MODISA_sst.csv",row.names=FALSE,col.names=!fe,sep=";",dec=",",append=fe) | |
# close connection | |
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
Hello, Thanks for sharing your script
I try to extract data from several years and compare it with a buoy, so data MODISA vs data buoy, same data period, and frequency of data
I reduce de lonmax, lonmin, latmax, latmin until they are at the same point (pixel). It is an easy way to do this?
I take *sst4.nc data so I change the variable to "sst4"
but the analysis stop in the 3th file (see below)
[1] "Processing file 1 from 355 AQUA_MODIS.20150101_20150108.L3m.8D.SST4.x_sst4.nc"
[1] "Processing file 2 from 355 AQUA_MODIS.20150109_20150116.L3m.8D.SST4.x_sst4.nc"
[1] "Processing file 3 from 355 AQUA_MODIS.20150117_20150124.L3m.8D.SST4.x_sst4.nc"
Error in data.frame(rep(as.integer(year, nrow(dat.varSAtmp))), rep(as.integer(month, :
arguments imply differing number of rows: 1, 0
Could you help me to resolve this, please?