1 | import netCDF4 as nc |
---|
2 | import sys |
---|
3 | import math |
---|
4 | import numpy as np |
---|
5 | |
---|
6 | def from_reduced(N,M): |
---|
7 | "N elements from south to north and N elements around equator " |
---|
8 | hmax = 2*math.pi/N |
---|
9 | hmin = hmax/2 |
---|
10 | nlon = N |
---|
11 | cells_lon = [] |
---|
12 | cells_lat = [] |
---|
13 | for i in range(M/2): |
---|
14 | lat1 = 180.0/M*i |
---|
15 | lat2 = 180.0/M*(i+1) |
---|
16 | y = math.sin(lat1*math.pi/180) |
---|
17 | r = math.cos(lat1*math.pi/180) |
---|
18 | h = 2.0*r/nlon |
---|
19 | reduce_nlon = (h < hmin) and (i > 0) and (nlon > 4) |
---|
20 | if reduce_nlon: |
---|
21 | nlon = nlon/2 |
---|
22 | for j in range(nlon): |
---|
23 | lon1 = 360.0*j/nlon |
---|
24 | lon2 = 360.0*(j+1)/nlon |
---|
25 | bounds_lon = [lon1, lon1, lon2, lon2] |
---|
26 | bounds_lat = [lat1, lat2, lat2, lat1] |
---|
27 | if reduce_nlon: |
---|
28 | bounds_lon.append((lon1+lon2)/2) |
---|
29 | bounds_lat.append(lat1) |
---|
30 | else: # close by netCDF convention |
---|
31 | bounds_lon.append(bounds_lon[0]) |
---|
32 | bounds_lat.append(bounds_lat[0]) |
---|
33 | # northern hemisphere |
---|
34 | cells_lon.append(bounds_lon) |
---|
35 | cells_lat.append(bounds_lat) |
---|
36 | # southern hemisphere |
---|
37 | cells_lon.append(bounds_lon) |
---|
38 | cells_lat.append(list(-np.array(bounds_lat))) # convert to array to negate elementwise |
---|
39 | return np.array(cells_lon), np.array(cells_lat) |
---|
40 | |
---|
41 | for N in [64, 128, 256, 512]: |
---|
42 | filename = "reduced" + str(N) + ".nc" |
---|
43 | |
---|
44 | print "Generating: N =", N |
---|
45 | lon, lat = from_reduced(N*2,N) |
---|
46 | |
---|
47 | print lon.shape[0], "cells -> writing as ", filename |
---|
48 | |
---|
49 | f = nc.Dataset(filename,'w') |
---|
50 | f.createDimension('n_vert', 5) |
---|
51 | f.createDimension('n_cell', lon.shape[0]) |
---|
52 | |
---|
53 | var = f.createVariable('lat', 'd', ('n_cell')) |
---|
54 | var.setncattr("long_name", "latitude") |
---|
55 | var.setncattr("units", "degrees_north") |
---|
56 | var.setncattr("bounds", "bounds_lat") |
---|
57 | var[:] = np.zeros(lon.shape[0]) |
---|
58 | var = f.createVariable('lon', 'd', ('n_cell')) |
---|
59 | var.setncattr("long_name", "longitude") |
---|
60 | var.setncattr("units", "degrees_east") |
---|
61 | var.setncattr("bounds", "bounds_lon") |
---|
62 | var[:] = np.zeros(lon.shape[0]) |
---|
63 | |
---|
64 | var = f.createVariable('bounds_lon', 'd', ('n_cell','n_vert')) |
---|
65 | var[:] = lon |
---|
66 | var = f.createVariable('bounds_lat', 'd', ('n_cell','n_vert')) |
---|
67 | var[:] = lat |
---|
68 | var = f.createVariable('val', 'd', ('n_cell')) |
---|
69 | var.setncattr("coordinates", "lon lat") |
---|
70 | var[:] = np.arange(lon.shape[0]) |
---|
71 | f.close() |
---|