0

I am trying to create a netCDF file from an excel sheet that contains precipitation data. The excel file is organized with latitude and longitude in columns: a column with latitudes, a column with longitudes both in EPSG:3006 - SWEREF99 TM. The subsequent columns contain the data for each day and I am reconstructing this into a netCDF file. The code works, but when opening the netCDF file the projection is completely wrong. I have followed the answer given in this post, but without getting the projection right (both when opening it in QGIS and running gdalinfo on the file).


nc_file = os.path.join(outpath, 'gridded_precip.nc')
nc = Dataset(nc_file, 'w', format='NETCDF4')

calendar = 'gregorian'
units = 'days since '+ str(start_date)

tim_dim = nc.createDimension('time', len(time))
lat_dim = nc.createDimension('lat', 9)
long_dim = nc.createDimension('lon', 11) 

time_var = nc.createVariable(varname='time', dimensions=('time',), datatype='float64')
time_var[:] = netCDF4.date2num(times, units=units, calendar=calendar)

latitude_var = nc.createVariable(varname='lat', dimensions=('lat',),datatype='f4')
latitude_var.units = 'meters'
latitude_var.axis = 'Y'
latitude_var[:] = lat.mean(axis=1)

longitude_var = nc.createVariable(varname='lon', dimensions=('lon',), datatype='f4')
longitude_var.units = 'meters'
longitude_var.axis = 'X'
longitude_var[:] = long.mean(axis=0)

precip_var = nc.createVariable('precipitation_mm_day','f4', ('time', 'lat', 'lon'),)
precip_var[:] = data
precip_var.units = 'mm/day'
precip_var.setncattr('grid_mapping', 'spatial_ref')

# Add crs to dataset
crs = nc.createVariable('spatial_ref', 'i4')
project_crs = 'PROJCS["SWEREF99 TM",GEOGCS["SWEREF99",DATUM["SWEREF99",SPHEROID["GRS 1980",6378137,298.257222101],TOWGS84[0,0,0,0,0,0,0]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4619"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",15],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","3006"]]'
crs.spatial_ref = project_crs
nc.close() 

Any ideas of what might be wrong?
A second but less important question: Both the longitude and latitude data are in 2D matrices, but I could only add the longitude and latitude data as vectors (the reason why a I am taking the average in the code)). Is there a way to have 2D latitude and longitude data in a netCDF?

1
  • Create your data array on a grid with dimensions (time,y,x) and the lon and lat arrays on grids with dimensions (y,x). Add a coordinates attribute to the data array naming the lon and lat arrays. See section 5.2 of the CF Convention for an example (which has a level dimension rather than a time dim), cfconventions.org/Data/cf-conventions/cf-conventions-1.11/…
    – R. Schmunk
    Commented Mar 15 at 21:34

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.