Open anthoni-p opened 3 years ago
Hi Peter,
It looks like TRENDY data is now just slightly over R's integer size limit. That is bad luck. Unfortunately I haven't implemented reading Fields per layer. I never quite had the reason, but now there is one.. It shouldn't be too hard but will take a bit of time.
Is this urgent for you? I am on holiday right now.
Hi Peter,
I have implemented a "layer" argument to getField() which allows selection of individual variables inside the netCDF file. I think this is a useful thing. But I just realised that doesn't solve the problem you had above, because it doesn't allow selecting individual PFT indices.
So, I propose adding an additional argument, something like 'layer.dim.indices' which would allow subsetting on the PFT (or whatever) dimension. Will take a bit more time to implement however...
@anthoni-p I guess this issue hasn't gone away? I never implemented the 'layer.dim.index' idea. I can see a way to solve this though, basically read the netCDF is separate slices for each element in the layer axis (in this case your 26 PFTs). That will keep the size below the maximum integers and keep the memory footprint down. It might even make the code faster, but I am not sure....
If you are still in interested, can you supply me with a test file?
getting an error when trying to read a netcdf with per PFT data, not sure where it comes from.
Can only parts of such a file get read in, e.g. PFT index 1, how can we specify that in a getfield?
here ncdump info of the files, yearly data, so not really that many time steps: netcdf LPJ-GUESS_S3_cVegpft { dimensions: PFT = 26 ; longitude = 720 ; latitude = 360 ; time = UNLIMITED ; // (321 currently) variables: int PFT(PFT) ; PFT:long_name = "plant functional type" ; double longitude(longitude) ; longitude:units = "degrees_east" ; longitude:long_name = "longitude" ; double latitude(latitude) ; latitude:units = "degrees_north" ; latitude:long_name = "latitude" ; double time(time) ; time:units = "days since 1700-01-01" ; time:long_name = "time" ; time:calendar = "noleap" ; float cVegpft(time, latitude, longitude, PFT) ; cVegpft:units = "kg m-2" ; cVegpft:_FillValue = -99999.f ; cVegpft:long_name = "Vegtype level Carbon in Vegetation" ;