xarray to terra with reticulate

xarray-to-terra-with-reticulate.R

library(reticulate)
library(terra)

py_require(c("xarray", "netcdf4", "numpy"))

xr <- import("xarray")
np <- import("numpy")

url <- "https://psl.noaa.gov/thredds/dodsC/Datasets/noaa.oisst.v2.highres/sst.day.mean.2022.nc"

dataset <- xr$open_dataset(url)

# Select the sea surface temperature variable
# In this dataset, the variable is named "sst"
sst_var <- dataset$get("sst")

# Let's get a single time slice (first day in the dataset)
sst_day1 <- sst_var$sel(time=dataset$time[0])

# Convert to numpy array 
array_data <- sst_day1$values

# Get coordinate information
lons <- dataset$coords$get("lon")$values
lats <- dataset$coords$get("lat")$values

# Convert numpy array to R matrix
r_array <- py_to_r(array_data)

# Create a terra SpatRaster from the R array
rast <- rast(r_array)

# Set the extent based on longitude and latitude coordinates
# Note: This dataset uses 0-360 longitude, so we adjust to -180 to 180
lon_values <- py_to_r(lons)
lat_values <- py_to_r(lats)

# Handle the case where longitudes are in 0-360 range instead of -180 to 180
if (max(lon_values) > 180) {
  # Adjust coordinates to -180 to 180 range
  lon_values[lon_values > 180] <- lon_values[lon_values > 180] - 360
}

# Set the extent
ext(rast) <- c(min(lon_values), max(lon_values), 
               min(lat_values), max(lat_values))

# Set the CRS to WGS84
crs(rast) <- "EPSG:4326"

# Plot the sea surface temperature map
plot(rast, main="Sea Surface Temperature (Day 1 of 2022)", 
     col = hcl.colors(100, "RdBu", rev = TRUE))

# Save the raster if needed
# writeRaster(rast, "sst_day1_2022.tif", overwrite=TRUE)

# Clean up
dataset$close()