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