Flipped raster map with no CRS and unable to mask with shapefile

I am trying to obtain the average temperature for different provinces in Ecuador on a specific date (2018-11-10).

I obtained the global temperature data from the CPC Global Unified Temperature website which comes in a .nc file for each year. Thus, I downloaded the data for 2018. However, when I extract the temperature from the array, the produced raster does not have a CRS crs : NA , extent : 0, 1, 0, 1 (xmin, xmax, ymin, ymax), the map is flipped, and I cannot project the raster CRS to that of the Ecuador shapefile in order to mask/crop the raster.

Below is the code and error I am getting. The data is available on my drive if you want to test the code.

# Load ncdf4 package
# install.packages('ncdf4')
library(ncdf4)

# Open the NetCDF file
tmax_2018 <- nc_open('data/temperature/tmax.2018.nc')

# Read the temperature variable
tmax <- ncvar_get(tmax_2018, "tmax")
time <- ncvar_get(tmax_2018, "time")
lon <- ncvar_get(tmax_2018, "lon")
lat <- ncvar_get(tmax_2018, "lat")

# Time is set in number of hours from '1900-1-1'
head(time)
as.Date(head(time)/24, origin = '1900-1-1')

> head(time)
[1] 1034376 1034400 1034424 1034448 1034472 1034496
> as.Date(head(time)/24, origin = '1900-1-1')
[1] "2018-01-01" "2018-01-02" "2018-01-03" "2018-01-04" "2018-01-05" "2018-01-06"

# Extract the temperature raster corresponding to '2018-11-10'
# Load raster package
library(raster)
tmax_Nov10 <- raster(tmax[, , which(as.Date(time/24, origin = '1900-1-1') == '2018-11-10')])

# Plot the global temperatures on '2018-11-10'
plot(tmax_Nov10)

> print(tmax_Nov10)
class      : RasterLayer 
dimensions : 720, 360, 259200  (nrow, ncol, ncell)
resolution : 0.002777778, 0.001388889  (x, y)
extent     : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      : layer 
values     : -38.12744, 44.14123  (min, max)

> tmax_Nov10 <- projectRaster(tmax_Nov10, '+proj=longlat +datum=WGS84 +no_defs')
Error in if (x@srs != "") { : missing value where TRUE/FALSE needed

As you can see the map is flipped, I have tried using the flip function but did not work. In the following code, I load the Ecuador shapefile but I am unable to mask/crop the raster.

# Load sf package
library(sf)

# Load Ecuador shapefiles
ecuador_shp <- st_read('data/ecu_shapefiles/ECU_adm1.shp')
ecuador_shp <- st_simplify(ecuador_shp, preserveTopology = TRUE, dTolerance = 100)

> st_geometry(ecuador_shp)
Geometry set for 24 features 
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -92.00854 ymin: -5.017001 xmax: -75.20007 ymax: 1.681093
Geodetic CRS:  WGS 84
First 5 geometries:
MULTIPOLYGON (((-78.55927 -2.552271, -78.5629 -...
MULTIPOLYGON (((-79.08675 -1.152232, -79.09587 ...
MULTIPOLYGON (((-79.08027 -2.288164, -79.08334 ...
MULTIPOLYGON (((-78.49581 1.18538, -78.49534 1....
MULTIPOLYGON (((-78.75167 -1.442297, -78.76531 ...

# Get the CRS
st_crs(ecuador_shp, as.Text = TRUE)

# ---- Mask the raster object with ecu_cantons_sp ----
# ecuador_shp to a Spatial object
ecuador_shp <- as(ecuador_shp, "Spatial")

# Transform tmax_Nov10 CRS to match ecuador_shp CRS (does not work as well)
tmax_Nov10 <- projectRaster(tmax_Nov10, st_crs(ecu_cantons, as.Text = TRUE))

> tmax_Nov10 <- projectRaster(tmax_Nov10, st_crs(ecu_cantons, as.Text = TRUE))
Error in if (x@srs != "") { : missing value where TRUE/FALSE needed

Solution

Instead of using ncdf4 to read the data, you can use raster or better terra functions, which get it in the right order with the right coordinates.

library(terra)
d = terra::rast("./tmax.2018.nc")
plot(d[[1]])

d is now a stack of 365 raster layers:

> d
class       : SpatRaster 
dimensions  : 360, 720, 365  (nrow, ncol, nlyr)
resolution  : 0.5, 0.5  (x, y)
extent      : 0, 360, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 
source      : tmax.2018.nc 
varname     : tmax (Daily Maximum Temperature) 
names       : tmax_1, tmax_2, tmax_3, tmax_4, tmax_5, tmax_6, ... 
unit        :   degC,   degC,   degC,   degC,   degC,   degC, ... 
time        : 2018-01-01 to 2018-12-31 UTC

and there's a time function:

> time(d)[1:10]
 [1] "2018-01-01 UTC" "2018-01-02 UTC" "2018-01-03 UTC" "2018-01-04 UTC"
 [5] "2018-01-05 UTC" "2018-01-06 UTC" "2018-01-07 UTC" "2018-01-08 UTC"
 [9] "2018-01-09 UTC" "2018-01-10 UTC"

You might want to use ncdf4 if there's any other useful metadata in the file but most of what you need gets read by terra:rast,