[3944] | 1 | #!/usr/bin/env python3 |
---|
| 2 | # -*- coding: utf-8 -*- |
---|
| 3 | #--------------------------------------------------------------------------------# |
---|
| 4 | # This file is part of the PALM model system. |
---|
| 5 | # |
---|
| 6 | # PALM is free software: you can redistribute it and/or modify it under the terms |
---|
| 7 | # of the GNU General Public License as published by the Free Software Foundation, |
---|
| 8 | # either version 3 of the License, or (at your option) any later version. |
---|
| 9 | # |
---|
| 10 | # PALM is distributed in the hope that it will be useful, but WITHOUT ANY |
---|
| 11 | # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR |
---|
| 12 | # A PARTICULAR PURPOSE. See the GNU General Public License for more details. |
---|
| 13 | # |
---|
| 14 | # You should have received a copy of the GNU General Public License along with |
---|
| 15 | # PALM. If not, see <http://www.gnu.org/licenses/>. |
---|
| 16 | # |
---|
[4481] | 17 | # Copyright 1997-2020 Leibniz Universitaet Hannover |
---|
[3944] | 18 | #--------------------------------------------------------------------------------# |
---|
| 19 | # |
---|
| 20 | # Current revisions: |
---|
| 21 | # ----------------- |
---|
| 22 | # |
---|
[4023] | 23 | # |
---|
[3944] | 24 | # Former revisions: |
---|
| 25 | # ----------------- |
---|
| 26 | # $Id: palm_csd_tools.py 3773 2019-03-01 08:56:57Z maronga $ |
---|
[4250] | 27 | # Bugfix: no bridges were generated in buildings_3d at the east and north boundary |
---|
| 28 | # |
---|
| 29 | # 3773 2019-03-01 08:56:57Z maronga |
---|
[4023] | 30 | # Bugfix: missing conversion to integer |
---|
| 31 | # |
---|
| 32 | # |
---|
| 33 | # 3773 2019-03-01 08:56:57Z maronga |
---|
[4021] | 34 | # Unspecified changes |
---|
[3944] | 35 | # |
---|
| 36 | # 3726 2019-02-07 18:22:49Z maronga |
---|
| 37 | # Index bound bugfix |
---|
| 38 | # |
---|
| 39 | # 3668 2019-01-14 12:49:24Z maronga |
---|
| 40 | # Minor changes |
---|
| 41 | # |
---|
| 42 | # 3567 2018-11-27 13:59:21Z maronga |
---|
| 43 | # Initial revisions |
---|
| 44 | # |
---|
| 45 | # |
---|
| 46 | # Description: |
---|
| 47 | # ------------ |
---|
| 48 | # Support routines for palm_csd |
---|
| 49 | # |
---|
| 50 | # @Author Bjoern Maronga (maronga@muk.uni-hannover.de) |
---|
| 51 | # |
---|
| 52 | #------------------------------------------------------------------------------# |
---|
| 53 | import numpy as np |
---|
| 54 | from scipy.interpolate import interp2d |
---|
| 55 | |
---|
| 56 | def blend_array_2d(array1,array2,radius): |
---|
| 57 | # Blend over the parent and child terrain height within a radius of 50 px |
---|
| 58 | |
---|
| 59 | gradient_matrix = np.copy(array1) |
---|
| 60 | gradient_matrix[:,:] = 1.0 |
---|
| 61 | gradient_px = 50 |
---|
| 62 | |
---|
| 63 | for i in range(0,gradient_px): |
---|
| 64 | gradient_matrix[:,i] = float(i)/float(gradient_px) |
---|
| 65 | |
---|
| 66 | |
---|
| 67 | for i in range(len(gradient_matrix[0,:])-gradient_px,len(gradient_matrix[0,:])): |
---|
| 68 | gradient_matrix[:,i] = float(len(gradient_matrix[0,:])-i)/float(gradient_px) |
---|
| 69 | |
---|
| 70 | |
---|
| 71 | for j in range(0,gradient_px): |
---|
| 72 | for i in range(0,len(gradient_matrix[0,:])): |
---|
| 73 | if gradient_matrix[j,i] == 1.0: |
---|
| 74 | gradient_matrix[j,i] = float(j)/float(gradient_px) |
---|
| 75 | else: |
---|
| 76 | gradient_matrix[j,i] = (gradient_matrix[j,i] + float(j)/float(gradient_px))/2.0 |
---|
| 77 | |
---|
| 78 | |
---|
| 79 | for j in range(len(gradient_matrix[:,0])-gradient_px,len(gradient_matrix[:,0])): |
---|
| 80 | for i in range(0,len(gradient_matrix[0,:])): |
---|
| 81 | if gradient_matrix[j,i] == 1.0: |
---|
| 82 | gradient_matrix[j,i] = (len(gradient_matrix[:,0])-j)/float(gradient_px) |
---|
| 83 | else: |
---|
| 84 | gradient_matrix[j,i] = (gradient_matrix[j,i] + (len(gradient_matrix[:,0])-j)/float(gradient_px))/2.0 |
---|
| 85 | |
---|
| 86 | array_blended = array1 * gradient_matrix + (1.0 - gradient_matrix ) * array2 |
---|
| 87 | |
---|
| 88 | return array_blended |
---|
| 89 | |
---|
| 90 | |
---|
| 91 | def interpolate_2d(array,x1,y1,x2,y2): |
---|
| 92 | tmp_int2d = interp2d(x1,y1,array,kind='linear') |
---|
| 93 | array_ip = tmp_int2d(x2.astype(float), y2.astype(float)) |
---|
| 94 | |
---|
| 95 | return array_ip |
---|
| 96 | |
---|
| 97 | |
---|
| 98 | def bring_to_palm_grid(array,x,y,dz): |
---|
| 99 | |
---|
| 100 | # Bring the parent terrain height to the child grid |
---|
| 101 | k_tmp = np.arange(0,max(array.flatten())+dz*2,dz) |
---|
| 102 | k_tmp[1:] = k_tmp[1:] - dz * 0.5 |
---|
| 103 | for l in range(0,len(x)): |
---|
| 104 | for m in range(0,len(y)): |
---|
| 105 | for k in range(1,len(k_tmp+1)): |
---|
| 106 | if k_tmp[k] > array[m,l]: |
---|
| 107 | array[m,l] = k_tmp[k-1] + dz * 0.5 |
---|
| 108 | break |
---|
| 109 | |
---|
| 110 | return array |
---|
| 111 | |
---|
| 112 | |
---|
| 113 | def make_3d_from_2d(array_2d,x,y,dz): |
---|
| 114 | |
---|
| 115 | k_tmp = np.arange(0,max(array_2d.flatten())+dz*2,dz) |
---|
| 116 | |
---|
| 117 | k_tmp[1:] = k_tmp[1:] - dz * 0.5 |
---|
| 118 | array_3d = np.ones((len(k_tmp),len(y),len(x))) |
---|
| 119 | |
---|
| 120 | for l in range(0,len(x)): |
---|
| 121 | for m in range(0,len(y)): |
---|
| 122 | for k in range(0,len(k_tmp)): |
---|
| 123 | if k_tmp[k] > array_2d[m,l]: |
---|
| 124 | array_3d[k,m,l] = 0 |
---|
| 125 | |
---|
| 126 | return array_3d.astype(np.byte), k_tmp |
---|
| 127 | |
---|
| 128 | |
---|
| 129 | def make_3d_from_bridges_2d(array_3d,array_2d,x,y,dz,width,fill): |
---|
| 130 | |
---|
[4250] | 131 | for l in range(0,len(x)): |
---|
| 132 | for m in range(0,len(y)): |
---|
[3944] | 133 | if array_2d[m,l] != fill: |
---|
[4250] | 134 | print(str(l) + "/" + str(m) + "/") |
---|
[4021] | 135 | k_min = max( int((array_2d[m,l] - width)/dz), 0 ) |
---|
[3944] | 136 | k_max = int(round(array_2d[m,l]/dz)) |
---|
| 137 | array_3d[k_min:k_max+1,m,l] = 1 |
---|
| 138 | |
---|
| 139 | |
---|
| 140 | return array_3d.astype(np.byte) |
---|
| 141 | |
---|
| 142 | |
---|
| 143 | def check_arrays_2(array1,array2,fill1,fill2): |
---|
| 144 | |
---|
| 145 | missing1 = np.where(array1 == fill1,1,0) |
---|
| 146 | missing2 = np.where(array2 == fill2,1,0) |
---|
| 147 | result = np.array_equal(missing1,missing2) |
---|
| 148 | |
---|
| 149 | return result |
---|
| 150 | |
---|
| 151 | def check_consistency_3(array1,array2,array3,fill1,fill2,fill3): |
---|
| 152 | |
---|
| 153 | tmp_array = np.where(array1 != fill1,1,0) + np.where(array2 != fill2,1,0) + np.where(array3 != fill3,-1,0) |
---|
| 154 | |
---|
| 155 | test = np.any(tmp_array.flatten() != 0) |
---|
| 156 | if test: |
---|
| 157 | print("soil_type array is not consistent!") |
---|
| 158 | print("max: " + str(max(tmp_array.flatten())) + ", min: " + str(min(tmp_array.flatten()))) |
---|
| 159 | else: |
---|
| 160 | print("soil_type array is consistent!") |
---|
| 161 | return tmp_array, test |
---|
| 162 | |
---|
| 163 | |
---|
| 164 | |
---|
| 165 | |
---|
| 166 | def check_consistency_4(array1,array2,array3,array4,fill1,fill2,fill3,fill4): |
---|
| 167 | |
---|
| 168 | tmp_array = np.where(array1 != fill1,1,0) + np.where(array2 != fill2,1,0) + np.where(array3 != fill3,1,0) + np.where(array4 != fill4,1,0) |
---|
| 169 | |
---|
| 170 | test = np.any(tmp_array.flatten() != 1) |
---|
| 171 | if test: |
---|
| 172 | print("*_type arrays are not consistent!") |
---|
| 173 | print("max: " + str(max(tmp_array.flatten())) + ", min: " + str(min(tmp_array.flatten()))) |
---|
| 174 | else: |
---|
| 175 | print("*_type arrays are consistent!") |
---|
| 176 | return tmp_array, test |
---|
| 177 | |
---|