How can I convert Sentinel-3 LST data in NetCDF format to Geotiff format with coordinates using R?

129 Views Asked by At

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: enter image description here

However, I get strange results when plotting and saving as a raster using the script R: enter image description here

The data can be downloaded from here.

I appreciate it if anyone can identify the mistakes I made and indicate any solution.

0

There are 0 best solutions below