How do I, correctly, reposition misplaced longitude and latitude in a RasterBrick from NETCDF file in R?

Question

I have a RasterBrick containing precipitation variable from TRMM (TMPA/3B43) Rainfall Estimate L3 1 month 0.25 degree x 0.25 degree V7 (TRMM_3B43. Extract from data is in the tif file.

This RasterBrick was generated from a NETCDF file using the following code

library(raster)

# set working directory    
setwd("C:/_data/TRMM_data/TRMM")

# import netcdf file into R using brick function in raster package
ppt.brick <- brick("TRMM_3B43_Aggregation.ncml.ncml.nc4")

# set coordinate resference system
crs(ppt.brick) = "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m"

The data ppt.brick had the information shown below:

> ppt.brick
class      : RasterBrick 
dimensions : 1440, 400, 576000, 252  (nrow, ncol, ncell, nlayers)
resolution : 0.25, 0.25  (x, y)
extent     : -50, 50, -180, 180  (xmin, xmax, ymin, ymax)
crs        : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs 
source     : C:/sijeh_data/TRMM_data/TRMM/TRMM_3B43_Aggregation.ncml.ncml.nc4 
names      : X1999.01.16, X1999.02.15, X1999.03.16, X1999.04.16, X1999.05.16, X1999.06.16, X1999.07.16, X1999.08.16, X1999.09.16, X1999.10.16, X1999.11.16, X1999.12.16, X2000.01.16, X2000.02.15, X2000.03.16, ... 
Date       : 1999-01-16, 2019-12-16 (min, max)
varname    : precipitation

plot(ppt.brick) produces the image below. But as you can see, latitudeplots on x-axis instead of y-axis and longitude on y-axis instead of x-axis .

The image should actually be as shown in from Image source

I have tried using the transpose function ppt.t <- t(ppt.brick) but this causes the RasterBrick to lose some information like the varname, date and time as shown below.

>ppt.t
class      : RasterBrick 
dimensions : 400, 1440, 576000, 252  (nrow, ncol, ncell, nlayers)
resolution : 0.25, 0.25  (x, y)
extent     : -180, 180, -50, 50  (xmin, xmax, ymin, ymax)
crs        : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs 
source     : C:/Users/saa815/AppData/Local/Temp/RtmpIHEl9R/raster/r_tmp_2021-05-19_141009_7792_98365.grd 
names      :  layer.1,  layer.2,  layer.3,  layer.4,  layer.5,  layer.6,  layer.7,  layer.8,  layer.9, layer.10, layer.11, layer.12, layer.13, layer.14, layer.15, ... 
min values :        0,        0,        0,        0,        0,        0,        0,        0,        0,        0,        0,        0,        0,        0,        0, ... 
max values : 1.837709, 1.792068, 2.171153, 1.309954, 1.687151, 2.159085, 3.079944, 2.278921, 1.547583, 1.631066, 2.687030, 1.929542, 1.488750, 1.950469, 1.307860, ...

Please, what is the best way to transpose a RasterBrick and still retain the information in the original data? It's also hard to tell if the transpose is done correctly because the axis that should be -50 become +50 after transpose. See image

I also need to mask the RasterBrick with a polygon afterward. So, I'll be grateful if you knew about that too.

Thank you

Answer 1

Here is a simple workflow that transposes and flips the data; and sets the correct extent. I downloaded a file from the source you point at, but these are HDF, not ncdf files. Otherwise all seems to be the same.

library(terra)
#terra version 1.2.13

f <- "3B43.20160701.7.HDF"
x <- rast(f)
#Warning message:
#[rast] unknown extent

y <- t(x)
z <- flip(y, "v")
ext(z) <- c(-180, 180, -50, 50)

z
#class       : SpatRaster 
#dimensions  : 400, 1440, 3  (nrow, ncol, nlyr)
#resolution  : 0.25, 0.25  (x, y)
#extent      : -180, 180, -50, 50  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#source      : memory 
#names       :           0,           1,           2 
#min values  : 0.000000000, 0.004333077, 0.000000000 
#max values  :   2.0519991,   0.5469698,  98.0000000 

plot(z, nc=1)

My understanding is that the first layer is "precipitation" (mm/hr), the second the "relative error", and the third the "gaugeRelativeWeighting".

If you have to do this multiple times, you can use a function like this (that also suppresses the warning; and only returns the first layer); and using |> for which you need R version 4.1

get_3B43 <- function(filename) {
   on.exit(options(warn=options()$warn))
   options(warn=-1)
   subf <- paste0('HDF4_SDS:UNKNOWN:"', filename, '":0"')
   x <- rast(subf) |> t() |> flip(direction="v")
   ext(x) <- c(-180, 180, -50, 50)
   time(x) <- as.Date(substr(filename, 6, 13), "%Y%m%d")
   names(x) <- "precipitation"
   units(x) <- "mm/hr" 
   x
}

r <- get_3B43(f)

r
#class       : SpatRaster 
#dimensions  : 400, 1440, 1  (nrow, ncol, nlyr)
#resolution  : 0.25, 0.25  (x, y)
#extent      : -180, 180, -50, 50  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#source      : memory 
#name        : precipitation 
#min value   :             0 
#max value   :      2.051999 
#unit        :         mm/hr 
#time        : 2016-07-01 

With the files you provided, I can do

ff <- list.files()
x <- lapply(ff, get_3B43)
x <- rast(x)
x
class       : SpatRaster 
dimensions  : 400, 1440, 12  (nrow, ncol, nlyr)
resolution  : 0.25, 0.25  (x, y)
extent      : -180, 180, -50, 50  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
sources     : memory  
              memory  
              memory  
              ... and 9 more source(s)
names       : preci~ation, preci~ation, preci~ation, preci~ation, preci~ation, preci~ation, ... 
min values  :           0,           0,           0,           0,           0,           0, ... 
max values  :    1.837709,    1.792068,    2.171153,    1.309954,    1.687151,    2.159085, ... 
unit        :       mm/hr,       mm/hr,       mm/hr,       mm/hr,       mm/hr,       mm/hr, ... 
time        : 1999-01-01 to 1999-12-01 
 

You may need to do this in batches (e.g. one year at a time, and then saving that to disk), perhaps like this

z <- writeCDF(x, "test.nc", varname="precipitation", unit="mm/hr")

Also, I do not know why you do this:

crs(ppt.brick) = "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m"

That is wrong, the CRS is lon/lat (this is not MODIS data).