raster and oceancolor L2 netcdf data

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
8 messages Options
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

raster and oceancolor L2 netcdf data

Warner, David
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

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: raster and oceancolor L2 netcdf data

Michael Sumner-2
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 <[hidden email]> 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
[hidden email]
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]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: raster and oceancolor L2 netcdf data

Warner, David
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 <[hidden email]> 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 <[hidden email]> 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
> [hidden email]
> 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

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: raster and oceancolor L2 netcdf data

Michael Sumner-2
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 <[hidden email]> 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 <[hidden email]>
> 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 <[hidden email]> 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
> [hidden email]
> 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]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: raster and oceancolor L2 netcdf data

Warner, David
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 <[hidden email]> 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 <[hidden email]> 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 <[hidden email]>
>> 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 <[hidden email]> 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
>> [hidden email]
>> 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
>
>


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

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: raster and oceancolor L2 netcdf data

Loïc Dutrieux-2
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 <[hidden email]> 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 <[hidden email]> 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 <[hidden email]>
>>> 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 <[hidden email]> 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
>>> [hidden email]
>>> 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
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: raster and oceancolor L2 netcdf data

Warner, David
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 <[hidden email]
> 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 <[hidden email]>
>> 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 <[hidden email]> 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 <[hidden email]>
>>>> 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 <[hidden email]> 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
>>>> [hidden email]
>>>> 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
> [hidden email]
> 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]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: raster and oceancolor L2 netcdf data

Michael Treglia
Hi All,

Just saw this thread, and wanted to share some code I've got, in case it's
useful, though I'll give the disclaimer I haven't looked carefully at the
focal dataset here, though it sounds somewhat similar.

Basically, I'm working with modeled data for NY Harbor, across multiple
time points, with multiple depth horizons.  It's also not a regular grid...
so I need look at the data as points or polygons. (it looks like a
spiderweb across the area...)

I've got code for pulling the data from a remote NetCDF source, assigning
values to x/y points, and plotting it with base R here:
https://github.com/mltConsEcol/R_NYHOPS_NetCDF_Access/blob/master/R/Working_RCode_MLT.R


You can disregard Ln 91 onward for now, but the rest might be of some use.
You can assign the data, with the x/y coords to a spatial (or simple
features) object - just havent included the code here yet (and sorry, in
the midst of travel so I'm appending this right now).

The next step will involve pulling x coords and y coords for polygons
associated with the points, which are represented in arrays (see lines
~66-71 of the code to get a sense of the structure), and then formatting
those such that they can be convert to spatial polygons/simple features
objects...  In the mean-time though, I've just used Voronoi polygons to
roughly approximate represent the areas associated with the points for
visualization (not included in the code, just done quickly in QGIS).

Again - hope this is of some use and not just noise... I'm all ears for
better ways to deal w/ the data too.

Best,
Mike




On Wed, Apr 12, 2017 at 6:02 PM, Warner, David <[hidden email]> wrote:

> 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 <
> [hidden email]
> > 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 <[hidden email]>
> >> 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 <[hidden email]> 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 <[hidden email]>
> >>>> 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 <[hidden email]> 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
> >>>> [hidden email]
> >>>> 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
> > [hidden email]
> > 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]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Loading...