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

Wed Apr 12 22:12:47 CEST 2017

I've done some work recently on ocean color L2 binning/mapping, see the 
discussion on the ocean color forum 
https://oceancolor.gsfc.nasa.gov/forum/oceancolor/topic_show.pl?tid=6550 
and the code in the gist (which I'll update). It's not a R solution, but 
could be useful still.

Cheers,
Loïc

On 12/04/17 12:05, Warner, David wrote:
> Mike
> I had not really thought about order of operations to be honest.  I just
> noticed early on when I was attempting to use raster approach that the data
> were not mapped as hoped or orthorectified.  I certainly don't need to
> remap before calculating chlor-a on a daily basis as all the bands I need
> for a single day are aligned (if not mapped the way I wish).  In the end I
> do need the data correctly mapped as I want to do matchups with data
> collected with an LRAUV.
>
> I am planning on using locally calibrated coefficients.  I will check out
> your package!  I initially wanted to use L3 data but I and a colleague
> determined that there was for some reason poor agreement between the L3
> data and our in situ matchup data even though at L2 there is good
> agreement.  This colleague has typically done the heavy lifting using ENVI,
> which I don't have and would rather not learn if what I want to do can be
> done in R.
>
> It looks like I can create a raster with vect2rast.SpatialPoints() from the
> plotKML package quite easily but the default settings for cell.size lead to
> loss of data (I think).  You can set a cell.size but I am not sure if it
> works correctly without having multiple values per cell or not.  Or what it
> does if you have multiple values per cell.  There is some functionality
> that allows you to pick the first, last, the min, the max, or the mean if
> there are multiple  values for the same grid cell but I can't get that to
> work without Saga GIS.
>
> Cheers and thanks,
> Dave
>
> On Wed, Apr 12, 2017 at 8:57 AM, Michael Sumner <mdsumner at gmail.com> wrote:
>
>> 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
>>
>>
>
>