| 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() |
|---|
