[3567] | 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 | # |
---|
| 17 | # Copyright 1997-2018 Leibniz Universitaet Hannover |
---|
| 18 | #--------------------------------------------------------------------------------# |
---|
| 19 | # |
---|
| 20 | # Current revisions: |
---|
| 21 | # ----------------- |
---|
| 22 | # |
---|
| 23 | # |
---|
| 24 | # Former revisions: |
---|
| 25 | # ----------------- |
---|
| 26 | # $Id: palm_csd_tools.py 3726 2019-02-07 18:22:49Z maronga $ |
---|
[3726] | 27 | # Index bound bugfix |
---|
| 28 | # |
---|
| 29 | # 3668 2019-01-14 12:49:24Z maronga |
---|
[3668] | 30 | # Minor changes |
---|
| 31 | # |
---|
| 32 | # 3567 2018-11-27 13:59:21Z maronga |
---|
[3567] | 33 | # Initial revisions |
---|
| 34 | # |
---|
| 35 | # |
---|
| 36 | # Description: |
---|
| 37 | # ------------ |
---|
| 38 | # Support routines for palm_csd |
---|
| 39 | # |
---|
| 40 | # @Author Bjoern Maronga (maronga@muk.uni-hannover.de) |
---|
| 41 | # |
---|
| 42 | #------------------------------------------------------------------------------# |
---|
| 43 | import numpy as np |
---|
| 44 | from scipy.interpolate import interp2d |
---|
| 45 | |
---|
| 46 | def blend_array_2d(array1,array2,radius): |
---|
| 47 | # Blend over the parent and child terrain height within a radius of 50 px |
---|
| 48 | |
---|
| 49 | gradient_matrix = np.copy(array1) |
---|
| 50 | gradient_matrix[:,:] = 1.0 |
---|
| 51 | gradient_px = 50 |
---|
| 52 | |
---|
| 53 | for i in range(0,gradient_px): |
---|
| 54 | gradient_matrix[:,i] = float(i)/float(gradient_px) |
---|
| 55 | |
---|
| 56 | |
---|
| 57 | for i in range(len(gradient_matrix[0,:])-gradient_px,len(gradient_matrix[0,:])): |
---|
| 58 | gradient_matrix[:,i] = float(len(gradient_matrix[0,:])-i)/float(gradient_px) |
---|
| 59 | |
---|
| 60 | |
---|
| 61 | for j in range(0,gradient_px): |
---|
| 62 | for i in range(0,len(gradient_matrix[0,:])): |
---|
| 63 | if gradient_matrix[j,i] == 1.0: |
---|
| 64 | gradient_matrix[j,i] = float(j)/float(gradient_px) |
---|
| 65 | else: |
---|
| 66 | gradient_matrix[j,i] = (gradient_matrix[j,i] + float(j)/float(gradient_px))/2.0 |
---|
| 67 | |
---|
| 68 | |
---|
| 69 | for j in range(len(gradient_matrix[:,0])-gradient_px,len(gradient_matrix[:,0])): |
---|
| 70 | for i in range(0,len(gradient_matrix[0,:])): |
---|
| 71 | if gradient_matrix[j,i] == 1.0: |
---|
| 72 | gradient_matrix[j,i] = (len(gradient_matrix[:,0])-j)/float(gradient_px) |
---|
| 73 | else: |
---|
| 74 | gradient_matrix[j,i] = (gradient_matrix[j,i] + (len(gradient_matrix[:,0])-j)/float(gradient_px))/2.0 |
---|
| 75 | |
---|
| 76 | array_blended = array1 * gradient_matrix + (1.0 - gradient_matrix ) * array2 |
---|
| 77 | |
---|
| 78 | return array_blended |
---|
| 79 | |
---|
| 80 | |
---|
| 81 | def interpolate_2d(array,x1,y1,x2,y2): |
---|
[3668] | 82 | tmp_int2d = interp2d(x1,y1,array,kind='linear') |
---|
| 83 | array_ip = tmp_int2d(x2.astype(float), y2.astype(float)) |
---|
[3567] | 84 | |
---|
[3668] | 85 | return array_ip |
---|
[3567] | 86 | |
---|
| 87 | |
---|
| 88 | def bring_to_palm_grid(array,x,y,dz): |
---|
| 89 | |
---|
| 90 | # Bring the parent terrain height to the child grid |
---|
| 91 | k_tmp = np.arange(0,max(array.flatten())+dz*2,dz) |
---|
| 92 | k_tmp[1:] = k_tmp[1:] - dz * 0.5 |
---|
| 93 | for l in range(0,len(x)): |
---|
| 94 | for m in range(0,len(y)): |
---|
[3726] | 95 | for k in range(1,len(k_tmp+1)): |
---|
[3567] | 96 | if k_tmp[k] > array[m,l]: |
---|
[3726] | 97 | array[m,l] = k_tmp[k-1] |
---|
[3567] | 98 | break |
---|
| 99 | |
---|
| 100 | return array |
---|
| 101 | |
---|
| 102 | |
---|
| 103 | def make_3d_from_2d(array_2d,x,y,dz): |
---|
| 104 | |
---|
| 105 | k_tmp = np.arange(0,max(array_2d.flatten())+dz*2,dz) |
---|
| 106 | |
---|
| 107 | k_tmp[1:] = k_tmp[1:] - dz * 0.5 |
---|
| 108 | array_3d = np.ones((len(k_tmp),len(y),len(x))) |
---|
| 109 | |
---|
[3726] | 110 | for l in range(0,len(x)): |
---|
| 111 | for m in range(0,len(y)): |
---|
| 112 | for k in range(0,len(k_tmp)): |
---|
[3567] | 113 | if k_tmp[k] > array_2d[m,l]: |
---|
| 114 | array_3d[k,m,l] = 0 |
---|
| 115 | |
---|
| 116 | return array_3d.astype(np.byte), k_tmp |
---|
| 117 | |
---|
| 118 | |
---|
| 119 | def make_3d_from_bridges_2d(array_3d,array_2d,x,y,dz,width,fill): |
---|
| 120 | |
---|
| 121 | for l in range(0,len(x)-1): |
---|
| 122 | for m in range(0,len(y)-1): |
---|
| 123 | if array_2d[m,l] != fill: |
---|
| 124 | k_min = max( int(array_2d[m,l] - width)/dz, 0 ) |
---|
| 125 | k_max = int(round(array_2d[m,l]/dz)) |
---|
| 126 | array_3d[k_min:k_max+1,m,l] = 1 |
---|
| 127 | |
---|
| 128 | |
---|
| 129 | return array_3d.astype(np.byte) |
---|
| 130 | |
---|
| 131 | |
---|
| 132 | def check_arrays_2(array1,array2,fill1,fill2): |
---|
| 133 | |
---|
| 134 | missing1 = np.where(array1 == fill1,1,0) |
---|
| 135 | missing2 = np.where(array2 == fill2,1,0) |
---|
| 136 | result = np.array_equal(missing1,missing2) |
---|
| 137 | |
---|
| 138 | return result |
---|
| 139 | |
---|
[3668] | 140 | def check_consistency_3(array1,array2,array3,fill1,fill2,fill3): |
---|
[3567] | 141 | |
---|
[3668] | 142 | tmp_array = np.where(array1 != fill1,1,0) + np.where(array2 != fill2,1,0) + np.where(array3 != fill3,-1,0) |
---|
| 143 | |
---|
| 144 | test = np.any(tmp_array.flatten() != 0) |
---|
| 145 | if test: |
---|
| 146 | print("soil_type array is not consistent!") |
---|
| 147 | print("max: " + str(max(tmp_array.flatten())) + ", min: " + str(min(tmp_array.flatten()))) |
---|
| 148 | else: |
---|
| 149 | print("soil_type array is consistent!") |
---|
| 150 | return tmp_array, test |
---|
| 151 | |
---|
| 152 | |
---|
| 153 | |
---|
| 154 | |
---|
[3567] | 155 | def check_consistency_4(array1,array2,array3,array4,fill1,fill2,fill3,fill4): |
---|
| 156 | |
---|
| 157 | 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) |
---|
| 158 | |
---|
| 159 | test = np.any(tmp_array.flatten() != 1) |
---|
| 160 | if test: |
---|
| 161 | print("*_type arrays are not consistent!") |
---|
| 162 | print("max: " + str(max(tmp_array.flatten())) + ", min: " + str(min(tmp_array.flatten()))) |
---|
| 163 | else: |
---|
| 164 | print("*_type arrays are consistent!") |
---|
| 165 | return tmp_array, test |
---|
| 166 | |
---|