rxtractogon to raster ?

[ingomiller]

Hi,

I downloaded monthly SST data for a polygon using rxtractogon(). So far so good. What I'm struggling with is converting the output into a raster object as I need to add to a raster stack of other variables. This might be a silly question but I would really appreciate advise :)

Cheers

[rmendels]

@ingomiller It should be pretty easy (famous last words no?). Send me your code that did an extract and I will add on how to import to raster - this way I can test what I say before I send it.

[ingomiller]

Thanks for the quick response!

Here's the code to download one month of SST data for my study area:

SST_info <- rerddap::info('jplMURSST41', url = 'https://coastwatch.pfeg.noaa.gov/erddap/') SST_info str(SST_info)

SST_nc <- rxtractogon( dataInfo = SST_info, parameter = "sst", xcoord = c(136, 160), ycoord = c(-30, -5), tcoord = c("2023-12-01", "2023-12-31"), xName = "longitude", yName = "latitude", tName = "time", verbose = TRUE, cache_remove = TRUE )

I tried using tidy_grid() to create a data.frame and then transform it into a raster using raster::rasterFromXYZ(), but this did not work.

Once the transformation into a raster works, I have to write a loop as I need the monthly SST data since 2010.

Cheers

[rmendels]

@ingomiller some quick questions, since there are some bugs in your code and improvements can be made. Frst is there any reason you are using rxtractogon() as it doesn't look like you have defined a polygon. If it is a bounding box, rxtracto_3D() will be more efficient.

Second, do you really need the .01 degree resolution? That will produce a very large dataset over those years. If you can live with 0.25 degrees resolution there is a new MUR product that will be easier to deal with

[ingomiller]

Ok thanks for clarification. I used rxtracto_3D() as well. If possible I would like to stick with 0.01 resolution product as my other variables have a similar resolution.

[rmendels]

@ingomiller do you want individual rasters, a raster stack or a brick.

[ingomiller]

a raster brick preferably

[rmendels]

@ingomiller try the following, noting:

1. it uses rxtracto_3D()
2. the parameter name is "analyzed_sst"
3. conversion to brick is quick and dirty and likely there are better ways to do this (and do double check the result though at least it gets warmer as things go north:

SST_info <- rerddap::info('jplMURSST41', url = 'https://coastwatch.pfeg.noaa.gov/erddap/')
SST_nc <- rxtracto_3D(
dataInfo = SST_info,
parameter = "analysed_sst",
xcoord = c(136, 160),
ycoord = c(-30, -5),
tcoord = c("2023-12-01", "2023-12-31"),
xName = "longitude",
yName = "latitude",
tName = "time"
)
sst <- SST_nc$analysed_sst

# Coordinates and extents
xmn <- 136
xmx <- 160
ymn <- -30
ymx <- -5

# Prepare an empty list to hold raster layers
layers <- vector("list", 31)

# Create each layer from the 3D array and add it to the list
for (i in 1:31) {
  # Create a raster with the appropriate dimensions and extent
  r <- raster(nrows=2501, ncols=2401, xmn=xmn, xmx=xmx, ymn=ymn, ymx=ymx)

  # Assign the values from the 3D array, transposing the matrix as required
  values(r) <- c(t(sst[,,i]))

  # Add the raster to the list
  layers[[i]] <- r
}

# Convert the list of rasters into a brick
raster_brick <- brick(layers)

# Optionally, name the layers in the raster brick
names(raster_brick) <- paste0("Time_", 1:31)

# Check the raster brick
print(raster_brick)
plot(raster_brick[[1]]) # Example plot of the first layer

This produces (to see image look in Preview mode):

[rmendels]

@ingomiller - i think the raster part is wrong. Give me a minute and I will have I believe a simpler and correct way of doing it.

[rmendels]

@ingomiller Okay here is an easier way, using rerddap to get the data. Now the plot looks correct, and the code is simple:

library(raster)
library(rerddap)
sst_nc <- griddap(SST_info,
         longitude = c(136, 160),
         latitude = c(-30, -5),
         fields = "analysed_sst",
         time =  c("2023-12-01", "2023-12-31"),
         read = FALSE
         )
sst_nc_file <- sst_nc$summary$filename
sst_brick <- brick(sst_nc_file)
sst_brick
class      : RasterBrick 
dimensions : 2501, 2401, 6004901, 31  (nrow, ncol, ncell, nlayers)
resolution : 0.01, 0.01  (x, y)
extent     : 135.995, 160.005, -30.005, -4.995  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : 84937163baf76c392888c8b58279e940.nc 
names      : X2023.12.01.09.00.00, X2023.12.02.09.00.00, X2023.12.03.09.00.00, X2023.12.04.09.00.00, X2023.12.05.09.00.00, X2023.12.06.09.00.00, X2023.12.07.09.00.00, X2023.12.08.09.00.00, X2023.12.09.09.00.00, X2023.12.10.09.00.00, X2023.12.11.09.00.00, X2023.12.12.09.00.00, X2023.12.13.09.00.00, X2023.12.14.09.00.00, X2023.12.15.09.00.00, ... 
Date/time  : 2023-12-01 09:00:00, 2023-12-31 09:00:00 (min, max)
varname    : analysed_sst 
plot(sst_brick[[1]])

produces (and like what i get in ERDDAP directly:

[rmendels]

@ingomiller Did the last post solve your problem? I want to close this issue.

Thanks

[ingomiller]

Thank you so much! This does exactly what I needed and totally solved my issue :D Thanks heaps! Cheers