0

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)

enter image description here

> 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
Alonso Quijano
  • 73
  • 6

1 Answers1

2

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]])

enter image description here

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,

Spacedman
  • 92,590
  • 12
  • 140
  • 224
  • 1
    And perhaps use `r <- rotate(d)` to get the standard -180 to 180 longitude representation. – Robert Hijmans Jun 08 '23 at 14:12
  • Thanks. Because it was a .nc file, I didn't know I could use the `raster` or even the `terra` packages. I tried to replicate the map using your code, but the plots are not the same. Did you added some more code before plotting? Additionally, can I use the same logic with `terra` for masking the raster with the shapefiles that I have got as it can be done with `raster`. I am more familiar with `raster`, but I am looking forward to learning `terra` as well. Thanks. – Alonso Quijano Jun 08 '23 at 15:34
  • No, that's what I get with those three lines from a clean R session. What's different? The only other thing I've done is resize my graphics window. You can do masking with `terra`, if you get stuck doing that ask a new question! – Spacedman Jun 08 '23 at 17:25
  • When running the code you mentioned this is the output: `class : SpatRaster dimensions : 360, 720, 365 (nrow, ncol, nlyr) resolution : 1, 1 (x, y) extent : 0, 720, 0, 360 (xmin, xmax, ymin, ymax) coord. ref. : source : tmax.2018.nc names : tmax.2018_1, tmax.2018_2, tmax.2018_3, tmax.2018_4, tmax.2018_5, tmax.2018_6, ...` The CRS is still empty. Just in case, did you use the file I provided or downloaded another one? Thanks – Alonso Quijano Jun 13 '23 at 23:50
  • The data I use can be obtained from [here](https://psl.noaa.gov/repository/entry/show?entryid=synth%3Ae570c8f9-ec09-4e89-93b4-babd5651e7a9%3AL2NwY19nbG9iYWxfdGVtcC90bWF4LjIwMTgubmM%3D) – Alonso Quijano Jun 14 '23 at 00:04
  • UPDATE: I used R Studio Cloud to check if it was my R's problem and the map produced is the same as at yours. I tried uninstalling my terra package and installing it again but I can't get your map. :( – Alonso Quijano Jun 14 '23 at 00:15
  • UPDATE2: I uninstalled the terra package and installed it again from the github source (not cran) and I could get your map. Thanks! – Alonso Quijano Jun 14 '23 at 00:23
  • I tried to do the masking but I had some problems with the CRS, I asked a new question [here](https://stackoverflow.com/questions/76476464/unable-to-mask-spatial-object-even-after-transforming-the-crs-of-the-shapefile-t). Sorry for bothering and thanks for helping! – Alonso Quijano Jun 14 '23 at 18:34