Search code examples
rspatialr-rasterterra

rast() from Terra package won't extract spatial information form NetCDF file: [rast] unknown extent


I have been trying to load a NetCDF with the rast function from terra but I get the following warning message:

> SM_rast <- rast("Soil_moisture_v3.nc", "soil_moisture") 
Warning message:
[rast] unknown extent

The corresponding SpatRaster has no crs, wrong resolution and wrong coordinates:

> SM_rast
class       : SpatRaster 
dimensions  : 2160, 4320, 172  (nrow, ncol, nlyr)
resolution  : 1, 1  (x, y)
extent      : 0, 4320, 0, 2160  (xmin, xmax, ymin, ymax)
coord. ref. :  
source      : Soil_moisture_v3.nc 
names       : Soil_~_v3_1, Soil_~_v3_2, Soil_~_v3_3, Soil_~_v3_4, Soil_~_v3_5, Soil_~_v3_6, ... 

You can download the NetCDF file here.

Following the answer to this question I tried adding drivers="NETCDF" to the rast() function, but it couldn't even read the file:

> SM_rast <- rast("Soil_moisture_v3.nc", "soil_moisture", drivers="NETCDF") 
Error: [rast] cannot open this file as a SpatRaster: Soil_moisture_v3.nc
In addition: Warning message:
`Soil_moisture_v3.nc' not recognized as a supported file format. (GDAL error 4) 

When converting the SpatRaster to a dataframe, all cells containing NAs are missing and the coordinate columns are wrong:

> df_rast <- terra::as.data.frame(SM_rast[[1]], xy = TRUE)
> summary(df_rast)
       x                y          Soil_moisture_v3_1
 Min.   :   2.5   Min.   : 414.5   Min.   :0.02000   
 1st Qu.:1525.5   1st Qu.: 937.5   1st Qu.:0.07621   
 Median :2403.5   Median :1240.5   Median :0.17518   
 Mean   :2340.6   Mean   :1197.1   Mean   :0.20558   
 3rd Qu.:2841.5   3rd Qu.:1440.5   3rd Qu.:0.30225   
 Max.   :4317.5   Max.   :1960.5   Max.   :0.91196  

If I do the same using brick() from the raster package, everything works fine:

> SM_brick <- brick("data-raw/TROPOMI/Soil_moisture_v3.nc", varname = "soil_moisture")
> SM_brick
class      : RasterBrick 
dimensions : 2160, 4320, 9331200, 172  (nrow, ncol, ncell, nlayers)
resolution : 0.0833333, 0.0833333  (x, y)
extent     : -180, 179.9998, -90, 89.99993  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : Soil_moisture_v3.nc 
names      : X2018.02.02, X2018.02.10, X2018.02.18, X2018.02.26, X2018.03.06, X2018.03.14, X2018.03.22, X2018.03.30, X2018.04.07, X2018.04.15, X2018.04.23, X2018.05.01, X2018.05.09, X2018.05.17, X2018.05.25, ... 
Date       : 2018-02-02, 2021-10-24 (min, max)
varname    : soil_moisture 

> df_brick <- terra::as.data.frame(SM_brick[[1]], xy = TRUE)
> summary(df_brick)
       x                    y              X2018.02.02     
 Min.   :-179.95833   Min.   :-89.95834   Min.   :0        
 1st Qu.: -89.97920   1st Qu.:-44.97919   1st Qu.:0        
 Median :  -0.00008   Median : -0.00004   Median :0        
 Mean   :  -0.00008   Mean   : -0.00004   Mean   :0        
 3rd Qu.:  89.97905   3rd Qu.: 44.97911   3rd Qu.:0        
 Max.   : 179.95818   Max.   : 89.95826   Max.   :1        
                                          NA's   :8181031 
 

I still need to use terra because it's faster and has more useful functions not present in raster.


Solution

  • It all looks good to me when I do:

    library(terra)
    # terra 1.7.20
    rast("Soil_moisture_v3.nc")
    #class       : SpatRaster 
    #dimensions  : 2160, 4320, 172  (nrow, ncol, nlyr)
    #resolution  : 0.0833333, 0.0833333  (x, y)
    #extent      : -180, 179.9998, -90, 89.99993  (xmin, xmax, ymin, ymax)
    #coord. ref. : lon/lat WGS 84 (EPSG:4326) 
    #source      : Soil_moisture_v3.nc 
    #varname     : soil_moisture 
    #names       : soil_~ure_1, soil_~ure_2, soil_~ure_3, soil_~ure_4, soil_~ure_5, soil_~ure_6, ... 
    #unit        :       m3/m3,       m3/m3,       m3/m3,       m3/m3,       m3/m3,       m3/m3, ... 
    #time (days) : 2018-02-02 to 2021-10-24 
    
    gdal()
    #[1] "3.6.2"
    
    g <- gdal(drivers=T)
    g[g$name=="netCDF", ]
    #     name   type        can   vsi                  long.name
    #88 netCDF raster read/write FALSE Network Common Data Format
    

    If you do not see this, that may be because you have another version of terra, and/or GDAL, and/or because of your operating system and how you installed "terra".

    From the contents of your zipfile I can see that you are using a Mac. The problem may be related to the GDAL that comes with the binary version on CRAN. You may try to install from R-Universe instead with

    install.packages('terra', repos='https://rspatial.r-universe.dev')
    

    Or, if that does not work, install from source code. See the instructions.

    There is also this work-around:

    library(raster)
    # add zero to force the values out of the file.
    # that takes a while
    b <- brick("Soil_moisture_v3.nc") + 0
    r <- rast(b)