Skip to content

Raster: Can't read raster-created NCDF back in

2 messages · Mark Payne, Michael Sumner

Mark Payne
# 
Hi,

I have a rasterBrick object that I have created through a series of
manipulations and written to disk using writeRaster(x,format="CDF"). In
another, independent script, I then need to read that netcdf file back in.
The following commands work:
but these don't:
Error in ncvar_get_inner(ncid2use, varid2use, nc$var[[li]]$missval,
addOffset,  :
  Error: variable has 2 dims, but start has 3 entries.  They must match!
Error in ncvar_get_inner(ncid2use, varid2use, nc$var[[li]]$missval,
addOffset,  :
  Error: variable has 2 dims, but start has 3 entries.  They must match!
which is just strange. Does anyone have any ideas what might be going on
here?

The file is available here (its 60kb)
http://bit.ly/1TDEJHZ

The netcdf header and my session info follow.

Cheers,

Mark




netcdf Bluefin_OISST_clim {
dimensions:
    longitude = 176 ;
    latitude = 80 ;
variables:
    double longitude(longitude) ;
        longitude:units = "degrees_east" ;
        longitude:long_name = "longitude" ;
    double latitude(latitude) ;
        latitude:units = "degrees_north" ;
        latitude:long_name = "latitude" ;
    float SST(latitude, longitude) ;
        SST:_FillValue = -3.4e+38 ;
        SST:missing_value = -3.4e+38 ;
        SST:long_name = "SST" ;
        SST:projection = "+proj=longlat +datum=WGS84 +ellps=WGS84
+towgs84=0,0,0" ;
        SST:projection_format = "PROJ.4" ;
        SST:min = 0.331762666431522 ;
        SST:max = 16.0566702570234 ;

// global attributes:
        :Conventions = "CF-1.4" ;
        :created_by = "R, packages ncdf4 and raster (version 2.5-2)" ;
        :date = "2016-05-25 11:47:53" ;
data:

 longitude = 308.125, 308.375, 308.625, 308.875, 309.125, 309.375, 309.625,
    309.875, 310.125, 310.375, 310.625, 310.875, 311.125, 311.375, 311.625,
    311.875, 312.125, 312.375, 312.625, 312.875, 313.125, 313.375, 313.625,
    313.875, 314.125, 314.375, 314.625, 314.875, 315.125, 315.375, 315.625,
    315.875, 316.125, 316.375, 316.625, 316.875, 317.125, 317.375, 317.625,
    317.875, 318.125, 318.375, 318.625, 318.875, 319.125, 319.375, 319.625,
    319.875, 320.125, 320.375, 320.625, 320.875, 321.125, 321.375, 321.625,
    321.875, 322.125, 322.375, 322.625, 322.875, 323.125, 323.375, 323.625,
    323.875, 324.125, 324.375, 324.625, 324.875, 325.125, 325.375, 325.625,
    325.875, 326.125, 326.375, 326.625, 326.875, 327.125, 327.375, 327.625,
    327.875, 328.125, 328.375, 328.625, 328.875, 329.125, 329.375, 329.625,
    329.875, 330.125, 330.375, 330.625, 330.875, 331.125, 331.375, 331.625,
    331.875, 332.125, 332.375, 332.625, 332.875, 333.125, 333.375, 333.625,
    333.875, 334.125, 334.375, 334.625, 334.875, 335.125, 335.375, 335.625,
    335.875, 336.125, 336.375, 336.625, 336.875, 337.125, 337.375, 337.625,
    337.875, 338.125, 338.375, 338.625, 338.875, 339.125, 339.375, 339.625,
    339.875, 340.125, 340.375, 340.625, 340.875, 341.125, 341.375, 341.625,
    341.875, 342.125, 342.375, 342.625, 342.875, 343.125, 343.375, 343.625,
    343.875, 344.125, 344.375, 344.625, 344.875, 345.125, 345.375, 345.625,
    345.875, 346.125, 346.375, 346.625, 346.875, 347.125, 347.375, 347.625,
    347.875, 348.125, 348.375, 348.625, 348.875, 349.125, 349.375, 349.625,
    349.875, 350.125, 350.375, 350.625, 350.875, 351.125, 351.375, 351.625,
    351.875 ;

 latitude = 71.875, 71.625, 71.375, 71.125, 70.875, 70.625, 70.375, 70.125,
    69.875, 69.625, 69.375, 69.125, 68.875, 68.625, 68.375, 68.125, 67.875,
    67.625, 67.375, 67.125, 66.875, 66.625, 66.375, 66.125, 65.875, 65.625,
    65.375, 65.125, 64.875, 64.625, 64.375, 64.125, 63.875, 63.625, 63.375,
    63.125, 62.875, 62.625, 62.375, 62.125, 61.875, 61.625, 61.375, 61.125,
    60.875, 60.625, 60.375, 60.125, 59.875, 59.625, 59.375, 59.125, 58.875,
    58.625, 58.375, 58.125, 57.875, 57.625, 57.375, 57.125, 56.875, 56.625,
    56.375, 56.125, 55.875, 55.625, 55.375, 55.125, 54.875, 54.625, 54.375,
    54.125, 53.875, 53.625, 53.375, 53.125, 52.875, 52.625, 52.375, 52.125 ;
}
R version 3.2.2 (2015-08-14)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 14.04.1 LTS

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
LC_TIME=en_NZ.UTF-8
 [4] LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_NZ.UTF-8
LC_MESSAGES=en_US.UTF-8
 [7] LC_PAPER=en_NZ.UTF-8       LC_NAME=C
LC_ADDRESS=C
[10] LC_TELEPHONE=C             LC_MEASUREMENT=en_NZ.UTF-8
LC_IDENTIFICATION=C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] lubridate_1.5.6 stringr_1.0.0   raster_2.5-2    sp_1.2-1

loaded via a namespace (and not attached):
 [1] magrittr_1.5       rgdal_1.1-1        rsconnect_0.4.1.11
parallel_3.2.2
 [5] tools_3.2.2        Rcpp_0.12.2        stringi_1.0-1
ncdf4_1.15
 [9] grid_3.2.2         lattice_0.20-33

  
  


  


    

    

      
      

        


  

        

      Michael Sumner
    

    

      
      #
    

  


  

      
On Wed, 25 May 2016 at 19:59 Mark R Payne <markpayneatwork at gmail.com> wrote:

        

            

      
      
      
    

        
I don't know the proper fix, but it works if you use a RasterLayer rather
than a 1-layer RasterBrick:

 crop(b[[1]],extent(320,340,55,60))

It's straightforwardly a bug IMO, due to the logic of a 1-layer degenerate
dimension - and one for the the raster maintainer.  If you can look at the
source of raster:::.readRowsBrickNetCDF you'll probably find it if for a
suggested fix or patch.

Cheers, Mike.