CODEWITHSUNDEEP
rrastertemperaturegeotiffsentinel3
ammaciel
ammaciel

Reputation: 13

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

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.

Upvotes: 0

Views: 230

Answers (1)

Corbjn
Corbjn

Reputation: 286

The problem you're experiencing is, that your target grid has a higher resolution than your input data. When rasterising using terra::rasterize() values are aggregated if multiple cell values from the source fit into the same cell of the target grid. But think about the opposite case, where you have no cell value from source. In this case NA will be written to the target grid. I'm not sure how the SNAP tool is handling the geocoding but I suspect it is done via GCP and affine transformation

Upvotes: 0

Related Questions