[R-sig-Geo] Error in netCDF file: cells are not equally spaced
Edzer Pebesma
edzer@pebe@m@ @end|ng |rom un|-muen@ter@de
Thu Aug 22 18:27:47 CEST 2019
Thanks!
stars::read_stars (GDAL) won't read coordinates or coordinate reference
system from this file, confirmed by running gdalinfo on it.
stars::read_ncdf (installing stars from github) reads it like this:
library(stars)
# Loading required package: abind
# Loading required package: sf
# Linking to GEOS 3.7.1, GDAL 2.4.2, PROJ 5.2.0
(r0 = read_ncdf("pr_Amon_GFDL-ESM2G_rcp45_r1i1p1_204101-204512.nc"))
# no 'var' specified, using pr
# other available variables:
# average_DT, average_T1, average_T2, lat, lon, bnds, time, time_bnds,
lat_bnds, lon_bnds
# No projection information found in nc file.
# Coordinate variable units found to be degrees,
# assuming WGS84 Lat/Lon.
# stars object with 3 dimensions and 1 attribute
# attribute(s):
# pr [kg/m^2/s]
# Min. :0.000e+00
# 1st Qu.:6.035e-06
# Median :1.700e-05
# Mean :2.799e-05
# 3rd Qu.:3.575e-05
# Max. :1.068e-03
# dimension(s):
# from to offset delta refsys point
# lon 1 144 NA NA +proj=longlat +datum=WGS8... NA
# lat 1 90 NA NA +proj=longlat +datum=WGS8... FALSE
# time 1 60 NA NA PCICt FALSE
# values
# lon [0,2.5),...,[357.5,360) [x]
# lat [-90,-88.98876),...,[88.98876,90) [y]
# time [2041-01-01,2041-02-01),...,[2045-12-01,2046-01-01)
st_get_dimension_values(r0, "lat") %>% diff %>% table
# .
# 1.51685393258427 2.02247191011234 2.02247191011236 2.02247191011237
# 2 1 74 12
where essentially only the first and last intervals differ from the rest.
So, yes, latitudes are not regular. Also, time is PCICt, and uses a 365
day calendar without leap years:
st_get_dimension_values(r0, "time") %>% diff %>% table
# .
# 28 30 31
# 5 20 34
On 8/22/19 8:55 PM, Frederico Faleiro wrote:
> Dear list,
> sorry about the file. The files from the CMIP5 are public, but need a
> registration first.
> You can access the same file in any of these links:
> https://www.sendspace.com/file/famgz3 and
> https://drive.google.com/file/d/1L1T03iQ6LU1yYRYeRypxwMB3E7fLjkJI/view?usp=sharing
> .
>
> In this meantime I will try the Mike' advices.
>
> Best regards,
>
> Frederico <https://www.sendspace.com/file/famgz3>
>
> On Thu, Aug 22, 2019 at 3:41 PM Frederico Faleiro <fvfaleiro using gmail.com>
> wrote:
>
>> Dear list,
>>
>> sorry about the file. The files from the CMIP5 are public, but need a
>> registration first. I am sending the same file attached.
>>
>> In this meantime I will try the Mike' advices.
>>
>> Best regards,
>>
>> Frederico
>>
>>
>>
>>
>>
>> On Thu, Aug 22, 2019 at 5:03 AM Michael Sumner <mdsumner using gmail.com> wrote:
>>
>>> raster does the wrong thing here in my opinion, the right outcome is to
>>> give a grid in index-extent (0, ncol, 0, nrow) and with no projection
>>> metadata (crs). There will be coordinate arrays in this file, and they need
>>> to be handled as though they are data (with local, x*y dependent values in
>>> every cell).
>>>
>>> With brick() you can proceed to process the data with no problem (using
>>> stopIfNotEqualSpaced = FALSE), but you can't rely on the extent being
>>> relevant, or any function that works in geographic space to do the right
>>> thing. To map a layer of these data you might try
>>>
>>> grd <- raster(yourfile, stopIfNotEqualSpaced = FALSE)
>>> coords <- brick(raster(yourfile, varname = "lon"),
>>> raster(yourfile, varname = "lat")) ## but only
>>> you will know the names of these variables values "lon", "lat" (might be
>>> "TLON", "TLAT" ? )
>>>
>>> In quadmesh there is a way to plot in geographic space that's not
>>> impossibly slow (possibly the coords will need a re-orientation transpose
>>> ...):
>>>
>>> quadmesh::mesh_plot(grd, coords = coords)
>>>
>>> It's likely that will "map" it properly, but CMIP is likely to wrap
>>> around an odd longitude (or something), and so cropping to your area and/or
>>> choosing a good project for the crs argument to mesh_plot is probably
>>> needed. You can investigate with plot(coords[[1]]) and plot(coords[[2]])
>>> to get a sense of their layout.
>>>
>>> In the stars package this is all internalized, with
>>> stars::read_stars(yourfile) catching all the information required as much
>>> as possible, and with plotting methods - but variable choice and actual
>>> results vary widely, and depend a lot on very specific details. (angstroms
>>> package has similar helpers for the approach but is unlikely to help here
>>> without access to the file to try)
>>>
>>> To help us guess further you can use the "ncdump" output, raster will
>>> print this with
>>>
>>> print(raster("yourfile"))
>>>
>>> and from that we could provide better guesses at variable names and
>>> subsetting etc.
>>>
>>> But, as Edzer asked, nothing is as good as having the file - any chance
>>> you can share it? (Personally I think the world rather needs more R
>>> attention on climate model output.)
>>>
>>> Cheers, Mike.
>>>
>>> On Thu, Aug 22, 2019 at 1:53 PM Frederico Faleiro <fvfaleiro using gmail.com>
>>> wrote:
>>>>
>>>> Dear list,
>>>>
>>>> I am using some netCDF files from the CMIP5 climate models in raster
>>>> package, but I am having some issues with one model.
>>>> The netCDF file from the GFDL-ESM2G model (e.g.
>>>>
>>> http://aims3.llnl.gov/thredds/fileServer/css03_data/cmip5/output1/NOAA-GFDL/GFDL-ESM2G/rcp45/mon/atmos/Amon/r1i1p1/v20120412/pr/pr_Amon_GFDL-ESM2G_rcp45_r1i1p1_204101-204512.nc
>>>> )
>>>> have the error message as in bellow example.
>>>>
>>>> #Example
>>>> library(raster)
>>>> s1 <- stack("pr_Amon_GFDL-ESM2G_rcp45_r1i1p1_204101-204512.nc")
>>>> Error in .rasterObjectFromCDF(x, type = objecttype, band = band, ...) :
>>>> cells are not equally spaced; you should extract values as points
>>>>
>>>> # I check some solutions that recomend force read the file with the
>>>> argument "stopIfNotEqualSpaced = F" as bellow.
>>>> sf <- stack("pr_Amon_GFDL-ESM2G_rcp45_r1i1p1_204101-204512.nc",
>>>> *stopIfNotEqualSpaced
>>>> = F*)
>>>> bf <- brick("pr_Amon_GFDL-ESM2G_rcp45_r1i1p1_204101-204512.nc",
>>>> *stopIfNotEqualSpaced
>>>> = F*)
>>>>
>>>> I performed some tests and only "works" with brick. However I did not
>>> find
>>>> any solution to check where is the problem and fix it in the file.
>>>>
>>>> How can I check if it is really an issue and fix it?
>>>>
>>>> sessionInfo()
>>>> R version 3.5.1 (2018-07-02)
>>>> Platform: x86_64-w64-mingw32/x64 (64-bit)
>>>> Running under: Windows 10 x64 (build 18362)
>>>>
>>>> Matrix products: default
>>>>
>>>> locale:
>>>> [1] LC_COLLATE=Portuguese_Brazil.1252 LC_CTYPE=Portuguese_Brazil.1252
>>>> [3] LC_MONETARY=Portuguese_Brazil.1252 LC_NUMERIC=C
>>>> [5] LC_TIME=Portuguese_Brazil.1252
>>>>
>>>> attached base packages:
>>>> [1] stats graphics grDevices utils datasets methods base
>>>>
>>>> other attached packages:
>>>> [1] raster_2.8-19 sp_1.3-1
>>>>
>>>> loaded via a namespace (and not attached):
>>>> [1] compiler_3.5.1 rgdal_1.4-3 Rcpp_1.0.1 codetools_0.2-15
>>>> ncdf4_1.16.1
>>>> [6] grid_3.5.1 lattice_0.20-35
>>>>
>>>> Best regards,
>>>>
>>>> Frederico
>>>>
>>>> [[alternative HTML version deleted]]
>>>>
>>>> _______________________________________________
>>>> R-sig-Geo mailing list
>>>> R-sig-Geo using r-project.org
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>>
>>>
>>> --
>>> Michael Sumner
>>> Software and Database Engineer
>>> Australian Antarctic Division
>>> Hobart, Australia
>>> e-mail: mdsumner using gmail.com
>>>
>>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081
-------------- next part --------------
A non-text attachment was scrubbed...
Name: pEpkey.asc
Type: application/pgp-keys
Size: 2472 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20190822/3957cc46/attachment.bin>
More information about the R-sig-Geo
mailing list