dankelley / oce

R package for oceanographic processing
http://dankelley.github.io/oce/
GNU General Public License v3.0
141 stars 42 forks source link

oce should read amsr satellite data #839

Closed dankelley closed 8 years ago

dankelley commented 8 years ago

@clayton33 and I have made a start on this. A question is whether to have it as a separate class, which would mean two satellite-related classes: landsat and amsr (to coin a name). Another option might be to have a larger class that handles both types, but we have had landsat for a while and I don't feel it is right to rename it. Also, the satellites are really quite different, so putting them all into one class is a bit like putting adp and adv into one class.

So, the proposal -- on which I'd like to get comments from @richardsc is whether the class name amsr makes sense, or whether he thinks we need to drop landsat and make something that has two varieties.

richardsc commented 8 years ago

I'm in favour of keeping a separate class, but the truth is that I know very little about AMSR data or methods.

This is a bit separate, but some time ago I helped someone read data from the MODIS satellite, and put my working in an Rpub at:

https://rpubs.com/clarkrichards/40319

dankelley commented 8 years ago

amsr is nice (in a way) because there are no headers at all. Things seem to be fixed, for a given satellite. I guess you can get the data the day after the satellite passes over. It's very nice to know about. We'll likely have a snapshot to show here by tomorrow sometime.

dankelley commented 8 years ago

This will be done in a new branch called amsr, which may very well not build at any time. (The idea is to permit me to work in bits and pieces, without harming the develop branch by leaving it in a broken state.)

CaraWilson commented 8 years ago

Can I make a comment on this? I am a satellite oceanographer, so I primarily work with satellite data, and my recent interest in the BioArgo data is a bit of an anomalous dip below the surface for me. Two questions/comments: 1) Why just AMSR data? I'm assuming you are referring to the SST data from AMSR data (it can also measure wind speed, rain fall and a few other parameters) and being a microwave measurement it can 'see' through clouds resulting in a more complete image. The drawback of microwave SST however is that the spatial resolution is lower than IR measurements (which can't see though clouds) and it can't make measurements too close to land. There is a slew of satellite data out there, and it seems a shame for you to spend much time coding to read just one type of data.
2) Are you familiar with ERDDAP? It's a web application 
and a web service designed to provide easier access to datasets for both people and machines which was developed at my lab, and our ERDDAP has many different satellite products on it (in addition to a lot of other data as well). See
http://coastwatch.pfeg.noaa.gov/erddap. While we have the AMSR-1 data on our ERDDAP (you can put AMSR into the search box on the top right of the ERDDAP page), we don't have the current AMSR data on ERDDAP, only because we serve other satellite SST products that blend AMSR with other SST products which are a better product in our opinion, like the MUSST products.
There is also an RERDDAP package being developed that reads data from ERDDAP into R.
Developing something to read satellite data from our ERDDAP site would be a better solution IMO that reading one specific AMSR file, but I also don't know the motivation behind this request.

dankelley commented 8 years ago

We (@clayton33 and I) are mostly doing it because a colleague at BIO suggested the dataset. Also, it's easy to decode ... a whole lot easier than landsat! I'll check into the website you mentioned.

One criterion is that it cover the past few years, and that it cover the geographical region called "Flemish Cap", near Newfoundland. Naturally, we would like better resolution and accuracy if we can get it!!

This is all just a few hours old.

The code is basically written, but not checked in to Github yet. It's easy to write because there is no header to the files, so it's just a matter of chopping out byte sequences, making matrices, and scaling from byte to physical values.

If I can find another that covers our time-space area, and if it's as easy to read as amsr, then I think I return to my earlier idea of making the function be called read.satellite ... we could abstract the existing read.landsat out from that.

THANKS very much for this really useful comment!!!

CaraWilson commented 8 years ago

Almost all our satellite data are global and go back as far as the dataset exists. The following url will generate an image around Newfoundland of the daily 1-km SST data set made by JPL, which goes back to June 2010:

http://coastwatch.pfeg.noaa.gov/erddap/griddap/jplG1SST.png?SST[(2010-06-09T12:00:00Z)][(38):(55)][(-66.):(-45)]&.draw=surface&.vars=longitude|latitude|SST&.colorBar=|||||&.bgColor=0xffccccff

The [(38):(55)][(-66.):(-45)] part of the url designates the latitude and longitude, and changing those values in the url will change the dimensions of the map generated.

Another GHRSST product, which is smoother, and goes back to 2002, is here: http://coastwatch.pfeg.noaa.gov/erddap/griddap/jplMURSST41.png?analysed_sst[(2016-01-19T09:00:00Z)][(38):(55)][(-66):(-45)]&.draw=surface&.vars=longitude|latitude|analysed_sst&.colorBar=|||||&.bgColor=0xffccccff

Changing the "png" in the url to "graph" will take you to a page where you can modify how the graph is displayed, and generate the corresponding url. Changing the "png" to "nc" will download a netcdf file of the data to your computer. There are also multiple other formats types available (.asc,.csv, .kml, .mat ect). A graph format will be displayed on your screen, a file format will be automatically downloaded.

These data can be read into R by changing the "png" to "nc" and removing the graphing information, i.e. the "&" and everything afterwards:


library(ncdf4)
test<-download.file('http://coastwatch.pfeg.noaa.gov/erddap/griddap/jplG1SST.nc?SST[(2010-06-09T12:00:00Z)][(38):(55)][(-66.):(-45)]',destfile="test.nc",mode='wb')
datafileID<-nc_open('test.nc') 

parameter<-as.vector(ncvar_get(datafileID))
longitude<-ncvar_get(datafileID, varid="longitude")
latitude<-ncvar_get(datafileID, varid="latitude")
datatime<-ncvar_get(datafileID, varid="time")
datatime<-as.Date(as.POSIXlt(datatime,origin='1970-01-01',tz= "GMT"))

SST<-array(parameter,dim=c(length(longitude),length(latitude),length(datatime)))

nc_close(datafileID)
dankelley commented 8 years ago

@CaraWilson - that first link fails for me. Apart from that, all I can say is I don't know whether this is a late Christmas or an early birthday for me! This coastwatch.pfeg.noaa.gov system seems to be fantastic!

CaraWilson commented 8 years ago

@dankelley - ha ha - I'm glad you like it! We are quite proud of ERDDAP, which was developed at my lab (ERD) by Bob Simons. It originally stood for Environmental Research Division Data Access Program, but Bob decided to de-acronym it about a year or so ago. There are about 50 others installed around the world. Instructions on how to install one are on the ERDDAP page, and Bob is always willing to help people out if they need help get one installed. Oceans Network Canada has one: http://dap.onc.uvic.ca/erddap/index.html

I think there might be something going on with that dataset, as the first link doesn't work for me either at the moment.

dankelley commented 8 years ago

The amsr branch (commit 2a24978c7065212f6f1e6384e6158fb54dbae921) now handles the following code, producing the graph shown. So this is a start. Possibly this will be useful even if the erddap solution proves to be better for our particular application.

PS. it strikes me that I should write an accessor, perhaps named [["SST"], that would average the day and night, to try to fill in some of those gaps.

library(oce)
## try(source("~/src/oce/R/amsr.R"))
d <- read.amsr("f34_20160102v7.2.gz")
## Test accessors (sensible for temperature?)
median(d[["SSTDay"]], na.rm=TRUE)
median(d[["SSTNight"]], na.rm=TRUE)
## Test summary (are they in right units?)
summary(d)

## Plot default and named channels
if (!interactive()) png("839a.png", pointsize=9)
par(mfrow=c(2,1))
plot(d) # default SSTDay
plot(d, 'SSTNight')
if (!interactive()) dev.off()

839a

CaraWilson commented 8 years ago

The problem with the jplG1SST dataset that I referenced above apparently is due to a change in the files on the JPL end. The power of ERDDAP is that you don't have to have all of the datasets at one place, but can rather reserve remotely served datasets, which is the case with this dataset. The jplG1SST dataset is currently (1/21/16) on the list of unavailable datasets (as of the last major dataset reload) at http://coastwatch.pfeg.noaa.gov/erddap/status.html then look for "n Datasets failed to load". We have contacted jpl to correct the problem, and it might take them a few days.

So this was an unfortunate timing for this example, as it was working earlier today...

CaraWilson commented 8 years ago

The issue with the jplG1SST dataset on our ERDDAP server has been partially corrected. The dates are bad in all the files since Jan. 13, so those files have been removed, but the earlier files (back to 2010) are available. Using last instead of the date in the url will produce a graph with the most recent data, i.e. (so if the link below produces a graph with data past 1/13/16 the issue has been resolved):

http://coastwatch.pfeg.noaa.gov/erddap/griddap/jplG1SST.png?SST%5B(last)%5D%5B(38):(55)%5D%5B(-66.):(-45)%5D&.draw=surface&.vars=longitude%7Clatitude%7CSST&.colorBar=%7C%7C%7C%7C%7C&.bgColor=0xffccccff

It's very embarrassing to have a broken link when trying to sell the merits of a system! But we are serving almost 1000 datasets on ERDDAP, and there will be occasional problems with some of them. Just my dumb luck that it would occur with this dataset just as I was posting a link using it...

dankelley commented 8 years ago

I am closing this because it (as defined in its title) has been addressed, with new capabilities being provided in the develop branch. I will check into the datasets mentioned by @CaraWilson and if I see a reasonable path to reading such data, I will open a new issue (and refer back here for Cara's comments).

CaraWilson commented 8 years ago

@dankelley - in addition to the amrs data you might wan to look at the geostationary SST data. Its IR data, so can't see through clouds, but because its on a geostationary satellite it gets hourly measurements. There are still cloud problems, but the spatial resolution is better than amrs. See

http://coastwatch.pfeg.noaa.gov/erddap/griddap/erdGAssta1day.graph?sst[(last)][(0.0)][(38):(55)][(295):(315)]&.draw=surface&.vars=longitude|latitude|sst&.colorBar=|||||&.bgColor=0xffccccff