Suggestions for irregular grid netCDF data?

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

Suggestions for irregular grid netCDF data?

Michael Treglia
​Hi all,

I've been trying to work with a somewhat complex netCDF dataset, and wanted
to see if anybody has suggestions.

The data are for marine variables in the form of a remote netCDF served
through OPENDAP. One of the complexities is that it comes in an irregular
grid, as fine as 25m and as coarse as 400m (coarser as you move further
from the coast) - this variable spacing makes it tough to work with in the
raster package, which seems like it would be the easiet.  Is there a
straightforward, relatively simple way to work with a dataset like this,
getting entire layers at a time? The data have a slew of different
variables, some of which also have time steps and different depth horizons.

So far, I’ve been able to connect to the data using the ncdf4 package in R,
and load simple, 3 dimensional data (e.g., depth). I can also force layers
to load using the raster package, adding the argument
‘stopIfNotEqualSpaced=FALSE', though that obviously has issues with the
irregular grid (though if it worked, would make it easier to pull
information by date and depth horizons). I'm having a tough time getting
data via the ncdf4 package for specific dates/depth horizons though. [note:
I end up having to use R in a linux environment for this, as I understand
the windows implementation for is more limited, and throws an ‘Unknown file
format’ error.]

The data are here:
http://colossus.dl.stevens-tech.edu/thredds/dodsC/LISS/Hindcast/all_nb_mon_81.nc,
and here’s the web-interface for the data:
http://colossus.dl.stevens-tech.edu/thredds/dodsC/LISS/Hindcast/all_nb_mon_81.nc.html.
I've included some sample R code below - and am totally open to alternative
ways that might be better for working with these data [still new-ish to
netCDFs, and have mainly dealt with much simpler, soil or climate datasets].

Thanks!
Mike


Sample R Code:
#Load and work w/ data in ncdf4
library(ncdf4)

nyhops.nc <- nc_open("
http://colossus.dl.stevens-tech.edu/thredds/dodsC/LISS/Hindcast/all_nb_mon_81.nc
")

#Get Point Data with Depth Measurements
depth.nc <- ncvar_get(nyhops.nc, 'depth')
lat.nc <- ncvar_get(nyhops.nc, 'lat')
long.nc <- ncvar_get(nyhops.nc, 'lon')

depth.df <- data.frame(lat=as.vector(lat.nc), long=as.vector(long.nc),
depth=as.vector(depth.nc)) #Can then be exported, or converted to spatial
object

#Trying to get the entirety of a single time point and depth horizon, but
not quite getting the start/count args - this gets a single data point it
seems
temp.nc <- ncvar_get(nyhops.nc, 'temp', start=c(10,10,1,1),
count=c(1,1,1,1))



#Load layers w/ raster package.
library(raster)
nyhops.url <- "
http://colossus.dl.stevens-tech.edu/thredds/dodsC/LISS/Hindcast/all_nb_mon_81.nc
"
depth <- raster(nyhops.url, varname="depth", lvar=1, level=1,
stopIfNotEqualSpaced=FALSE) #Forcing it to read layer despite unequal
spacing of coords
plot(depth) #Plot looks distorted, as expected given irregular grid
spacing. [I realize the data aren't projected at this point]

        [[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: Suggestions for irregular grid netCDF data?

Michael Sumner-2
Hi there,

my general advice is

* use raster to read the data layers (including the coordinates, which in
this context are just data)
* treat each layer "index space", i.e. setExtent(r, extent(0, nrow(x), 0,
ncol(x)) (or possibly transpose for some NetCDF files)
* use raster's cell-index logic to do extractions and aggregations in that
index space, it's neat because cropping inherits the parent georeference so
it's "global"
* use resampling to put real-world data objects into that index space for
extractions (you have to worry about the resolution and the boundary to do
that)

I have a few R packages that do this for ROMS and CMIP5 that can be
generally applied, but there not in great shape atm.

This rough document shows some of what you can do, it was written for a
particular purpose I had exploring leaflet which is on the backburner for
now, so it doesn't really go anywhere.

http://rpubs.com/cyclemumner/roms0

I'm happy to help you get what you need though, we deal with scads of this
kind of data and it doesn't fit any good mould except the "code it
yourself" club. It's not a huge amount of code but the "many coordinate
systems" can be bewildering and it's really specifically bound to the
intricacies of raster and sp and rgdal and ncdf4.

I'll try to look at your examples sometime soon.

Cheers, Mike.


On Thu, 2 Mar 2017 at 14:53 Michael Treglia <[hidden email]> wrote:

> ​Hi all,
>
> I've been trying to work with a somewhat complex netCDF dataset, and wanted
> to see if anybody has suggestions.
>
> The data are for marine variables in the form of a remote netCDF served
> through OPENDAP. One of the complexities is that it comes in an irregular
> grid, as fine as 25m and as coarse as 400m (coarser as you move further
> from the coast) - this variable spacing makes it tough to work with in the
> raster package, which seems like it would be the easiet.  Is there a
> straightforward, relatively simple way to work with a dataset like this,
> getting entire layers at a time? The data have a slew of different
> variables, some of which also have time steps and different depth horizons.
>
> So far, I’ve been able to connect to the data using the ncdf4 package in R,
> and load simple, 3 dimensional data (e.g., depth). I can also force layers
> to load using the raster package, adding the argument
> ‘stopIfNotEqualSpaced=FALSE', though that obviously has issues with the
> irregular grid (though if it worked, would make it easier to pull
> information by date and depth horizons). I'm having a tough time getting
> data via the ncdf4 package for specific dates/depth horizons though. [note:
> I end up having to use R in a linux environment for this, as I understand
> the windows implementation for is more limited, and throws an ‘Unknown file
> format’ error.]
>
> The data are here:
>
> http://colossus.dl.stevens-tech.edu/thredds/dodsC/LISS/Hindcast/all_nb_mon_81.nc
> ,
> and here’s the web-interface for the data:
>
> http://colossus.dl.stevens-tech.edu/thredds/dodsC/LISS/Hindcast/all_nb_mon_81.nc.html
> .
> I've included some sample R code below - and am totally open to alternative
> ways that might be better for working with these data [still new-ish to
> netCDFs, and have mainly dealt with much simpler, soil or climate
> datasets].
>
> Thanks!
> Mike
>
>
> Sample R Code:
> #Load and work w/ data in ncdf4
> library(ncdf4)
>
> nyhops.nc <- nc_open("
>
> http://colossus.dl.stevens-tech.edu/thredds/dodsC/LISS/Hindcast/all_nb_mon_81.nc
> ")
>
> #Get Point Data with Depth Measurements
> depth.nc <- ncvar_get(nyhops.nc, 'depth')
> lat.nc <- ncvar_get(nyhops.nc, 'lat')
> long.nc <- ncvar_get(nyhops.nc, 'lon')
>
> depth.df <- data.frame(lat=as.vector(lat.nc), long=as.vector(long.nc),
> depth=as.vector(depth.nc)) #Can then be exported, or converted to spatial
> object
>
> #Trying to get the entirety of a single time point and depth horizon, but
> not quite getting the start/count args - this gets a single data point it
> seems
> temp.nc <- ncvar_get(nyhops.nc, 'temp', start=c(10,10,1,1),
> count=c(1,1,1,1))
>
>
>
> #Load layers w/ raster package.
> library(raster)
> nyhops.url <- "
>
> http://colossus.dl.stevens-tech.edu/thredds/dodsC/LISS/Hindcast/all_nb_mon_81.nc
> "
> depth <- raster(nyhops.url, varname="depth", lvar=1, level=1,
> stopIfNotEqualSpaced=FALSE) #Forcing it to read layer despite unequal
> spacing of coords
> plot(depth) #Plot looks distorted, as expected given irregular grid
> spacing. [I realize the data aren't projected at this point]
>
>         [[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: Suggestions for irregular grid netCDF data?

Michael Treglia
Thanks Mike!

This is helpful, so far at least to see/think about different ways to work
with the data, as I wade through this structure and understand how
different R packages handle it.

I'll give things a go, and follow up, and will keep an eye out for any
other suggestions. And if you catch anything I'm missing if you go through
the example code, I'm all ears.  This isn't pressing on my end, but working
through this dataset (and similarly structured ones) will ultimately be
really useful.

Thanks again,
Mike T

Please pardon any typos, this message was sent from a mobile device.

On Mar 1, 2017 11:43 PM, "Michael Sumner" <[hidden email]> wrote:

> Hi there,
>
> my general advice is
>
> * use raster to read the data layers (including the coordinates, which in
> this context are just data)
> * treat each layer "index space", i.e. setExtent(r, extent(0, nrow(x), 0,
> ncol(x)) (or possibly transpose for some NetCDF files)
> * use raster's cell-index logic to do extractions and aggregations in that
> index space, it's neat because cropping inherits the parent georeference so
> it's "global"
> * use resampling to put real-world data objects into that index space for
> extractions (you have to worry about the resolution and the boundary to do
> that)
>
> I have a few R packages that do this for ROMS and CMIP5 that can be
> generally applied, but there not in great shape atm.
>
> This rough document shows some of what you can do, it was written for a
> particular purpose I had exploring leaflet which is on the backburner for
> now, so it doesn't really go anywhere.
>
> http://rpubs.com/cyclemumner/roms0
>
> I'm happy to help you get what you need though, we deal with scads of this
> kind of data and it doesn't fit any good mould except the "code it
> yourself" club. It's not a huge amount of code but the "many coordinate
> systems" can be bewildering and it's really specifically bound to the
> intricacies of raster and sp and rgdal and ncdf4.
>
> I'll try to look at your examples sometime soon.
>
> Cheers, Mike.
>
>
> On Thu, 2 Mar 2017 at 14:53 Michael Treglia <[hidden email]> wrote:
>
>> ​Hi all,
>>
>> I've been trying to work with a somewhat complex netCDF dataset, and
>> wanted
>> to see if anybody has suggestions.
>>
>> The data are for marine variables in the form of a remote netCDF served
>> through OPENDAP. One of the complexities is that it comes in an irregular
>> grid, as fine as 25m and as coarse as 400m (coarser as you move further
>> from the coast) - this variable spacing makes it tough to work with in the
>> raster package, which seems like it would be the easiet.  Is there a
>> straightforward, relatively simple way to work with a dataset like this,
>> getting entire layers at a time? The data have a slew of different
>> variables, some of which also have time steps and different depth
>> horizons.
>>
>> So far, I’ve been able to connect to the data using the ncdf4 package in
>> R,
>> and load simple, 3 dimensional data (e.g., depth). I can also force layers
>> to load using the raster package, adding the argument
>> ‘stopIfNotEqualSpaced=FALSE', though that obviously has issues with the
>> irregular grid (though if it worked, would make it easier to pull
>> information by date and depth horizons). I'm having a tough time getting
>> data via the ncdf4 package for specific dates/depth horizons though.
>> [note:
>> I end up having to use R in a linux environment for this, as I understand
>> the windows implementation for is more limited, and throws an ‘Unknown
>> file
>> format’ error.]
>>
>> The data are here:
>> http://colossus.dl.stevens-tech.edu/thredds/dodsC/LISS/
>> Hindcast/all_nb_mon_81.nc,
>> and here’s the web-interface for the data:
>> http://colossus.dl.stevens-tech.edu/thredds/dodsC/LISS/
>> Hindcast/all_nb_mon_81.nc.html.
>> I've included some sample R code below - and am totally open to
>> alternative
>> ways that might be better for working with these data [still new-ish to
>> netCDFs, and have mainly dealt with much simpler, soil or climate
>> datasets].
>>
>> Thanks!
>> Mike
>>
>>
>> Sample R Code:
>> #Load and work w/ data in ncdf4
>> library(ncdf4)
>>
>> nyhops.nc <- nc_open("
>> http://colossus.dl.stevens-tech.edu/thredds/dodsC/LISS/
>> Hindcast/all_nb_mon_81.nc
>> ")
>>
>> #Get Point Data with Depth Measurements
>> depth.nc <- ncvar_get(nyhops.nc, 'depth')
>> lat.nc <- ncvar_get(nyhops.nc, 'lat')
>> long.nc <- ncvar_get(nyhops.nc, 'lon')
>>
>> depth.df <- data.frame(lat=as.vector(lat.nc), long=as.vector(long.nc),
>> depth=as.vector(depth.nc)) #Can then be exported, or converted to spatial
>> object
>>
>> #Trying to get the entirety of a single time point and depth horizon, but
>> not quite getting the start/count args - this gets a single data point it
>> seems
>> temp.nc <- ncvar_get(nyhops.nc, 'temp', start=c(10,10,1,1),
>> count=c(1,1,1,1))
>>
>>
>>
>> #Load layers w/ raster package.
>> library(raster)
>> nyhops.url <- "
>> http://colossus.dl.stevens-tech.edu/thredds/dodsC/LISS/
>> Hindcast/all_nb_mon_81.nc
>> "
>> depth <- raster(nyhops.url, varname="depth", lvar=1, level=1,
>> stopIfNotEqualSpaced=FALSE) #Forcing it to read layer despite unequal
>> spacing of coords
>> plot(depth) #Plot looks distorted, as expected given irregular grid
>> spacing. [I realize the data aren't projected at this point]
>>
>>         [[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: Suggestions for irregular grid netCDF data?

Michael Treglia
Hi again,

Just to follow-up - got this working without too much issue... as of now, I
have the data values associated w/ the center lat/long of each grid cell.
There are also lat/long vertices for polygons associated with each grid
cell, but that'll take more working through it to get it together. However,
for now this largely suits my needs, and for now I can approximate the grid
cells using voronoi polygons for the grid cell centroids.

Here's some working code [I stop at creating dataframes, but going the next
steps for spatial objects is pretty easy - either sp or simple features
(still learning the latter)].

Thanks again!
Mike

#Sample Code
#Gets netCDF data for lat/long points
library(ncdf4)

#Connect to netCDF
nyhops.nc <- nc_open("
http://colossus.dl.stevens-tech.edu/thredds/dodsC/LISS/Hindcast/all_nb_mon_81.nc
")

##To get 'simple', 3D variale (e.g., depth)
#Load depth, lat, and long
depth.nc <- ncvar_get(nyhops.nc, 'depth')
lat.nc <- ncvar_get(nyhops.nc, 'lat')
long.nc <- ncvar_get(nyhops.nc, 'lon')

#Create dataframe
depth.df <- data.frame(lat=as.vector(lat.nc), long=as.vector(long.nc),
depth=as.vector(depth.nc))

#Remove NA values
depth.df <- depth.df[!is.na(depth.df$depth),]
#Can then convert to spatial object, export as table, etc.

##Example getting data involving a vertical horizon dimension and a date
dimension
#Fetch the temperature data; start is the starting position for x, y,
sigma, date; count is how many
# here, all rows & columns, and one horizon (horizon 3) and date (394)
temp.nc <- ncvar_get(nyhops.nc, 'temp', start=c(1,1,3,394),
count=c(147,452,1,1))

#Can Find out how many dimensions and entries:
nyhops.nc$var[['temp']]$varsize
#Output: [1] 147 452  11 394
#147 rows; 452 columns; 11 horizons; 394 dates/times

#Create dataframe w/ lat, long, temp
temp.df <- data.frame(cbind(long=as.vector(long.nc), lat=as.vector(lat.nc),
temp=as.vector(temp.nc)))

#remove NA rows for temp
temp.df <- temp.df[!is.na(temp.df$temp),] #Can conver to spatial object,
export as table, etc.

##Can plot, just as a non-spatial object for now.
#assign color ramp
rbPal <- colorRampPalette(c('blue',  'red'))

#Adds a column of color values based on the temp values
temp.df$Col <- rbPal(10)[as.numeric(cut(temp.df$temp,breaks = 10))]

plot(temp.df$long, temp.df$lat, col=temp.df$Col)


##Next step - use lon_vertices and lat_vertices to get polygons associated
with each point


On Thu, Mar 2, 2017 at 8:58 AM, Michael Treglia <[hidden email]> wrote:

> Thanks Mike!
>
> This is helpful, so far at least to see/think about different ways to work
> with the data, as I wade through this structure and understand how
> different R packages handle it.
>
> I'll give things a go, and follow up, and will keep an eye out for any
> other suggestions. And if you catch anything I'm missing if you go through
> the example code, I'm all ears.  This isn't pressing on my end, but working
> through this dataset (and similarly structured ones) will ultimately be
> really useful.
>
> Thanks again,
> Mike T
>
> Please pardon any typos, this message was sent from a mobile device.
>
> On Mar 1, 2017 11:43 PM, "Michael Sumner" <[hidden email]> wrote:
>
>> Hi there,
>>
>> my general advice is
>>
>> * use raster to read the data layers (including the coordinates, which in
>> this context are just data)
>> * treat each layer "index space", i.e. setExtent(r, extent(0, nrow(x), 0,
>> ncol(x)) (or possibly transpose for some NetCDF files)
>> * use raster's cell-index logic to do extractions and aggregations in
>> that index space, it's neat because cropping inherits the parent
>> georeference so it's "global"
>> * use resampling to put real-world data objects into that index space for
>> extractions (you have to worry about the resolution and the boundary to do
>> that)
>>
>> I have a few R packages that do this for ROMS and CMIP5 that can be
>> generally applied, but there not in great shape atm.
>>
>> This rough document shows some of what you can do, it was written for a
>> particular purpose I had exploring leaflet which is on the backburner for
>> now, so it doesn't really go anywhere.
>>
>> http://rpubs.com/cyclemumner/roms0
>>
>> I'm happy to help you get what you need though, we deal with scads of
>> this kind of data and it doesn't fit any good mould except the "code it
>> yourself" club. It's not a huge amount of code but the "many coordinate
>> systems" can be bewildering and it's really specifically bound to the
>> intricacies of raster and sp and rgdal and ncdf4.
>>
>> I'll try to look at your examples sometime soon.
>>
>> Cheers, Mike.
>>
>>
>> On Thu, 2 Mar 2017 at 14:53 Michael Treglia <[hidden email]> wrote:
>>
>>> ​Hi all,
>>>
>>> I've been trying to work with a somewhat complex netCDF dataset, and
>>> wanted
>>> to see if anybody has suggestions.
>>>
>>> The data are for marine variables in the form of a remote netCDF served
>>> through OPENDAP. One of the complexities is that it comes in an irregular
>>> grid, as fine as 25m and as coarse as 400m (coarser as you move further
>>> from the coast) - this variable spacing makes it tough to work with in
>>> the
>>> raster package, which seems like it would be the easiet.  Is there a
>>> straightforward, relatively simple way to work with a dataset like this,
>>> getting entire layers at a time? The data have a slew of different
>>> variables, some of which also have time steps and different depth
>>> horizons.
>>>
>>> So far, I’ve been able to connect to the data using the ncdf4 package in
>>> R,
>>> and load simple, 3 dimensional data (e.g., depth). I can also force
>>> layers
>>> to load using the raster package, adding the argument
>>> ‘stopIfNotEqualSpaced=FALSE', though that obviously has issues with the
>>> irregular grid (though if it worked, would make it easier to pull
>>> information by date and depth horizons). I'm having a tough time getting
>>> data via the ncdf4 package for specific dates/depth horizons though.
>>> [note:
>>> I end up having to use R in a linux environment for this, as I understand
>>> the windows implementation for is more limited, and throws an ‘Unknown
>>> file
>>> format’ error.]
>>>
>>> The data are here:
>>> http://colossus.dl.stevens-tech.edu/thredds/dodsC/LISS/Hindc
>>> ast/all_nb_mon_81.nc,
>>> and here’s the web-interface for the data:
>>> http://colossus.dl.stevens-tech.edu/thredds/dodsC/LISS/Hindc
>>> ast/all_nb_mon_81.nc.html.
>>> I've included some sample R code below - and am totally open to
>>> alternative
>>> ways that might be better for working with these data [still new-ish to
>>> netCDFs, and have mainly dealt with much simpler, soil or climate
>>> datasets].
>>>
>>> Thanks!
>>> Mike
>>>
>>>
>>> Sample R Code:
>>> #Load and work w/ data in ncdf4
>>> library(ncdf4)
>>>
>>> nyhops.nc <- nc_open("
>>> http://colossus.dl.stevens-tech.edu/thredds/dodsC/LISS/Hindc
>>> ast/all_nb_mon_81.nc
>>> ")
>>>
>>> #Get Point Data with Depth Measurements
>>> depth.nc <- ncvar_get(nyhops.nc, 'depth')
>>> lat.nc <- ncvar_get(nyhops.nc, 'lat')
>>> long.nc <- ncvar_get(nyhops.nc, 'lon')
>>>
>>> depth.df <- data.frame(lat=as.vector(lat.nc), long=as.vector(long.nc),
>>> depth=as.vector(depth.nc)) #Can then be exported, or converted to
>>> spatial
>>> object
>>>
>>> #Trying to get the entirety of a single time point and depth horizon, but
>>> not quite getting the start/count args - this gets a single data point it
>>> seems
>>> temp.nc <- ncvar_get(nyhops.nc, 'temp', start=c(10,10,1,1),
>>> count=c(1,1,1,1))
>>>
>>>
>>>
>>> #Load layers w/ raster package.
>>> library(raster)
>>> nyhops.url <- "
>>> http://colossus.dl.stevens-tech.edu/thredds/dodsC/LISS/Hindc
>>> ast/all_nb_mon_81.nc
>>> "
>>> depth <- raster(nyhops.url, varname="depth", lvar=1, level=1,
>>> stopIfNotEqualSpaced=FALSE) #Forcing it to read layer despite unequal
>>> spacing of coords
>>> plot(depth) #Plot looks distorted, as expected given irregular grid
>>> spacing. [I realize the data aren't projected at this point]
>>>
>>>         [[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
Loading...