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 <btup...@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 <elgab...@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 <btup...@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 <elgab...@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@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@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@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]] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo