[R-sig-Geo] raster and oceancolor L2 netcdf data

[Michael Sumner]

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
> 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>

        [[alternative HTML version deleted]]

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

-- 
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

	[[alternative HTML version deleted]]