CRS errors whed reading MODIS NetCDF file

Question

I am working with MODIS EVI 250m file in NetCDF format and trying to get the stack of EVI layers. So, I do:

# check nc-file for variables names
nc = nc_open(flist[1], write=FALSE, readunlim=TRUE, verbose=FALSE, auto_GMT=TRUE, suppress_dimvals=FALSE ) # open connection
nc # check attributes and var names
nc_close(nc) # close connection

and after

tmpin.s <- stack(flist[1], varname="_250m_16_days_EVI") 

and I see the error with CRS

> tmpin.s <- stack(flist[1], varname="_250m_16_days_EVI") 
[1] ">>>> WARNING <<<  attribute epsg_code is an 8-byte value, but R"
[1] "does not support this data type. I am returning a double precision"
[1] "floating point, but you must be aware that this could lose precision!"
[1] ">>>> WARNING <<<  attribute semi_major_axis is an 8-byte value, but R"
[1] "does not support this data type. I am returning a double precision"
[1] "floating point, but you must be aware that this could lose precision!"
Error in CRS(as.character(projection(crs))) : no colon in init= string

What am I doing wrong?

PS Here is the output of netCDF file

File D:/2019_R_ncr/MOD13Q1.006_250m_NCRfull.nc (NC_FORMAT_NETCDF4):

     3 variables (excluding dimension variables):
        byte crs[]   (Contiguous storage)  
[1] ">>>> WARNING <<<  attribute epsg_code is an 8-byte value, but R"
[1] "does not support this data type. I am returning a double precision"
[1] "floating point, but you must be aware that this could lose precision!"
[1] ">>>> WARNING <<<  attribute semi_major_axis is an 8-byte value, but R"
[1] "does not support this data type. I am returning a double precision"
[1] "floating point, but you must be aware that this could lose precision!"
            grid_mapping_name: latitude_longitude
            _CoordinateAxisTypes: GeoX GeoY
            epsg_code: 4326
            horizontal_datum_name: WGS84
            semi_major_axis: 6378137
            inverse_flattening: 298.257223563
            longitude_of_prime_meridian: 0
        short _250m_16_days_EVI[lon,lat,time]   (Chunking: [492,144,26])  (Compression: shuffle,level 4)
            _FillValue: -3000
            coordinates: time lat lon
            grid_mapping: crs
            valid_min: -2000
            valid_max: 10000
            long_name: 250m 16 days EVI
            units: EVI
            scale_factor_err: 0
            add_offset_err: 0
            calibrated_nt: 5
            scale_factor: 1e-04
            add_offset: 0
        int _250m_16_days_VI_Quality[lon,lat,time]   (Chunking: [389,121,21])  (Compression: shuffle,level 4)
            _FillValue: 65535
            coordinates: time lat lon
            grid_mapping: crs
            valid_min: 0
            valid_max: 65534
            long_name: 250m 16 days VI Quality
            units: bit field
            Legend: 
     Bit Fields Description (Right to Left): 
    [0-1] : MODLAND_QA [2 bit range] 
         00: VI produced, good quality 
         01: VI produced, but check other QA 
         10: Pixel produced, but most probably cloudy 
         11: Pixel not produced due to other reasons than clouds 
    [2-5] : VI usefulness [4 bit range] 
         0000: Highest quality 
         0001: Lower quality  
         0010..1010: Decreasing quality 
         1100: Lowest quality 
         1101: Quality so low that it is not useful 
         1110: L1B data faulty 
         1111: Not useful for any other reason/not processed 
    [6-7] : Aerosol quantity [2 bit range] 
         00: Climatology 
         01: Low 
         10: Average 
         11: High (11) 
    [8] : Adjacent cloud detected; [1 bit range] 
         1: Yes 
         0: No 
    [9] : Atmosphere BRDF correction performed [1 bit range] 
         1: Yes 
         0: No 
    [10] : Mixed clouds  [1 bit range] 
         1: Yes 
         0: No 
    [11-13] : Land/Water Flag [3 bit range] 
         000: Shallow ocean 
         001: Land (Nothing else but land) 
         010: Ocean coastlines and lake shorelines 
         011: Shallow inland water 
         100: Ephemeral water 
         101: Deep inland water 
         110: Moderate or continental ocean 
         111: Deep ocean 
    [14] : Possible snow/ice [1 bit range] 
         1: Yes 
         0: No 
    [15] : Possible shadow [1 bit range] 
         1: Yes 
         0: No 


     3 dimensions:
        time  Size:415
            standard_name: time
            axis: T
            calendar: julian
            units: days since 2000-01-01 00:00:00
        lat  Size:2298
            standard_name: latitude
            units: degrees_north
            _CoordinateAxisType: GeoY
            axis: Y
        lon  Size:7377
            standard_name: longitude
            units: degrees_east
            _CoordinateAxisType: GeoX
            axis: X

    7 global attributes:
        _NCProperties: version=1|netcdflibversion=4.4.1.1|hdf5libversion=1.8.18
        title: MOD13Q1.006 for aid0001
        Conventions: CF-1.6
        institution: Land Processes Distributed Active Archive Center (LP DAAC)
        source: AppEEARS v2.28
        references: See README.txt
        history: See README.txt

Answer 1

...looks like you need to define your crs within stack();

x<-stack('yourfile.nc', 
#               varname = '_250m_16_days_EVI', crs=CRS("+init=epsg:4326"))

if this doesnt work, try to process your data in modis-native crs

'+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs'

with AppEARS https://lpdaacsvc.cr.usgs.gov/appeears/ and reproject it with projectRaster() also check your Date-Format:

class(getZ(x))

Answer 2

You are not doing anything wrong. This is a bug in raster. I think this is now fixed in the development version (https://github.com/rspatial/raster).

You could try what consumere suggests

x <- stack(flist[i], varname = '_250m_16_days_EVI', crs=CRS("+init=epsg:4326"))

But that may fail also as the file provides a crs, or

b <- brick(flist[1], varname="_250m_16_days_EVI") 

But that will probably give the same error as well.

You can try with the terra package instead:

library(terra)
r <- rast(flist[1], "_250m_16_days_EVI") 

And continue with terra or return to raster like this:

b <- brick(r + 0)

(adding zero to remove the link to the .nc file).