I have Sentinel-3 SLSTR LST data stored in different NetCDF files. I know how to extract the Land Surface Temperature, LST, from the Sentinal-3 through the SNAP, but, I would like to use R to work with this kind of data. The LST variable is stored in the "LST_in.nc" file. The latitude and longitude are in another file "geodetic_in.nc". Thus, I would like to convert Sentinel-3 LST NetCDF to GeoTiff format.
Here are my attempts:
dir <- "/home/user/S3A_SL_2_LST____20221125T125014_20221125T125314_20221126T214452_0179_092_266_3060_PS1_O_NT_004.SEN3/"
output_raster = "20221126T214452_0179_092_266_3060_PS1_O_NT_004"
# Loading libraries.
library(ncdf4)
library(raster)
library(dplyr)
library(ggplot2)
library(terra)
# Creating filepath names
climate_filepath <- paste0(dir, "LST_in.nc")
cart_filepath <- paste0(dir, "geodetic_in.nc")
# Reading them in using nc_open
nc <- nc_open(climate_filepath)
cord <- nc_open(cart_filepath)
# All three files have a 1200 x 1500 cell matrix. Thus, I collapsed the matrix, and bound them all into a dataframe:
latitude <- ncvar_get(cord, "latitude_in") %>% as.vector()
longitude <- ncvar_get(cord, "longitude_in") %>% as.vector()
lst <- ncvar_get(nc, "LST") %>% as.vector()
LST_DF = data.frame(lon = longitude,
lat = latitude,
LST = lst) %>%
#Convert from Kelvin to Celcius
dplyr::mutate(LST = LST - 273.15)
# I converted all the variables in a data frame to a matrix
LST_DF_matrix <- data.matrix(LST_DF, rownames.force = NA)
colnames(LST_DF_matrix) <- c('X', 'Y', 'Z')
head(LST_DF_matrix)
# Set up a raster geometry, using terra package
e <- ext(apply(LST_DF_matrix[,1:2], 2, range))
# Set up the raster. I choose this ncol and nrow, however, I don't know if it is correct.
# The dimension of the data is 1200, 1500, 1800000 (nrow, ncol, ncell) with the 1 km grid resolution. In another attempt, I opened each NetCDF as a raster too.
r <- rast(e, ncol=1000, nrow=1000)
# I apply rasterize
x <- rasterize(LST_DF_matrix[, 1:2], r, LST_DF_matrix[,3]) #, fun=mean)
plot(x) #, col = rev(rainbow(25)))
# Set CRS
crs(x) = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
# Saving to a GeoTIFF
writeRaster(x = x, filename = paste0(dir, output_raster, "_v3.tif"), overwrite=TRUE)
Using SNAP tool generated this raster, already with projection applied:
However, I get strange results when plotting and saving as a raster using the script R:
The data can be downloaded from here.
I appreciate it if anyone can identify the mistakes I made and indicate any solution.