Skip to content
Prev 25594 / 29559 Next

raster and oceancolor L2 netcdf data

You can't georeference these data without remapping the data, essentially
treating the pixels as points. They have no natural regular grid form,
except possibly a unique satellite-perspective one. The data are in 2D
array form, but they have explicit "geolocation arrays", i.e. a longitude
and latitude for every cell and not based on a regular mapping.

R does not have tools for this directly from these data, though it can be
treated as a resampling or modelling problem.
You can use raster to get at the values of the locations easily enough,
here's a couple of plotting options in case it's useful:

u <- "
https://github.com/dmwarn/Tethys/blob/master/A2016244185500.L2_LAC_OC.x.nc?raw=true
"
f <- basename(f)
download.file(u, f, mode = "wb")

library(raster)
## use raster to treat as raw point data, on long-lat locations
rrs <- raster(f, varname = "geophysical_data/Rrs_412")
longitude <- raster(f, varname = "navigation_data/longitude")
latitude <- raster(f, varname = "navigation_data/latitude")

## plot in grid space, and add the geolocation space as a graticule
plot(rrs)
contour(longitude, add = TRUE)
contour(latitude, add = TRUE)

## raw scaling against rrs values
scl <- function(x) (x - min(x, na.rm = TRUE))/diff(range(x, na.rm = TRUE))
plot(values(longitude), values(latitude), pch = ".", col =
topo.colors(56)[scl(values(rrs)) * 55 + 1])

## ggplot
library(ggplot2)
d <- data.frame(x = values(longitude), y = values(latitude), rrs =
values(rrs))
ggplot(d, aes(x = x, y = y, colour = rrs)) + geom_point(pch = ".")

## might as well discard the missing values (depends on the other vars in
the file)
d <- d[!is.na(d$rrs), ]
ggplot(d, aes(x = x, y = y, colour = rrs)) + geom_point(pch = 19, cex = 0.1)

There are some MODIS and GDAL based packages that might be of use, but I
haven't yet seen any R tool that does this remapping task at scale. (I
believe the MODIS tools and the best warping tools in GDAL use thin-plate
spline techniques).

 Some applications would use the observations as points (i.e. the ocean
colour L3 bins as a daily aggregate from L2) and others use a direct
remapping of the data as an array, for say high-resolution sea ice imagery.

You might not need to do anything particularly complicated though, what's
the goal?

Cheers, Mike.
On Wed, Apr 12, 2017, 20:06 Warner, David <dmwarner at usgs.gov> wrote:
Greetings all

I am trying to develop R code for processing L2 data (netcdf v4 files) from
the Ocean Biology Processing Group.

The data file I am working with to develop the code is at
https://github.com/dmwarn/Tethys/blob/master/A2016244185500.L2_LAC_OC.x.nc
and represents primarily Lake Michigan in the United States.  The data were
extracted from a global dataset by the oceancolor L1 and L2 data server,
not by me.

I have been using the code below to try to get the data into R and ready
for processing but am having problems with dimensions and/or
orthorectification.  The

#extent of the scene obtained using nc_open and ncvar_get
nc <- nc_open('A2016214184500.L2_LAC_OC.x.nc')
lon <- ncvar_get(nc, "navigation_data/longitude")
lat <- ncvar_get(nc, "navigation_data/latitude")
minx <- min(lon)
maxx <- max(lon)
miny <- min(lat)
maxy <- max(lat)

#create extent object
myext <- extent(-90.817, -81.92438, 40.46493, 47.14244)

#create raster
rrs.412 <- raster('A2016214184500.L2_LAC_OC.x.nc', var
="geophysical_data/Rrs_412" ,
                  ext=myext)
rrs.412
class       : RasterLayer
dimensions  : 644, 528, 340032  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 0.5, 528.5, 0.5, 644.5  (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : /Users/dmwarner/Documents/MODIS/OC/
A2016214184500.L2_LAC_OC.x.nc
names       : Remote.sensing.reflectance.at.412.nm
zvar        : geophysical_data/Rrs_412

In spite of having tried to assign an extent, the raster extent is in rows
and columns.

Further, plotting the raster reveals that it is flipped on x axis and
somewhat rotated relative to what it should look like.  Even when flipped,
it is still not orthorectified.

How do I get the raster to have the correct extent and also get it
orthorectified?
Thanks,
Dave Warner

--
David Warner
Research Fisheries Biologist
U.S.G.S. Great Lakes Science Center
1451 Green Road
Ann Arbor, MI 48105
734-214-9392 <(734)%20214-9392>


_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo