1

I want to create a netCDF-file that has a defined projection ("WGS 84 / UTM zone 32N") with the netcdf4 package in R. The input float array has the size of x=880 and y=760. So i define the two dimensions:

ext = extent(raster_stack.ras)
res = 50
x = seq(ext[1]+res/2,ext[2]-res/2,50)
y = seq(ext[3]+res/2,ext[4]-res/2,50)

dimx = ncdim_def(name = 'x', longname = 'x coordinate of projection', units='m', vals = as.double(x))
dimy = ncdim_def(name = 'y', longname = 'y coordinate of projection', units='m', vals = as.double(y))

The extent gives the xmin, xmax, ymin and ymax of the input raster and define the columns and rows of the array. To get the midpoint of one pixel it's necessary to add a half of the pixelsize (final pixelsize should 50 m). Finally with dimx and dimy the dimensions of the netCDF are getting produced.

Then I define the variable as follows:

varz = ncvar_def(name = out.layername, longname = out.layername, units = 'percent', dim = list(dimx,dimy), missval = -99, prec = 'float')

outnc = nc_create(out.name,varz,force_v4=TRUE)    
ncvar_put(outnc,varz,soil_moisture)
nc_close(outnc)    

If I create the netCDF there is no projection for the file. But i can't find any documentation how to define the projection or set it in a variable.

To have an example for the netCDF Header, I converted a geoTIFF with gdal_translate to netCDF and get the following header. There you could see the defined projection. How could I declare it in R with netCDF4 ?

 2 variables (excluding dimension variables):
    char transverse_mercator[]   
        grid_mapping_name: transverse_mercator
        longitude_of_central_meridian: 9
        false_easting: 5e+05
        false_northing: 0
        latitude_of_projection_origin: 0
        scale_factor_at_central_meridian: 0.9996
        longitude_of_prime_meridian: 0
        semi_major_axis: 6378137
        inverse_flattening: 298.257223563
        spatial_ref: PROJCS["WGS 84 / UTM zone 32N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",9],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32632"]]
        GeoTransform: 550000 50 0 5740000 0 -50 
    float Band1[x,y]   
        long_name: GDAL Band Number 1
        _FillValue: -99
        grid_mapping: transverse_mercator
 2 dimensions:
    x  Size:880
        standard_name: projection_x_coordinate
        long_name: x coordinate of projection
        units: m
    y  Size:760
        standard_name: projection_y_coordinate
        long_name: y coordinate of projection
        units: m
1

Given that start with RasterStack you might consider writeRaster. Example data:

library(raster)
r <- raster(crs="+proj=utm +zone=32 +datum=WGS84", xmn=0, xmx=10, ymn=0, ymx=10)
values(r) <- 1:ncell(r)
s <- stack(r, r/2)

This will write your RasterStack to a ncdf file:

x <- writeRaster(s, filename="test.nc") 

It includes the CRS info. But not in the CF-1.4 standard way, but rather in the more useful proj.4 notation. See

print(x)

But perhaps that does not matter to you.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.