Web Content Display

Web Content Display (CoastWatch User Forum)

The NOAA CoastWatch Forum provides information in several categories.  The information is relevant to CoastWatch products and software (https://coastwatch.noaa.gov)  and vary in technical details.  For questions and issues not addressed here,  start a thread in the appropriate category or contact the CoastWatch HelpDesk (coastwatch.info@noaa.gov).  Please note our Privacy Policy (https://www.noaa.gov/privacy.html) and that this is a public forum.  All information is voluntarily and anonymous submissions are accepted pending review prior to posting. 

Forums

Back

RE: R script questions

U
Anonymous, modified 3 Years ago.

R script questions

Hi - I am trying to use this example to plot an SST timeseries: https://github.com/CoastWatch-WestCoast/r_code/blob/master/timeseries_satellite_data.md and I have modified the script up to this point: 

#set coordinates for Sitka area
xcoord<-c(-140, -130)
ycoord<-c(50,60)

##Format Box Coordinates for cosmetics, to make a nice map title
ttext<-paste(paste(abs(xcoord), collapse="-"),"W, ", paste(ycoord, collapse="-"),"N")

# Use rerddap to get information about the dataset
url="https://coastwatch.pfeg.noaa.gov/erddap/"
dataInfo <- rerddap::info('jplMURSST41', url=url)

# Display the dataset metadata
dataInfo

# Extract the parameter name from the metadata in dataInfo
parameter <- "analysed_sst"

# Set the altitude coordinate to zero
#zcoord <- 0.

# Extract the beginning and ending dates of the dataset from the metadata in dataInfo
global <- dataInfo$alldata$NC_GLOBAL
tt <- global[ global$attribute_name %in% c('time_coverage_end','time_coverage_start'), "value", ]

# Populate the time vector with the time_coverage_start from dataInfo
# Use the "last" option for the ending date
tcoord <- c(tt[2],"last")

# Extract the timeseries data using rxtracto_3D
sstMUR<-rxtracto_3D(dataInfo,
                    parameter=parameter,
                    tcoord=tcoord,
                    xcoord=xcoord,
                    ycoord=ycoord)

but I am getting the following error: 

 0s 0s[1] "error in trying to download the subset"                                                  
[1] "check your settings"
$x
<ERDDAP info> jplMURSST41 
 Base URL: https://coastwatch.pfeg.noaa.gov/erddap 
 Dataset Type: griddap 
 Dimensions (range):  
     time: (2002-06-01T09:00:00Z, 2022-01-25T09:00:00Z) 
     latitude: (-89.99, 89.99) 
     longitude: (-179.99, 180.0) 
 Variables:  
     analysed_sst: 
         Units: degree_C 
     analysis_error: 
         Units: degree_C 
     mask: 
     sea_ice_fraction: 
         Units: 1 

$url
[1] "https://coastwatch.pfeg.noaa.gov/erddap/"

$longitude
[1] -140 -130

$latitude
[1] 50 60

$time
[1] "2002-06-01T09:00:00Z" "2022-01-25T09:00:00Z"

$fields
[1] "analysed_sst"

$read
[1] FALSE

[1] "stopping execution  - will return what has been downloaded so far"
[1] "There was an error in the url call, perhaps a time out. See message on screen and URL called"
 0s
> 

I am not sure if I can make office hours but I will try - I figured I would post here just in case. Thanks!

-Rani Gaddam

MA
Melanie Abecassis, modified 3 Years ago.

RE: R script questions

Youngling Posts: 65 Join Date: 6/12/17 Recent Posts

Hi Rani,

The error messages are not always super useful, but what is happening here is that you are requesting too much data.

This a daily product, and you are requesting data for every day from June 2002 to now, which is a lot.

One way you can test this is by trying a much smaller range of dates: 

tcoord <- c(tt[2],'2002-06-15')
> sstMUR<-rxtracto_3D(dataInfo,
+                     parameter=parameter,
+                     tcoord=tcoord,
+                     xcoord=xcoord,
+                     ycoord=ycoord)

that worked for me.

If you request too much data, the server times out before it can finish the request.

Try playing with the date range to see how much you ask at once before it gives you an error. If you need such a long time series, you might need to explore weekly or monthly composites.

Hope this helps!

Melanie

CW
Cara Wilson, modified 3 Years ago.

RE: R script questions

Youngling Posts: 41 Join Date: 6/12/17 Recent Posts

I see Melanie and I responded at basically the same time, saying the same thing! 

CW
Cara Wilson, modified 3 Years ago.

RE: R script questions

Youngling Posts: 41 Join Date: 6/12/17 Recent Posts

Hi Rani - 

I think the problem is that you are asking for too much data. The MUR dataset is 1 km resolution, and asking for daily data at that resolution for a long time results in a lot of data! The limit is 2 GB per request.  For looking at 20 year timeseries daily data is usually overkill, so I suggest using the monthly data instead.  If you are still having an issue try testing the request with just a year of data by setting tcoord to something like c("2020-01-01", "2020-12-31"). 

I am having serious technical issues this week and currently do not have access to a computer with R so I cant test it directly myself. 

Cara 

MA
Melanie Abecassis, modified 3 Years ago.

RE: R script questions

Youngling Posts: 65 Join Date: 6/12/17 Recent Posts

Or you can chunk your request. ask for a month at a time, or maybe a whole year if it lets you.

U
Anonymous, modified 3 Years ago.

RE: R script questions

Thanks so much - that worked! I may limit to monthly since the period of time I am looking at is 2015-2017, so doing a month at a time may be a lot. Thanks again! - Rani

U
Anonymous, modified 3 Years ago.

RE: R script questions

OK - hoping I can make office hours in a bit but for now my graph is empty - likely it has something to do with the settings in the plot that are actually set up for chl...anyway this is what I have: ## Plot sstMUR
plot(as.Date(sstMUR$time), sstMUR$avg, 
     type='b', bg="blue", pch=21, xlab="", cex=.7,
     xlim=as.Date(c("2015-07-01","2015-07-31")),
     ylim=c(0,3),
     ylab="sst", main=ttext)

OK - hope I can sign on after 3 or so!!

U
Anonymous, modified 3 Years ago.

RE: R script questions

Hi - I have managed to compare 2 of the SST datasets but wanted to add the interpolated one and keep getting errors - this is the first one that has an altitude and I think I am assigning it wrong? I have tried zlev and zcoord but nothing seems to work:

 dataInfo
<ERDDAP info> ncdcOisst21Agg 
 Base URL: https://coastwatch.pfeg.noaa.gov/erddap 
 Dataset Type: griddap 
 Dimensions (range):  
     time: (1981-09-01T12:00:00Z, 2022-01-10T12:00:00Z) 
     zlev: (0.0, 0.0) 
     latitude: (-89.875, 89.875) 
     longitude: (0.125, 359.875) 
 Variables:  
     anom: 
         Units: degree_C 
     err: 
         Units: degree_C 
     ice: 
         Units: 1 
     sst: 
         Units: degree_C 
> # Set the altitude coordinate to zero
> zcoord <- 0.
> # Extract the beginning and ending dates of the dataset from the metadata in dataInfo
> global <- dataInfo$alldata$NC_GLOBAL
> tt <- global[ global$attribute_name %in% c('time_coverage_end','time_coverage_start'), "value", ]
> # Populate the time vector with the time_coverage_start from dataInfo
> # Use the "last" option for the ending date
> tcoord <- c('2014-09-01','2017-12-31')
> # Extract the timeseries data using rxtracto_3D
> sstOI<-rxtracto_3D(dataInfo,
+                      parameter=parameter,
+                      tcoord=tcoord,
+                      xcoord=xcoord,
+                      ycoord=ycoord,
+                      zcoord=zcoord)
[1] "Requested coordinate names do no match dataset coordinate names"
[1] "Requested coordinate names: longitude" "Requested coordinate names: latitude" 
[3] "Requested coordinate names: altitude"  "Requested coordinate names: time"     
[1] "Dataset coordinate names: time"      "Dataset coordinate names: zlev"     
[3] "Dataset coordinate names: latitude"  "Dataset coordinate names: longitude"

CW
Cara Wilson, modified 3 Years ago.

RE: R script questions

Youngling Posts: 41 Join Date: 6/12/17 Recent Posts

The problem here is that this dataset has an altitude dimension, but it is called zlev, not altitude, which is the default name for this dimension.  You just need to specify the name of that dimension by including zName='zlev' in the call

 

Cara 

U
Anonymous, modified 3 Years ago.

RE: R script questions

Hi Cara! I tried this but am still getting errors - I know I am missing something silly...I do have a graph with 2 datasets so not a big deal but it would be interesting to also see this one...

zName='zlev'
zcoord <- 0.

# Extract the beginning and ending dates of the dataset from the metadata in dataInfo
global <- dataInfo$alldata$NC_GLOBAL
tt <- global[ global$attribute_name %in% c('time_coverage_end','time_coverage_start'), "value", ]

# Populate the time vector with the time_coverage_start from dataInfo
# Use the "last" option for the ending date
tcoord <- c('2017-01-01','2017-12-31')

# Extract the timeseries data using rxtracto_3D
sstOI<-rxtracto_3D(dataInfo,
                     parameter=parameter,
                     tcoord=tcoord,
                     xcoord=xcoord,
                     ycoord=ycoord,
                     zcoord=zcoord)

MA
Melanie Abecassis, modified 3 Years ago.

RE: R script questions

Youngling Posts: 65 Join Date: 6/12/17 Recent Posts

For this, here's what worked for me:

url = "http://coastwatch.pfeg.noaa.gov/erddap/"
dataInfo <- rerddap::info('ncdcOisst21Agg',url=url)


> dataInfo
<ERDDAP info> ncdcOisst21Agg 
Base URL: http://coastwatch.pfeg.noaa.gov/erddap/ 
  Dimensions (range):  
  time: (1981-09-01T12:00:00Z, 2022-01-10T12:00:00Z) 
zlev: (0.0, 0.0) 
latitude: (-89.875, 89.875) 
longitude: (0.125, 359.875) 
Variables:  
  anom: 
  Units: degree_C 
err: 
  Units: degree_C 
ice: 
  Units: 1 
sst: 
  Units: degree_C 

tcoord <- c('2017-01-01','2017-12-31')
xcoord<-c(-140, -130)
ycoord<-c(50,60)
zcoord=c(0,0)
parameter <- "sst"

sstOI<-rxtracto_3D(dataInfo,
                   parameter=parameter,
                   tcoord=tcoord,
                   xcoord=xcoord,
                   ycoord=ycoord,
                   zName='zlev',
                   zcoord=zcoord)
 

U
Anonymous, modified 3 Years ago.

RE: R script questions

Oh great - I was able to extract the time series and calculate the averages but I am getting the following error when I try to plot - I can see the sstMUR and sstNOAA but not the sstOI - is it because I need to account for the altitude/zlev? 


> # Plot sstMUR
> plot(as.Date(sstMUR$time), sstMUR$avg, 
+      type='l', col="blue", pch=21, xlab="", cex=.7,
+      xlim=as.Date(c("2017-01-01","2017-12-31")),
+      ylim=c(0,20),
+      ylab="daily mean sst", main=ttext)
> 

> ## Add NOAA and OI data
> points(as.Date(sstNOAA$time), sstNOAA$avg, type='l', col="red", pch=21,cex=.7)
> points(as.Date(sstOI$time), sstOI$avg, type='l', col="black", pch=21,cex=.7)
Error in xy.coords(x, y) : 'x' and 'y' lengths differ

 

 

MA
Melanie Abecassis, modified 3 Years ago.

RE: R script questions

Youngling Posts: 65 Join Date: 6/12/17 Recent Posts

Hi Rani,

Got it!

because of that zlev extra varialble, the structure of sstOI$sst is different than the others.

dim(sstOI$sst)
[1]  41  41   1 365

You can compare it to sstMUR$sst and see how it differs.

So to get the average along the time axis, you need to change that line a bit:

sstOI$avg <- apply(sstOI$sst, c(4),function(x) mean(x,na.rm=TRUE))

The plotting should then work.

CW
Cara Wilson, modified 3 Years ago.

RE: R script questions

Youngling Posts: 41 Join Date: 6/12/17 Recent Posts

Rani - 

Alternatively you can drop the extraneous dimension: 

sstOI$sst <- drop(sstOI$sst)

which will make it a 3d matrix and then you don't have to change the line that does the averaging 

 

Cara

U
Anonymous, modified 3 Years ago.

RE: R script questions

Hi Cara, I tried this but got an error - maybe it dropped too much (since there are only 2 fields now)?

> sstOI$sst <- drop(sstOI$sst)
> dim(sstOI$sst)
[1]   2 365
> sstOI$avg <- apply(sstOI$sst, c(3),function(x) mean(x,na.rm=TRUE))
Error in apply(sstOI$sst, c(3), function(x) mean(x, na.rm = TRUE)) : 
  'MARGIN' does not match dim(X)

The other fix worked, but I figured I would try both - thanks! - Rani

U
Anonymous, modified 3 Years ago.

RE: R script questions

Thanks Melanie! That worked! I am going to try the other way too just to see, but I was able to make the 3rd graph! Looks pretty different from the other 2...interesting...

CW
Cara Wilson, modified 3 Years ago.

RE: R script questions

Youngling Posts: 41 Join Date: 6/12/17 Recent Posts

You have to specifiy Zname='zlev' in the function as Melanie shows in her post. ie 

sstOI<-rxtracto_3D(dataInfo,
                     parameter=parameter,
                     tcoord=tcoord,
                     xcoord=xcoord,
                     ycoord=ycoord,
                     zcoord=zcoord,

                      zName='zlev')

U
Anonymous, modified 3 Years ago.

RE: R script questions

Hi Cara! I did do that - and was able to extract the timeseries - here is a bigger chunk if it helps: 

# Extract the timeseries data using rxtracto_3D
sstOI<-rxtracto_3D(dataInfo,
                     parameter=parameter,
                     tcoord=tcoord,
                     xcoord=xcoord,
                     ycoord=ycoord,
                     zName='zlev',
                     zcoord=zcoord)

# Now spatially average the data into a timeseries
sstMUR$avg <- apply(sstMUR$analysed_sst, c(3),function(x) mean(x,na.rm=TRUE))
sstNOAA$avg <- apply(sstNOAA$CRW_SST, c(3),function(x) mean(x,na.rm=TRUE))
sstOI$avg <- apply(sstOI$sst, c(3),function(x) mean(x,na.rm=TRUE))

# Now temporally average the data into one map 
sstMUR$avgmap <- apply(sstMUR$analysed_sst,c(1,2),function(x) mean(x,na.rm=TRUE))
sstNOAA$avgmap <- apply(sstNOAA$CRW_SST,c(1,2),function(x) mean(x,na.rm=TRUE))
sstOI$avgmap <- apply(sstOI$sst,c(1,2),function(x) mean(x,na.rm=TRUE))


# Plot sstMUR
plot(as.Date(sstMUR$time), sstMUR$avg, 
     type='l', col="blue", pch=21, xlab="", cex=.7,
     xlim=as.Date(c("2017-01-01","2017-12-31")),
     ylim=c(0,20),
     ylab="daily mean sst", main=ttext)

## Add NOAA and OI data
points(as.Date(sstNOAA$time), sstNOAA$avg, type='l', col="red", pch=21,cex=.7)
points(as.Date(sstOI$time), sstOI$avg, type='l', col="black", pch=21,cex=.7)

 

Thanks! - Rani
 

U
Anonymous, modified 3 Years ago.

RE: R script questions

Hi, I am running the "Extract data within a boundary" part of the R-tutorial  but I keep getting the error about $ being invalid, e.g. here:

> sanctchl1$chla <- sanctchl1$chla[, , 2]

Error in sanctchl1$chla : $ operator is invalid for atomic vectors

 

MA
Melanie Abecassis, modified 3 Years ago.

RE: R script questions

Youngling Posts: 65 Join Date: 6/12/17 Recent Posts

Hi Rani,

Did you run the whole code without any changes?

I just did, and didn't get any issue. I can't replicate the error you got.

Melanie

U
Anonymous, modified 3 Years ago.

RE: R script questions

Hi Melanie - the chl post (from the tutorial) wasn't me! I am still trying to get the interpolated data to work for me but may just leave it out for now! - Rani

MA
Melanie Abecassis, modified 3 Years ago.

RE: R script questions

Youngling Posts: 65 Join Date: 6/12/17 Recent Posts

Oh ok, sorry. Someone added to your thread but didn't sign then. Hopefully they'll see the reply.

Reply to Main Thread
Quick Reply