source: palm/trunk/SCRIPTS/palm_csd_files/palm_csd_tools.py @ 4250

Last change on this file since 4250 was 4250, checked in by maronga, 3 years ago

bugfix in palm_csd (calculation of bridge objects)

File size: 5.6 KB
Line 
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 3773 2019-03-01 08:56:57Z maronga $
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
30# Bugfix: missing conversion to integer
31#
32#
33# 3773 2019-03-01 08:56:57Z maronga
34# Unspecified changes
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#------------------------------------------------------------------------------#
53import numpy as np
54from scipy.interpolate import interp2d
55
56def 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
91def 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   
98def 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
113def 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
129def make_3d_from_bridges_2d(array_3d,array_2d,x,y,dz,width,fill):
130     
131      for l in range(0,len(x)):
132         for m in range(0,len(y)):
133            if array_2d[m,l] != fill:
134               print(str(l) + "/" + str(m) + "/")
135               k_min = max( int((array_2d[m,l] - width)/dz), 0 )
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
143def 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
151def 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
166def 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
Note: See TracBrowser for help on using the repository browser.