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

Warner, David dmwarner at usgs.gov
Thu Apr 13 00:02:58 CEST 2017 

Thanks Loïc!  I will read this.  I am not completely wedded to doing what I
want in R either.  Python would be ok but I have to learn more about
programming with Python.

Cheers,
Dave



On Wed, Apr 12, 2017 at 4:12 PM, Loïc Dutrieux <loic.dutrieux at conabio.gob.mx
> wrote:

> I've done some work recently on ocean color L2 binning/mapping, see the
> discussion on the ocean color forum https://oceancolor.gsfc.nasa.g
> ov/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
>>>
>>>
>>>
>>
>>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>


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

	[[alternative HTML version deleted]]

More information about the R-sig-Geo mailing list