Thursday, 9 January 2020

raster - Set projection (netCDF) with R-ncdf4-Package


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



No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...