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
| import os, math, time, sys, glob, logging
import xarray as xr
import numpy as np
from scipy import stats as st
logging.basicConfig(level=logging.WARNING)
in_static = xr.open_dataset("in_static.nc")
lcz_dataset = xr.open_dataset("lcz.nc")
print(in_static)
print(lcz_dataset)
latCell = in_static["latCell"][:].values*57.2957795
lonCell = in_static["lonCell"][:].values*57.2957795
nCells = latCell.shape[0]
dcEdge = in_static["dcEdge"][:].values
edgesOnCell = in_static["edgesOnCell"][:].values
nEdges = dcEdge.shape[0]
maxEdges = edgesOnCell.shape[1]
ivgtyp = in_static["ivgtyp"][:].values
print("="*40)
print("nCells, nEdges, maxEdges = ", nCells, nEdges, maxEdges)
print("latCell = ", latCell.shape)
print("lonCell = ", lonCell.shape)
print("dcEdge = ", dcEdge.shape)
print("edgesOnCell = ", edgesOnCell.shape)
print("ivgtyp = ", ivgtyp.shape)
print("min/max = ", np.min(latCell), np.max(latCell), np.min(lonCell), np.max(lonCell))
lat = lcz_dataset["lat"][:].values
lon = lcz_dataset["lon"][:].values
Band2 = lcz_dataset["Band2"][:].values
nlat = lat.shape[0]
nlon = lon.shape[0]
min_lat, max_lat = np.min(lat), np.max(lat)
min_lon, max_lon = np.min(lon), np.max(lon)
dlatlon = np.abs(lat[1] - lat[0])*111
print("="*40)
print("nlat, nlon = ", nlat, nlon)
print("lat = ", lat.shape)
print("lon = ", lon.shape)
print("min_lat, max_lat = ", min_lat, max_lat)
print("min_lon, max_lon = ", min_lon, max_lon)
print("diff lat/lon = ", lat[1:5]-lat[0:4], lon[1:5]-lon[0:4])
print("Band2 = ", Band2.shape)
print("="*40)
tmp_count_urban = 0
tmp_count_urban_change = 0
for px in range(nCells):
if(px%5000 == 0):
print(round(px/nCells*100), " %")
if(ivgtyp[px] == 13 or ivgtyp[px] == 32):
tmp_count_urban = tmp_count_urban + 1
if( (min_lat <= latCell[px]) and (max_lat >= latCell[px]) and (min_lon <= lonCell[px]) and (max_lon >= lonCell[px]) ):
index_lat = np.argmin(np.abs(lat - latCell[px]))
index_lon = np.argmin(np.abs(lon - lonCell[px]))
tmp_info = str(latCell[px] - lat[index_lat]) +" "+ str(lonCell[px] - lon[index_lon]) +" "+ str(Band2[index_lat, index_lon])
logging.info(tmp_info)
tmp_count_urban_change = tmp_count_urban_change +1
if ( Band2[index_lat, index_lon] < 11 or Band2[index_lat, index_lon] == 15 ):
tmp_dx = np.min(dcEdge[edgesOnCell[px,:]])/1000
tmp_num_cell = round((tmp_dx/dlatlon)/2)
tmp_ar1 = Band2[index_lat-tmp_num_cell:index_lat+tmp_num_cell, index_lon-tmp_num_cell:index_lon+tmp_num_cell].flatten()
tmp_ar2 = tmp_ar1
for py in [11,12,13,14,16,17]:
if np.where(tmp_ar2 == py)[0].size != 0:
tmp_ar2 = np.delete(tmp_ar2, np.where(tmp_ar2 == py)[0][:])
tmp_array = tmp_ar2 + 30
tmp = st.mode(tmp_array)[0]
if (len(tmp) == 0):
continue
if (int(tmp[0]) == 45):
in_static["ivgtyp"][px] = 41
else:
in_static["ivgtyp"][px] = int(tmp[0])
else:
in_static["ivgtyp"][px] = 32
else:
in_static["ivgtyp"][px] = 32
print("Total urban in MPAS mesh = ", tmp_count_urban)
print("Change urban in MPAS mesh = ", tmp_count_urban_change)
in_static.to_netcdf("out_static.nc")
|