Trouble reading a netcdf file and convert it into a SpatRaster-CodePudding

I am struggling to open a NetCDF file and convert it into a raster using R. The data is supposed to be on a regular grid of 25 km by 25 km. It contains sea ice concentration in the Arctic.

library(terra)
#> terra 1.5.21
library(ncdf4)

file <- "~/Downloads/data_sat_Phil_changt_grid/SIC_SMMR_month_2015.nc"

I am getting a warning about the extent not found.

r <- rast(file)
#> Error in R_nc4_open: No such file or directory
#> Warning: [rast] GDAL did not find an extent. Cells not equally spaced?

We can see that there is a problem with the coordinates/extent.

r
#> class       : SpatRaster 
#> dimensions  : 448, 304, 12  (nrow, ncol, nlyr)
#> resolution  : 0.003289474, 0.002232143  (x, y)
#> extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 
#> source      : SIC_SMMR_month_2015.nc:sic 
#> varname     : sic 
#> names       : sic_1, sic_2, sic_3, sic_4, sic_5, sic_6, ...

I can open the nc file with nc_open() and I see that the coordinates are present.

nc <- nc_open(file)
names(nc$var)
#> [1] "lat" "lon" "sic"

lat <- ncvar_get(nc, "lat")
lon <- ncvar_get(nc, "lon")

dim(lat)
#> [1] 304 448
dim(lon)
#> [1] 304 448
dim(r)
#> [1] 448 304  12

Is it possible to assemble this data (the SIC values and the coordinates) to create a SpatRaster? The nc file can be downloaded here:

And also from

p = cbind(as.vector(lon), as.vector(lat))
plot(p, cex=.1, xlab="lon", ylab="lat")

So what you need to find out, is which coordinate reference system (crs) is used, clearly some kind of polar crs. And what the extent of the data set is.

From the website you point to, I take it we can use:

crs(r) = " proj=stere  lat_0=90  lat_ts=70  lon_0=-45  k=1  x_0=0  y_0=0  a=6378273  b=6356889.449  units=m  no_defs"
ext(r) = c(-3850000, 3750000, -5350000, 5850000)
r
#class       : SpatRaster 
#dimensions  : 448, 304, 12  (nrow, ncol, nlyr)
#resolution  : 25000, 25000  (x, y)
#extent      : -3850000, 3750000, -5350000, 5850000  (xmin, xmax, ymin, ymax)
#coord. ref. :  proj=stere  lat_0=90  lat_ts=70  lon_0=-45  x_0=0  y_0=0  a=6378273  b=6356889.449  units=m  no_defs 
#source      : SIC_SMMR_month_2015.nc:sic 
#varname     : sic 
#names       : sic_1, sic_2, sic_3, sic_4, sic_5, sic_6, ...

The results look good:

g = geodata::gadm("Greenland", level=0, path=".")
gg = project(g, crs(r))
plot(r,1)
lines(gg)

But this is of course not a good way to do such things; the ncdf file should have contained all the metadata required.