[R-sig-Geo] Raster: Can't read raster-created NCDF back in

[Michael Sumner]

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

> 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:
>
> > b <- brick("Bluefin_OISST_clim.nc")
> > plot(b)
>
> but these don't:
>
>

> > crop(b,extent(320,340,55,60))
> 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!
> > readAll(b)
> 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!
> >
>
>
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.


> 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
> ;
> }
>
> > sessionInfo()
> 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
> >
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
-- 
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

	[[alternative HTML version deleted]]