0

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.

Can anyone identify the mistakes I made and indicate any solution?

1 Answer 1

0

Here is an approach using {stars} instead of {terra}. The code constructs a curvilinear grid (i.e a raster with cells having different sizes in x and y direction) and warps it to a regular grid (i.e. equal cell sizes) using st_warp(). You might change the parameters crs, cellsize, and threshold to your specific requirements.

library(stars)
library(sf)

coords <- read_ncdf("geodetic_in.nc", var = c("latitude_in", "longitude_in"))
lats <- coords$latitude_in
lons <- coords$longitude_in
lst <- read_ncdf("LST_in.nc", var = "LST")
lst <- lst$LST
lst <- st_as_stars(lst)
lst <- st_as_stars(lst, curvilinear = list(X1 = lons, X2 = lats))

lst_warp = st_warp(lst, crs = "EPSG:4326", cellsize = 0.01, threshold = 0.02)
write_stars(lst_warp, "lst_warped.tif")
1
  • Thank you very much for your answer, it worked out nicely! I applied it in my code and now GeoTIFF is generating perfectly without fail! I will study more about the stars package.
    – ammaciel
    May 29, 2023 at 14:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.