1

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
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))
| improve this answer | |
0

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

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.