stars::read_stars and raster::raster have different projections when reading the same .nc file

Miao Cai 0
Viewed 36

I'm trying to read a .nc raster file in R. The older function raster::raster() reads the data perfectly fine. I'd like to reproduce the results using a newer function stars::read_stars(), but somehow it does not work for me. The data CHAP_PM2.5_Y1K_2020_V4.nc can be downloaded from here (~6.1 MB). Below is a minimal reproducible example:

library(stars)
library(raster)

pm_raster = raster::raster('CHAP_PM2.5_Y1K_2020_V4.nc')
pm_stars = stars::read_stars('CHAP_PM2.5_Y1K_2020_V4.nc')
# Warning messages:
# 1: In CPL_read_gdal(as.character(x), as.character(options), as.character(driver),  :
#   GDAL Message 1: The dataset has several variables that could be identified as vector fields, 
# but not all share the same primary dimension. Consequently they will be ignored.
# 2: In CPL_read_gdal(as.character(x), as.character(options), as.character(driver),  :
#  GDAL Message 1: The dataset has several variables that could be identified as vector fields, 
# but not all share the same primary dimension. Consequently they will be ignored.

The reading the files looks good, but the problem is when I plot them, the figure read by stars::read_stars looks wrong:

plot(pm_raster, main = 'raster::raster')
plot(pm_stars, main = 'stars::read_stars')

It looks like the projection is wrong for stars::read_stars(), but I have no clue on this. Any suggestion or comment would be appreciated.

enter image description here

ncdf4 r-stars r-raster ggplot2 r
1 Answers

This is a non-answer but using slightly different read options. It looks like stars::read_ncdf is your better read option.

library(terra)
library(stars)
pm_terra <- rast('~/Downloads/CHAP_PM2.5_Y1K_2020_V4.nc')
pm_terra
class       : SpatRaster 
dimensions  : 3571, 6148, 1  (nrow, ncol, nlyr)
resolution  : 0.009999999, 0.01  (x, y)
extent      : 73.45216, 134.9322, 17.9697, 53.67971  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 
source      : CHAP_PM2.5_Y1K_2020_V4.nc 
varname     : PM2.5 
name        : PM2.5 
unit        : µg/m3 
plot(flip(pm_terra))

# and then using stars::read_ncdf
pm_stars_ncdf <- read_ncdf('~/Downloads/CHAP_PM2.5_Y1K_2020_V4.nc')
no 'var' specified, using PM2.5
other available variables:
 lat, lon
Will return stars object with 21954508 cells.
Warning message:
In .get_nc_projection(meta$attribute, rep_var, cv) :
  No projection information found in nc file.

pm_using_flip_terra

flip(terra) above

pm_using_stars::read_ncdf

stars::read_ncdf above.

userkHk30
Related