2

I'm working with MODIS Level 2 Chlorophyll data in R. Specifically, I'm analyzing Chlorophyll data offshore of the Channel Islands, particularly around San Nicolas Island (SNI). The data are publicly available, and the exact file can be downloaded here.

I've encountered an odd issue where the Chlorophyll data are inconsistent with the latitude / longitude locations within the .nc file.

I first extract information from the .nc file:

library("raster")
library("ncdf4")
library("rgdal")
library("sp")

# open .nc file
nc_data <- nc_open("T2020165182000.L2_LAC_OC.nc")

# extract chlorophyll data and latitude / longitude coordinates
chlor <- ncvar_get(nc_data, "geophysical_data/chlor_a")
lat <- as.data.frame(ncvar_get(nc_data, "navigation_data/latitude"))
long <- as.data.frame(ncvar_get(nc_data, "navigation_data/longitude"))

# create raster of chlorophyll data bounded by min and max of latitude / longitude
p1 <- raster(t(chlor),
xmn=min(long), xmx=max(long),
ymn=min(lat), ymx=max(lat))

plot(p1, zlim=c(0,10))

now I add coastlines (source of the coastline data included here to be thorough):

# download the data 
download.file("http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/physical/ne_10m_coastline.zip", destfile = 'coastlines.zip')

# unzip the file
unzip(zipfile = "coastlines.zip", exdir = 'ne-coastlines-10m')

# load the data 
coastlines <- readOGR("ne-coastlines-10m/ne_10m_coastline.shp")

# add coastlines to existing plot, and add lines to intersect at SNI
plot(coastlines, add=TRUE)
abline(v=-119.4992, h=33.2465, col="blue")

we can now see that the .nc file's latitude and longitude coordinates are spot on, i.e., the blue lines intersect precisely at where SNI really is (33.2465 N, -119.4992 W) though it's hard to see the small island on the figure.

BUT, crucially, the Chlorophyll data seem to be offset from their own latitude / longitude coordinates??

enter image description here

edited to note that the coastline layer isn't necessary here . . . it simply visually confirms that the .nc lat / long are "correct" -- the intersecting lines will correctly identify the location of SNI even if the only layer defined is:

p1 <- raster(t(chlor),
xmn=min(long), xmx=max(long),
ymn=min(lat), ymx=max(lat))

plot(p1, zlim=c(0,10))
abline(v=-119.4992, h=33.2465, col="blue")

In fact, if the xmin, xmax, etc., are alternatively not defined from the beginning, the Chlorophyll data still plot exactly as shown above, only with both axes defaulting to 0-1.

Furthermore, note that I've loaded this exact data file into SeaDAS, a NOAA program, and the coastlines and the Chlorophyll data all plot correctly, thus I don't think there's an issue with these data. SO, I think something is going awry with:

p1 <- raster(t(chlor))

Can someone point out what I'm doing wrong, or what I need to do?

1

You are treating the data as if it were a regular raster, but it is not. See:

library(ncdf4)
nc_data <- nc_open("T2020165182000.L2_LAC_OC.nc")
chlor <- ncvar_get(nc_data, "geophysical_data/chlor_a")
lat <- ncvar_get(nc_data, "navigation_data/latitude")
long <- ncvar_get(nc_data, "navigation_data/longitude")

xy <- cbind(as.vector(long), as.vector(lat))
# a sample for faster plotting. 
i <- seq(1, nrow(xy), 500)
plot(xy[i,])

enter image description here

You could fix this by rasterizing the points to a raster (and pick an appropriate resolution to avoid gaps)

library(raster)
r <- raster(extent(xy), res=.1)
x <- rasterize(xy, r, as.vector(chlor), mean)
plot(x)
lines(coastlines)

enter image description here

1
  • Thank you, Robert, I really appreciate your assistance! That solved it beautifully! :-)
    – zhrandell
    Jun 6 at 18:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.