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

Wed Apr 12 14:57:42 CEST 2017

Glad it's some help, but it sounds like you intend to calculate after
mapping (?) which is definitely not the right way to go. Calculate
chlorophyll and then map, that's how Seadas does it, even though the
remapping is the hard part. And apologies if I misread,  just checking.

I have two algorithms in my roc package on GitHub in case they help
understanding how the calcs get done. Presumably you'll have locally
calibrated parameters for a local algo?

If you want to aggregate into a local map I think it's fair to group-by
directly on L2 pixels coords and then sum into a geographic map, without
worrying about swath-as-image at all. We've road tested doing this but want
the entire southern ocean eventually so it needs a bit of unrelated
preparation for the raw files.

I'd be happy to explore an R solution off list if you're interested. L2 is
surprisingly easy and efficient in R via GDAL.

(This is also a good example for future workflows for the planned stars
package imo.)

Cheers, Mike

On Wed, Apr 12, 2017, 22:35 Warner, David <dmwarner at usgs.gov> wrote:

> Thanks Mike!
>
> The goal is to estimate daily chlorophyll via band ratio polynomial
> equation for hundreds of days of data (hundreds of data files).  Sounds
> like rather than finding a way to orthorectify in R I should learn to batch
> reproject using SeaDAS, which does produce a product that is in geotiff
> format, is orthorectified, and has readily mappable.  I was trying to avoid
> that as the help and documentation available for doing that seems much less
> abundant.  One file at a time is easy using the SeaDAS gui.
>
> Thanks very, very much for the other tricks.  Not surprisingly, ggplot2
> comes through again with plots that look right!
> Cheers,
> Dave
>
>
>
> On Wed, Apr 12, 2017 at 7:01 AM, Michael Sumner <mdsumner at gmail.com>
> wrote:
>
> 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
>
>
>
>
> --
> David Warner
> Research Fisheries Biologist
> U.S.G.S. Great Lakes Science Center
> 1451 Green Road
> Ann Arbor, MI 48105
> 734-214-9392
>
-- 
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

	[[alternative HTML version deleted]]