# -*- coding: cp1252 -*- # Define a simple class Point2D for XY-points, polygons are lists of Point2D objects import math, random def diameter_at_height(hx, treeheight, d13, b): """This function returns the diameter at height hx""" # check the sanity of hx, the height of interest if (hx < 0) or (hx > treeheight): print "diameter_at_height(), problem with hx!" return 0 # Do conversions to float to assure "correct" decimal math hx = float(hx); d13 = float(d13); treeheight = float(treeheight) # Solve the diameter at 20% height, using the d13 as reference x = (1-1.3/treeheight) fx = b[0]*x + b[1]*x**2 + b[2]*x**3 + b[3]*x**5 + b[4]*x**8+ b[5]*x**13 + b[6]*x**21 + b[7]*x**34 d2h=d13/fx # Now, solve for diameter dx at height hx x = (1-hx/treeheight) fx = b[0]*x + b[1]*x**2 + b[2]*x**3 + b[3]*x**5 + b[4]*x**8+ b[5]*x**13 + b[6]*x**21 + b[7]*x**34 return d2h*fx def height_for_diameter(dx, treeheight, d13, b): """Finds thru iteration height at which diameter is dx""" # Data type changes to float and initialization of needed vars h = float(treeheight); g = 1.0; step = 0.5 d13 = float(d13) dx = float(dx) # Check the sanity of the parameters if ((0 > dx) or (dx > 85.0) or (1.4 > h) or (h > 40.0)): print "d13:" , dx, " h:" , h, " ? arguments in height for diameter" return 0 # print dx, diameter_at_height(0.01, h, d13, b) if dx > diameter_at_height(0.01, h, d13, b): return 0 # Do a coarse 1-m step-search down from the top and stop when diameter # is greater than dx l = h / step # number of coarse steps hr = h # hr stores the current height if (l * step == h): hr = hr + 0.01 # start point for iteration slightly above top # Coarse search, down from the top, by 1 m step for i in range (1, int(l+1),1): hr = hr - step a = diameter_at_height(hr, h, d13, b) # if the diameter is greater than dx, stop coarse search if (a) > dx: break # Fine tuning, coarse solution (height) now stored in variable hr if a < dx: return 0 # There is no such diameter available! g = 1.0 # direction to go to has to be reversed, must return up the stem diff = 1 # assign value for accuracy variable while (abs(diff)>0.001) : # The wanted accuracy is diff, we continue until this acc. is reached step = step / 2 # reduce the step size at every round of loop hr = hr + g * step # compute new candidate position in the direction of g (-1 or 1) dr = diameter_at_height(hr, h, d13, b) # check diameter for hr diff = dx-dr # difference if (dr > dx): # check if we still are approching from below g = 1.0 if (dr <= dx): # if not, bo back down from above g = -1.0 return hr def log_volume(h1, h2, h, d13, b): """Computes the volume between heights h_low and h_high""" if (h1 > h2): (h1,h2) = (h2, h1) # swap the order if (h1 < 0 or h2 > h): print "Parameter problem in log_volume" return -1 B1 = b[0]; B2 = b[1]; B3 = b[2]; B4 = b[3] B5 = b[4]; B6 = b[5]; B7 = b[6]; B8 = b[7] x = (1-float(h1)/h) # Evaluate integral at low end V1 = 2./37*B2*B8*x**37+1./20*B8*B4*x**40+2./23*B1*B7*x**23+1./24*B6*B8*x**48+ \ 1./12*B7*B2*x**24+2./15*B1*B6*x**15+1./4*B4*B2*x**8+1./7*B5*B4*x**14+ \ 1./2*B1*B2*x**4+1./5*B5*B1*x**10+1./19*B3*B8*x**38+2.0/19*B6*B4*x**19+2./9*B3*B4*x**9 \ +1./3*B2*B3*x**6+1./28*B7*B8*x**56+1./8*B2*B6*x**16+2./25*B7*B3*x**25+1./18*B8*B1*x**36 \ +1./15*B7*B5*x**30+1./11*B6*B5*x**22+1./6*B3*B5*x**12+2./35*B6*B7*x**35+1./27*(B6**2+2*B4*B7)*x**27 \ +1./43*(2*B8*B5+B7**2)*x**43+1./7*(2*B4*B1+B3**2)*x**7+1./5*(2*B1*B3+B2**2)*x**5+1./69*B8**2*x**69 \ +1./17*(B5**2+2*B6*B3)*x**17+1./11*(B4**2+2*B5*B2)*x**11+1./3*B1**2*x**3 x = (1-float(h2)/h) # Evaluate integral at high end V2 = 2./37*B2*B8*x**37+1./20*B8*B4*x**40+2./23*B1*B7*x**23+1./24*B6*B8*x**48+ \ 1./12*B7*B2*x**24+2./15*B1*B6*x**15+1./4*B4*B2*x**8+1./7*B5*B4*x**14+ \ 1./2*B1*B2*x**4+1./5*B5*B1*x**10+1./19*B3*B8*x**38+2./19*B6*B4*x**19+2./9*B3*B4*x**9 \ +1./3*B2*B3*x**6+1./28*B7*B8*x**56+1./8*B2*B6*x**16+2./25*B7*B3*x**25+1./18*B8*B1*x**36 \ +1./15*B7*B5*x**30+1./11*B6*B5*x**22+1./6*B3*B5*x**12+2./35*B6*B7*x**35+1./27*(B6**2+2*B4*B7)*x**27 \ +1./43*(2*B8*B5+B7**2)*x**43+1./7*(2*B4*B1+B3**2)*x**7+1./5*(2*B1*B3+B2**2)*x**5+1./69*B8**2*x**69 \ +1./17*(B5**2+2*B6*B3)*x**17+1./11*(B4**2+2*B5*B2)*x**11+1./3*B1**2*x**3 x = (1-1.3/h) fx = b[0]*x + b[1]*x**2 + b[2]*x**3 + b[3]*x**5 + b[4]*x**8+ b[5]*x**13 + b[6]*x**21 + b[7]*x**34 d2h=d13/fx return (V1-V2) * (d2h/100.0)**2 * h * math.pi/4. def gimme_b(sp, b): if sp == 1: b.append(2.1288),b.append(-0.63157),b.append(-1.6082),b.append(2.4886) b.append(-2.4147),b.append(2.3619), b.append(-1.7539),b.append(1.0817) if sp == 2: b.append(2.3366), b.append(-3.2684) ,b.append(3.6513), b.append(-2.2608) b.append(0), b.append(2.1501),b.append(-2.7412),b.append(1.8876) if sp == 3: b.append(0.93838), b.append(4.106), b.append(-7.8517), b.append(7.8993) b.append(-7.5018), b.append(6.3863), b.append(-4.3918), b.append(2.1604) return b class Point2D: def __init__(self,x,y): self.x = float(x) self.y = float(y) # This function returns TRUE if point Px is inside polygon P, P is assumed to # have the first point to appear twice, first and last point in the list. def InsidePolygon(P, Px): """ Function says if Px is inside P[] consisting of N elements""" counter = 0 i = 0 xinters = 0.0 N = len(P) p1 = Point2D(P[0].x,P[0].y) for i in range (1, N+1,1): p2 = Point2D(P[i%N].x,P[i%N].y) if (Px.y > min(p1.y, p2.y)): if (Px.y <= max(p1.y, p2.y)): if (Px.x <= max(p1.x, p2.x)): if (p1.y != p2.y): xinters = (Px.y - p1.y) * (p2.x - p1.x) / (p2.y - p1.y) + p1.x if ((p1.x == p2.x) or (Px.x <= xinters)): counter = counter + 1 p1 = p2 if (counter%2 == 0): InsidePolygon = False else: InsidePolygon = True return InsidePolygon # *******************************************************************************’ # point distance to (line in) a polygon # Adopted from function: p_poly_dist() in Matlab's matlabcentral def point_to_polyline(p,Px): """Distance of Point2D Px to line segments in p[]""" # Lists that are needed to hold the points and comutations (vector algerba) A=range(len(p)-1);B=range(len(p)-1);C=range(len(p)-1) AB=range(len(A));vv=range(len(A));A_AB=range(len(A));B_AB=range(len(A)); xp =range(len(A)); yp = range(len(A)) # Differences of coordinates for i in range(len(p)-1): A[i] = -1.0 * (p[i+1].x - p[i].x) B[i] = 1.0 * (p[i+1].y - p[i].y) C[i] = p[i+1].y*p[i].x - p[i+1].x*p[i].y # Find the orthogonal projection of point (x,y) on each line segment for i in range(len(A)): AB[i] = 1/(A[i]**2+B[i]**2) vv[i] = A[i]*Px.x + B[i]*Px.y + C[i] for i in range(len(A)): A_AB[i] = A[i]*AB[i] B_AB[i] = B[i]*AB[i] # xp and yp hold the xy intersections for i in range(len(A)): xp[i] = Px.x - A_AB[i]*vv[i] yp[i] = Px.y - B_AB[i]*vv[i] # Find all cases where the projected point is inside the line segment, store indeces in idx[] idx = [] for i in range(len(xp)): for j in range(len(p)-1): if (xp[i] >= p[j].x and xp[i] <=p[j+1].x ) or (xp[i]>=p[j+1].x and xp[i] <= p[j].x): if (yp[i] >= p[j].y and yp[i] <=p[j+1].y ) or (yp[i]>=p[j+1].y and yp[i] <= p[j].y): idx.append(i) # Distance from point Px.x,Px.y to the vertices (end points), store in dv[] dv = range(len(p)-1) # We could also append the entries, here we know the length for i in range(len(p)-1): dv[i] = math.sqrt((p[i].x-Px.x)**2 + (p[i].y-Px.y)**2) if(len(idx)==0): # all projections are outside of polygon ribs, the minimum to vertices d = min(dv) else: # Distance from point Px.x, Px.y to the projection on ribs exists, compute them dp = [] for i in range(len(idx)): dp[i] = math.sqrt( (p[idx].x-Px.x)**2 + (p[idx].y-Px.y)**2) d = min(min(dv), min(dp)) return d class OutputVariables(): def __init__(self): self.Volume = None self.G = None self.Dg = None self.Hg = None class Simulation(): def __init__(self): self.N_simu = None self.N_plots = None self.Radius = None self.SpErr = None self.GaussVmodelSigmaMu = None # Sdev of between stand bias (volume model) self.GaussVmodelSigma = None # Sdev of within-stand noise (volume model self.Dist_mu = None self.Dist_sigma = None self.Dbh_mu = None self.Dbh_sigma = None self.Height_mu = None self.Height_sigma = None self.PolylineFilename = None self.TreeFilename = None self.OutputFilename = None class Tree(): def __init__(self, id, x, y, z, sp, status, height, dbh, hc): self.id = int(id) self.x = float(x) self.y = float(y) self.z = float(z) self.sp = int(sp) self.sp_estim = None self.status = int(status) self.dbh = float(dbh) self.dbh_estim = None self.h = float(height) self.h_estim = None self.hc = float(hc) def getV(self): # This method gives the volume of the stem, using a simple formula # "Smalian II" # dbh = cm, h = m return 0.4 * (self.dbh_estim/100.0) ** 2 * self.h_estim def input_trees(yet_another_filename): Trees = [] infile = open(yet_another_filename,'r') i = 0 for line in infile: i = i + 1 id, x, y,z,sp, status, dbh, height, hc = line.split(',') Trees.append(Tree(id, x, y, z, sp, status, height, dbh, hc)) print "I read ", i, "lines" infile.close() return Trees def input_simulation_parameters(filename): Simul=Simulation() infile = open(filename, 'r') for line in infile: splitlist = line.split() if splitlist[0]=="N_simu": Simul.N_simu = int(splitlist[1]) if splitlist[0]=="N_plots": Simul.N_plots = int(splitlist[1]) if splitlist[0]=="Radius": Simul.Radius = float(splitlist[1]) if splitlist[0]=="SpErr": Simul.SpErr = float(splitlist[1]) if splitlist[0]=="GaussVmodelSigmaMu": Simul.GaussVmodelSigmaMu = float(splitlist[1]) if splitlist[0]=="GaussVmodelSigma": Simul.GaussVmodelSigma = float(splitlist[1]) if splitlist[0]=="GaussDist": Simul.Dist_mu = float(splitlist[1]) Simul.Dist_sigma = float(splitlist[2]) if splitlist[0]=="GaussDBH": Simul.Dbh_mu = float(splitlist[1]) Simul.Dbh_sigma = float(splitlist[2]) if splitlist[0]=="GaussHeight": Simul.Height_mu = float(splitlist[1]) Simul.Height_sigma = float(splitlist[2]) if splitlist[0]=="PolylineFilename": Simul.PolylineFilename = splitlist[1] if splitlist[0]=="TreeFilename": Simul.TreeFilename = splitlist[1] if splitlist[0]=="OutputFilename": Simul.OutputFilename = splitlist[1] infile.close() return Simul def input_polyline(PolylineFilename): Polyline = [] i = 0 infile = open(PolylineFilename, 'r') for line in infile: splitlist = line.split(',') x = float(splitlist[0]) y = float(splitlist[1]) Polyline.append(Point2D(x,y)) i +=1 infile.close() print "Read", i, "points" return Polyline def volume_function(sp, d, h): if sp == 1: vol = 0.036089*d**2.01395*(0.99676)**2*h**2.07025*(h-1.3)**(-1.07209) if sp == 2: vol= 0.022927*d**1.91505*(0.99146)**d*h**2.82541*(h-1.3)**(-1.53547) if sp > 2: vol = 0.011197*d**2.10253*(0.986)**d*h**3.98519*(h-1.3)**(-2.659) return vol/1000. # m3 def measure_a_plot(Simul, TreeList, centre, measured_trees): # This function checks by distance measurements, which trees # belong to the plot. # For the trees that belong to the plot, it performs # sp, dbh and height measurements (with Gaussian carelessness) c = centre length_at_enter = len(measured_trees) for t in TreeList: real_distance = math.sqrt((t.x-c.x)**2 + (t.y-c.y)**2) observed_distance = real_distance + random.gauss(Simul.Dist_mu,Simul.Dist_sigma) if observed_distance <= Simul.Radius: # Tree inside, tree in the sample!!! # Let's measure species, assume broadleaved trees 100% but 90 % for pine/spruce if t.sp == 1 or t.sp == 2: # We have a pine or a spruce in reality a = random.random() if a > Simul.SpErr : # Make a mistake, make spruce pine and vice versa if t.sp == 1 : t.sp_estim = 2 elif t.sp == 2: t.sp_estim = 1 else: t.sp_estim = t.sp elif t.sp > 2 : t.sp_estim = 3 else: pass # dbh measurement t.dbh_estim = t.dbh + random.gauss(Simul.Dbh_mu,Simul.Dbh_sigma) # heigh measurement t.h_estim = t.h + random.gauss(Simul.Height_mu,Simul.Height_sigma) measured_trees.append(t) return measured_trees def gimme_a_valid_plot_centre(Simul, Polyline): # define a box that surrounds the forest MinP = Point2D(10**20,10**20) MaxP = Point2D(-10**20,-10**20) for p in Polyline: if p.x < MinP.x : MinP.x = p.x if p.y < MinP.y : MinP.y = p.y if p.x > MaxP.x : MaxP.x = p.x if p.y > MaxP.y : MaxP.y = p.y MinP.x -=1 MinP.y -=1 MaxP.x +=1 MaxP.x +=1 ok = False candidate = Point2D(0,0) while (ok == False): randx = MinP.x + random.random() * (MaxP.x - MinP.x) randy = MinP.y + random.random() * (MaxP.y - MinP.y) candidate = Point2D(randx, randy) #Let's check if candidate is inside InOut = InsidePolygon(Polyline, candidate) if InOut == True : Distance = point_to_polyline(Polyline, candidate) if Distance >= Simul.Radius: ok = True return candidate #Point2D object def Laasis(sp, d, h): if sp == 1: v = 0.036089*d**2.01395*(0.99676)**2*h**2.07025*(h-1.3)**(-1.07209) if sp == 2: v = 0.022927*d**1.91505*(0.99146)**d*h**2.82541*(h-1.3)**(-1.53547) if sp == 3: v = 0.011197*d**2.10253*(0.986)**d*h**3.98519*(h-1.3)**(-2.659) return v/1000.0 def list_mean(lista): tsum = 0 for t in lista: tsum += t return tsum/float(len(lista)) def list_sdev(lista): mean = list_mean(lista) tsum2 = 0 N = len(lista) for t in lista: tsum2 += t**2 return math.sqrt((1/float(N))*tsum2-mean**2) def Simulator(Simul, Polyline, Treelist): i = 0 VolList = [] VolTaperList = [] VolTaperSawList = [] VolTaperPulpList = [] BasalAreaList = [] DgList = [] HgList = [] while (i <= Simul.N_simu): # Stand cruising will start here measured_trees = [] # Stand-level volume model bias, relative ~0.05 BiasModel = random.gauss(0,Simul.GaussVmodelSigmaMu) for j in range(Simul.N_plots): # Here starts a plot measurement # Lets measure a plot, but where? centre = gimme_a_valid_plot_centre(Simul, Polyline) measured_trees = measure_a_plot(Simul, TreeList, centre, measured_trees) # Checking what the program does, so far, what are the contents of the list of measured trees TreesSampled = len(measured_trees) AreaCovered = math.pi * Simul.Radius ** 2 * Simul.N_plots Density = TreesSampled * 10000/AreaCovered # number of trees per hectare = 10000 m2 # print "Stand density", Density, " stems/ha" # ************************************************************************* # Volume, VolumeTaper, Basal area, Hg, Dg SumVol = 0; BasalSum = 0; gd_sum = 0; gh_sum = 0; SumVolTaper = 0 SumVolTaperSaw = 0; SumVolTaperPulp = 0 for t in measured_trees: if t.status < 15 : g = (t.dbh_estim/100.0) ** 2 * math.pi/4 gd = g*(t.dbh_estim/100.0) gh = g*(t.h_estim) # vol = t.getV() vol_saw_taper = 0; vol_pulp_taper = 0 if t.dbh_estim > 3 and t.h_estim > 2: vol = Laasis(t.sp_estim, t.dbh_estim, t.h_estim) # Add stand level bias, multiplicative term vol = vol * (1 + BiasModel) # Add between tree noise vol = vol + vol * random.gauss(0, Simul.GaussVmodelSigma) # vol_taper b = [] b = gimme_b (t.sp_estim, b) vol_taper = log_volume(0.1, t.h_estim, t.h_estim, t.dbh_estim, b) tx1 = height_for_diameter(18.0, t.h_estim, t.dbh_estim, b) if tx1 > 4.1 : vol_saw_taper = log_volume(0.1, tx1, t.h_estim, t.dbh_estim, b) tx2 = height_for_diameter(9.0, t.h_estim, t.dbh_estim, b) if tx2-tx1 > 3: vol_pulp_taper = log_volume(tx1, tx2, t.h_estim, t.dbh_estim, b) gd_sum += gd gh_sum += gh SumVol += vol BasalSum += g SumVolTaper += vol_taper SumVolTaperSaw += vol_saw_taper SumVolTaperPulp += vol_pulp_taper VolList.append(10000/AreaCovered * SumVol) VolTaperList.append(10000/AreaCovered * SumVolTaper) VolTaperSawList.append(10000/AreaCovered * SumVolTaperSaw) VolTaperPulpList.append(10000/AreaCovered * SumVolTaperPulp) BasalAreaList.append(10000/AreaCovered * BasalSum) DgList.append((gd_sum/BasalSum * 100)) HgList.append(gh_sum/BasalSum) capu = "Vol: %0.3f" % list_mean(VolList) capu += ", %0.3f" % list_sdev(VolList) # print capu i += 1 #Stand cruising ends here # Prosessoidaan nämä 5 listaa (alkiot) keskiarvoiksi ja hajonnoiksi. # Tulostus tietostoon: o = open(Simul.OutputFilename, "a") capu = "Vol: \t%7.1f" % list_mean(VolList) capu += " %7.1f m3/ha" % list_sdev(VolList) capu += " %7.1f %%" % ((list_sdev(VolList)/list_mean(VolList)) *100.0) capu1 = "VolTap: %7.1f" % list_mean(VolTaperList) capu1 += " %7.1f m3/ha" % list_sdev(VolTaperList) capu5 = "VolSaw: %7.1f" % list_mean(VolTaperSawList) capu5 += " %7.1f m3/ha" % list_sdev(VolTaperSawList) capu6 = "VolPul: %7.1f" % list_mean(VolTaperPulpList) capu6 += " %7.1f m3/ha" % list_sdev(VolTaperPulpList) capu2 = "G: \t%7.1f" % list_mean(BasalAreaList) capu2 += " %7.1f m2/ha" % list_sdev(BasalAreaList) capu3 = "Dg: \t%7.1f" % list_mean(DgList) capu3 += " %7.1f cm" % list_sdev(DgList) capu4 = "Hg: \t%7.1f" % list_mean(HgList) capu4 += " %7.1f m" % list_sdev(HgList) o.write(capu + "\n") o.write(capu1 + "\n") o.write(capu5 + "\n") o.write(capu6 + "\n") o.write(capu2 + "\n") o.write(capu3 + "\n") o.write(capu4 + "\n") print capu # print capu1 # print capu5 # print capu6 # print capu2 # print capu3 # print capu4 # print "" # ********* Main program ************************* Simul = input_simulation_parameters('c:/temp/simulation.txt') Polyline = input_polyline(Simul.PolylineFilename) TreeList = input_trees(Simul.TreeFilename) j = 0 RadList = [1,2,4,8,16] while(j < len(RadList)): Simul.N_plots = RadList[j] Simulator(Simul, Polyline, TreeList) j +=1