[R-sig-Geo] converting raster to netcdf .nc with keeping CRS information
Michael Sumner
md@umner @end|ng |rom gm@||@com
Wed Feb 13 06:11:46 CET 2019
I agree, I wouldn't try to bend NetCDF to work this way via raster. The
convention/s require parameters of the projection to be specified via
attributes and you could do that via RNetCDF or ncdf4, but it's a lot of
painful work.
It's not easy to create NetCDF from R, it's a simple as that - and
typically those files don't play well with projections - it's changed
recently but you'll see a lot of patchy partial support across different
softwares.
You can use rgdal or call out out to GDAL - but then round-tripping gives
warnings, because again raster has to unpack all the components of metadata
that GDAL created, and none of this is very systematic and clearly is
broken atm.
system('gdal_translate r001.nc r001_aux.nc -of NetCDF -a_srs
"+init=epsg:3412"')
## you could do this with stars too
rgdal::writeGDAL(as(r, "SpatialGridDataFrame"), "r001_rgdal.nc", drivername
= "NetCDF")
But, beware as in neither case is raster guaranteed to get the right
interpretation (THIS IS WRONG):
raster("r001_aux.nc")
raster("r001_rgdal.nc")
class : RasterLayer
dimensions : 1328, 1264, 1678592 (nrow, ncol, ncell)
resolution : 6250, 6250 (x, y)
extent : -3950000, 3950000, -3950000, 4350000 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=stere +lon_0=0 +x_0=0 +y_0=0 +lat_0=-90 +pm=0
+a=6378273 +rf=298.279411123064
data source : r001_rgdal.nc
names : GDAL.Band.Number.1
zvar : Band1
BTW, you called your X and Y "longitude" and "latitude" but this grid isn't
defined by those (eastings and northings is not a better name either, but
that's another problem).
Do you have to create new files? I find that constantly a fraught problem,
For usage it's way easier to wrap stuff up with R code and work around all
these problems.
Cheers, Mike.
On Wed, 13 Feb 2019 at 12:47 Ben Tupper <btupper using bigelow.org> wrote:
> Hi Ahmed,
>
> I can confirm that on my system I get the same result as you.
>
> Does your workflow confine you to using NetCDF? I can successfully
> write/read to the native raster format.
>
> r2 <- writeRaster(r, filename = "r001.grd", overwrite=TRUE)
> r2
> #class : RasterLayer
> #dimensions : 1328, 1264, 1678592 (nrow, ncol, ncell)
> #resolution : 6250, 6250 (x, y)
> #extent : -3950000, 3950000, -3950000, 4350000 (xmin, xmax, ymin,
> ymax)
> #coord. ref. : +proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0 +k=1 +x_0=0
> +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs
> #data source : /home/btupper/r001.grd
> #names : layer
> #values : 0, 100 (min, max)
>
>
> r3 <- raster("r001.grd")
> r3
> #class : RasterLayer
> #dimensions : 1328, 1264, 1678592 (nrow, ncol, ncell)
> #resolution : 6250, 6250 (x, y)
> #extent : -3950000, 3950000, -3950000, 4350000 (xmin, xmax, ymin,
> ymax)
> #coord. ref. : +proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0 +k=1 +x_0=0
> +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs
> #data source : /home/btupper/r001.grd
> #names : layer
> #values : 0, 100 (min, max)
>
> Ben
>
> > sessionInfo()
> R version 3.5.1 (2018-07-02)
> Platform: x86_64-redhat-linux-gnu (64-bit)
> Running under: CentOS Linux 7 (Core)
>
> Matrix products: default
> BLAS/LAPACK: /usr/lib64/R/lib/libRblas.so
>
> locale:
> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
> [9] LC_ADDRESS=C LC_TELEPHONE=C
> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
>
> attached base packages:
> [1] stats graphics grDevices utils datasets methods base
>
> other attached packages:
> [1] raster_2.8-19 sp_1.3-1 ncdf4_1.16
>
> loaded via a namespace (and not attached):
> [1] compiler_3.5.1 rgdal_1.3-6 Rcpp_1.0.0 codetools_0.2-15
> [5] grid_3.5.1 lattice_0.20-35
>
>
>
>
>
> > On Feb 12, 2019, at 9:19 AM, Ahmed El-Gabbas <elgabass using gmail.com> wrote:
> >
> > Thanks Ben for your suggestion.
> >
> > Adding *prj = TRUE *argument caused this error:
> > Error in ncdf4::ncvar_def(varname, varunit, list(xdim, ydim),
> > NAvalue(x), :
> > unused argument (prj = T)
> >
> > Ahmed
> >
> > On Tue, Feb 12, 2019 at 3:13 PM Ben Tupper <btupper using bigelow.org> wrote:
> >
> >> Hi,
> >>
> >> Have you tried using the 'prj' argument to writeRaster? I don't know
> that
> >> it will be the solution, but according to ...
> >>
> >>
> >>
> https://www.rdocumentation.org/packages/raster/versions/2.8-19/topics/writeRaster
> >>
> >> It is documented among the '...' arguments. Setting prj = TRUE will
> cause
> >> the crs to be written to the file.
> >>
> >> Cheers,
> >> Ben
> >>
> >> On Feb 12, 2019, at 8:34 AM, Ahmed El-Gabbas <elgabass using gmail.com>
> wrote:
> >>
> >> Hello all,
> >>
> >> I am having a problem while converting a raster object as NetCDF (.nc)
> >> file, with keeping the CRS information in the output file.
> >> Here is a reproducible code:
> >>
> >> require(raster)
> >> require(ncdf4)
> >> CurrTemp <- tempfile()
> >> download.file(url = "
> https://seaice.uni-bremen.de/data/amsre/asi_daygrid_swath/s6250/2003/feb/Antarctic/asi-s6250-20030214-v5.hdf",
> destfile = CurrTemp, mode = "wb", quiet = T)
> >> r <- raster(CurrTemp)
> >> r <- flip(r,2)
> >> extent(r) <- c(-3950000, 3950000, -3950000, 4350000)
> >> crs(r) <- "+proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0 +k=1 +x_0=0
> +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs"
> >> r# class : RasterLayer # dimensions : 1328, 1264, 1678592
> (nrow, ncol, ncell)# resolution : 6250, 6250 (x, y)# extent :
> -3950000, 3950000, -3950000, 4350000 (xmin, xmax, ymin, ymax)# coord. ref.
> : +proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378273
> +b=6356889.449 +units=m +no_defs # data source : in memory# names :
> layer # values : 0, 100 (min, max)
> >>
> >> So far the raster object reads well, with CRS information included.
> >> However, when I try to save it as .nc file, R prints "coord. ref. : NA",
> >> and the produced file does not contain the CRS information.
> >>
> >> writeRaster(r, filename = "O:/Ahmed/r001.nc", varname="IceConc",
> >> overwrite=TRUE, format="CDF",
> >> xname="Longitude", yname="Latitude")# class :
> RasterLayer # dimensions : 1328, 1264, 1678592 (nrow, ncol, ncell)#
> resolution : 6250, 6250 (x, y)# extent : -3950000, 3950000,
> -3950000, 4350000 (xmin, xmax, ymin, ymax)# coord. ref. : NA # data source
> : O:/Ahmed/r001.nc # names : IceConc # zvar : IceConc
> >>
> >> raster("O:/Ahmed/r001.nc")# class : RasterLayer # dimensions :
> 1328, 1264, 1678592 (nrow, ncol, ncell)# resolution : 6250, 6250 (x, y)#
> extent : -3950000, 3950000, -3950000, 4350000 (xmin, xmax, ymin,
> ymax)# coord. ref. : NA # data source : O:/Ahmed/r001.nc # names :
> IceConc # zvar : IceConc
> >>
> >> Any solution?
> >>
> >> N.B. I also sent the same question at stackoverflow
> >> <
> https://stackoverflow.com/questions/54593552/saving-r-raster-to-netcdf-nc-file-with-keeping-crs-information
> >,
> >> without a solution so far.
> >>
> >> Thanks
> >>
> >> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
> >> *Dr. Ahmed El-Gabbas,*
> >> Ocean Acoustics Lab,
> >> Alfred-Wegener-Institut
> >> Helmholtz-Zentrum für Polar und Meeresforschung
> >> <image.png>
> >> My Website: https://elgabbas.github.io
> >> _______________________________________________
> >> R-sig-Geo mailing list
> >> R-sig-Geo using r-project.org
> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >>
> >>
> >> Ben Tupper
> >> Bigelow Laboratory for Ocean Sciences
> >> 60 Bigelow Drive, P.O. Box 380
> >> East Boothbay, Maine 04544
> >> http://www.bigelow.org
> >>
> >> Ecological Forecasting: https://eco.bigelow.org/
> >>
> >>
> >>
> >>
> >>
> >>
> >
> > [[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
>
> Ben Tupper
> Bigelow Laboratory for Ocean Sciences
> 60 Bigelow Drive, P.O. Box 380
> East Boothbay, Maine 04544
> http://www.bigelow.org
>
> Ecological Forecasting: https://eco.bigelow.org/
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo using 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]]
More information about the R-sig-Geo
mailing list