Working with Aqua MODIS data in R
Data format:
Easiest format for starting are the Level 3 SMI files obtained through the browser at Ocean Color Web. There is very likely a way to automate retrieving data from within R, but I haven’t figured that out yet.
The files are HDF (version 4) files, which are not trivial to read in R, apparently (e.g. stackexchange). Fortunately there is a handy hdf-to-netcdf conversion tool at
http://hdfeos.org/software/h4cflib.php
Specifically h4tonccf_nc4 which can be made to work by doing
$ chmod 777 h4tonccf_nc4and then
$ h4tonccf_nc4 INPUTFILE.hdfto get a .nc file.
Reading in R
To load the converted netcdf file, we use the ncdf4 package, as well as the oce package (which has good routines for mapping and plotting).
install.packages('oce')
install.packages('ncdf4')As an example, let’s use the file A2014293.L3m_DAY_CHL_chlor_a_9km, which was unzipped from the downloaded file A2014293.L3m_DAY_CHL_chlor_a_9km.bz2 using the bunzip2 command, e.g.
$ bunzip2 A2014293.L3m_DAY_CHL_chlor_a_9km.bz2We convert it to netcdf
$ h4tonccf_nc4 A2014293.L3m_DAY_CHL_chlor_a_9km.bz2Note that the created file is called A2014293.nc. Now go to R for the rest. Load packages,
library(oce)## Loading required package: mapproj
## Loading required package: maps
## Loading required package: proj4library(ncdf4)and then load the netcdf file, and print its contents to summarize it:
nc <- nc_open('A2014293.nc')
nc## File A2014293.nc (NC_FORMAT_NETCDF4):
##
## 1 variables (excluding dimension variables):
## float l3m_data[Number_of_Columns,Number_of_Lines] (Chunking: [1440,720]) (Compression: shuffle,level 2)
## Fill: -32767
## Scaling: linear
## Scaling_Equation: (Slope*l3m_data) + Intercept = Parameter value
## Slope: 1
## Intercept: 0
## origname: l3m_data
## long_name: l3m_data
## scale_factor: 1
## add_offset: 0
## coordinates: Number_of_Lines Number_of_Columns
##
## 2 dimensions:
## Number_of_Columns Size:4320
## origname: longitude
## long_name: longitude
## coordinates: Number_of_Columns
## units: degrees_east
## standard_name: Longitude
## Number_of_Lines Size:2160
## origname: latitude
## long_name: latitude
## coordinates: Number_of_Lines
## units: degrees_north
## standard_name: Latitude
##
## 65 global attributes:
## Product_Name: A2014293.L3m_DAY_CHL_chlor_a_9km
## Sensor_Name: HMODISA
## Sensor:
## Title: HMODISA Level-3 Standard Mapped Image
## Data_Center:
## Station_Name:
## Station_Latitude: 0
## Station_Longitude: 0
## Mission:
## Mission_Characteristics:
## Sensor_Characteristics:
## Product_Type: O
## Processing_Version: 2013.1.1QL
## Software_Name: smigen
## Software_Version: 4.42
## Processing_Time: 2014294074631000
## Input_Files: A2014293.L3b_DAY_CHL.main
## Processing_Control: smigen par=A2014293.L3m_DAY_CHL_chlor_a_9km.param
## Input_Parameters: IFILE = A2014293.L3b_DAY_CHL.main|OFILE = A2014293.L3m_DAY_CHL_chlor_a_9km|PFILE = |PROD = chlor_a|PALFILE = DEFAULT|PROCESSING VERSION = 2013.1.1QL|MEAS = 1|STYPE = 0|DATAMIN = 0.000000|DATAMAX = 0.000000|LONWEST = -180.000000|LONEAST = 180.000000|LATSOUTH = -90.000000|LATNORTH = 90.000000|RESOLUTION = 9km|PROJECTION = RECT|GAP_FILL = 0|SEAM_LON = -180.000000|PRECISION=F
## L2_Flag_Names: ATMFAIL,LAND,HILT,HISATZEN,STRAYLIGHT,CLDICE,COCCOLITH,LOWLW,CHLWARN,CHLFAIL,NAVWARN,MAXAERITER,ATMWARN,HISOLZEN,NAVFAIL,FILTER,HIGLINT
## Period_Start_Year: 2014
## Period_Start_Day: 293
## Period_End_Year: 2014
## Period_End_Day: 293
## Start_Time: 2014293001508696
## End_Time: 2014294030006092
## Start_Year: 2014
## Start_Day: 293
## Start_Millisec: 908696
## End_Year: 2014
## End_Day: 294
## End_Millisec: 10806092
## Start_Orbit: -1
## End_Orbit: 66321
## Orbit: 66321
## Map_Projection: Equidistant Cylindrical
## Latitude_Units: degrees North
## Longitude_Units: degrees East
## Northernmost_Latitude: 90
## Southernmost_Latitude: -90
## Westernmost_Longitude: -180
## Easternmost_Longitude: 180
## Latitude_Step: 0.0833333358168602
## Longitude_Step: 0.0833333358168602
## SW_Point_Latitude: -89.9583358764648
## SW_Point_Longitude: -179.95832824707
## Data_Bins: 1468721
## Number_of_Lines: 2160
## Number_of_Columns: 4320
## Parameter: Chlorophyll a concentration
## Measure: Mean
## Units: mg m^-3
## Scaling: linear
## Scaling_Equation: (Slope*l3m_data) + Intercept = Parameter value
## Slope: 1
## Intercept: 0
## Data_Minimum: 0.00047999998787418
## Data_Maximum: 97.6829299926758
## Suggested_Image_Scaling_Minimum: 0.00999999977648258
## Suggested_Image_Scaling_Maximum: 20
## Suggested_Image_Scaling_Type: LOG
## Suggested_Image_Scaling_Applied: No
## _lastModified: 2014294074631000
## Conventions: CF-1.6
## institution: NASA/GSFC OBPGThe variables we want here are: l3m_data, Number_of_Columns (longitude), and Number_of_Lines (latitude). Extract these from the nc object:
chl <- ncvar_get(nc, 'l3m_data')
chl[chl < 0] <- NA # remove negative and missing values
lon <- ncvar_get(nc, 'Number_of_Columns')
lat <- ncvar_get(nc, 'Number_of_Lines')Now let’s plot it:
data(coastlineWorld)
## make a colormap that looks like the website
col <- colorRampPalette(c("purple", "#00007F", "blue",
"#007FFF", "cyan","#7FFF7F",
"yellow", "#FF7F00", "red", "#7F0000"))
imagep(lon, lat, log10(chl), col=col, zlim=log10(c(0.01, 20)), missingColor=1, zlab='log10(Chl)')
polygon(coastlineWorld[['longitude']], coastlineWorld[['latitude']], col='grey')