-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathcreate_grids.py
135 lines (114 loc) · 5.28 KB
/
create_grids.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
#!/usr/bin/env python
import xarray as xr
import numpy as np
#----------------------- bedmachine ---------------------------------
PROJSTRING= "+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
def add_lon_lat(ds, PROJSTRING, x='x', y='y', chunks={}):
""" add longitude and latitude as compute from the inverse projection
given in PROJSTRING
PARAMETERS:
-----------
ds: xarray.Dataset
PROJSTRING: str
"""
from pyproj import CRS, Transformer
# create the coordinate reference system
crs = CRS.from_proj4(PROJSTRING)
# create the projection from lon/lat to x/y
proj = Transformer.from_crs(crs.geodetic_crs, crs)
# make x,y 2d arrays
xx, yy = np.meshgrid(ds[x].values, ds[y].values)
# compute the lon/lat
lon, lat = proj.transform(xx, yy, direction='INVERSE')
# add to dataset
ds['lon'] = xr.DataArray(data=lon, dims=('y', 'x'))
ds['lat'] = xr.DataArray(data=lat, dims=('y', 'x'))
ds['lon'].attrs = dict(units='degrees_east')
ds['lat'].attrs = dict(units='degrees_north')
return ds
def create_bedmachine_xy():
"""recreate coordinates from bedmachine into xarray Dataset"""
# x = np.arange(-3333000,3333000+500,500, dtype=np.int32)
# y = np.arange(3333000,-3333000-500,-500, dtype=np.int32)
ds = xr.Dataset()
ds['x'] = xr.DataArray(data=x, dims=('x'))
ds['y'] = xr.DataArray(data=y, dims=('y'))
ds['x'].attrs = {'long_name': "Cartesian x-coordinate",
'standard_name': "projection_x_coordinate",
'units': "meter"}
ds['y'].attrs = {'long_name': "Cartesian y-coordinate",
'standard_name': "projection_y_coordinate",
'units': "meter"}
return ds
#----------------------- gebco --------------------------------------
# quite interestingly this function is unable to reproduce bitwise
# identical values from the original gebco grid, with differences O(10e-14)
def create_grid_gebco():
""" GEBCO has a 15" of arc regular lat/lon grid"""
increment = np.float64(1) / np.float64(240)
half = np.float64(1) / np.float64(2)
nx=360*240
ny=180*240
# longitude
lonmin = np.float64(-180.)
lonmax = np.float64(180.)
#lon_center = np.empty((nx), dtype=np.float64)
#lon_center[0] = lonmin + half * increment
#for k in range(1,nx):
# lon_center[k] = lon_center[0] + (k * increment)
lon_edges = np.empty((nx+1), dtype=np.float64)
lon_edges[0] = lonmin
for k in range(1,nx+1):
lon_edges[k] = lon_edges[0] + (k * increment)
lon_center = half * lon_edges[1:] + half * lon_edges[:-1]
# latitude
latmin = np.float64(-90.)
latmax = np.float64(90.)
lat_center = np.empty((ny), dtype=np.float64)
lat_center[0] = latmin + half * increment
for k in range(1,ny):
lat_center[k] = lat_center[0] + (k * increment)
ds = xr.Dataset()
ds['lon'] = xr.DataArray(data=lon_center, dims=('lon'))
ds['lat'] = xr.DataArray(data=lat_center, dims=('lat'))
ds['lon'].attrs = {'standard_name': "longitude",
'long_name': "longitude",
'units': "degrees_east"}
ds['lat'].attrs = {'standard_name': "latitude",
'long_name': "latitude",
'units': "degrees_north"}
return ds
def create_grid_gebco_arctic(ds, fileout):
""" subset gebco grid where bedmachine exists (and Arctic6 grid exists) """
# gebco_SO = ds.sel(lat=slice(58.5, 84.17),lon=slice(-90.9, 9.27))
gebco_SO = ds.sel(lat=slice(58.5, 84.17),lon=slice(-76.7, -4.6))
gebco_SO['lon'].attrs = dict(units='degrees_east')
gebco_SO['lat'].attrs = dict(units='degrees_north')
gebco_encoding = {'lon': {'dtype': 'float64'},
'lat': {'dtype': 'float64'}}
gebco_SO.to_netcdf(fileout, format='NETCDF3_64BIT',
engine='netcdf4', encoding=gebco_encoding)
#----------------------- main ---------------------------------------
if __name__ == '__main__':
# create gebco grid: not working yet
#gebco = create_grid_gebco()
#gebco_encoding = {'lon': {'dtype': 'float64'},
# 'lat': {'dtype': 'float64'}}
#gebco.to_netcdf('grid_gebco_30sec.nc', format='NETCDF3_64BIT',
# engine='netcdf4', encoding=gebco_encoding)
# subset gebco grid:
# gebco30 = xr.open_dataset('grid_gebco_30sec_original.nc')
# create_grid_gebco_antarctic(gebco30, 'grid_gebco_30sec_southof62.nc')
# Don't need to redo this.
gebco15 = xr.open_dataset('/import/AKWATERS/jrcermakiii/bathy/gebco/GEBCO_2020.nc')
create_grid_gebco_arctic(gebco15, 'grid_gebco_15sec_northof58.nc')
# create bedmachine grid:
# bedmachine = create_bedmachine_xy()
bedmachine = xr.open_dataset('/import/AKWATERS/kshedstrom/bathy/BedMachineGreenland-2017-09-20.nc')
bedmachine = add_lon_lat(bedmachine, PROJSTRING)
bedmachine_encoding = {'lon': {'dtype': 'float64', '_FillValue': 1.0e+15},
'lat': {'dtype': 'float64', '_FillValue': 1.0e+15},
'x': {'dtype': 'int32'},
'y': {'dtype': 'int32'}}
bedmachine.to_netcdf('grid_bedmachineGrn.nc', format='NETCDF3_64BIT',
engine='netcdf4', encoding=bedmachine_encoding)