Hi all,

A while ago I had to do a similar task for Lidar analysis.... So I
made a quick python program for it. The only problem is that it spits
out boxes and not nice round shapes. The code will create both GRASS
ASCII vector format or a csv file that can be used to subset las data.
Also, since I created the program specifically for our needs, the
column numbers are fixed.

Anyway, the code is attached. Hope it helps.

Cheers
Daniel

PS1 - I'm no professional coder so I'm certain there are flaws and
leftover lines in the code. But it worked for me...
PS2 - I have another version that will generate a shapefile with
ellipses as crowns but the code is too messy. I'd like to make it more
presentable before releasing it but, if you can stomach poorly writen
python, just let me know...


On Fri, Oct 25, 2013 at 8:33 AM, Nikos Alexandris
<[email protected]> wrote:
> On Friday 25 of October 2013 07:56:32 Matteo Mura wrote:
>> Dear all,
>>
>> I have a vector-points file containing the points of trees location in a 100
>> x 100 square, it is a permanent plot. The database of this vector contains
>> for each sampled tree the crown projection for each direction (N-S-E-W). I
>> need to convert the points to polygons which shape and size is determined
>> by the crown projection measures. In other therms I would like to have the
>> polygon file which describe the tree crown shape and dimension.
>>
>> Is it possible to do in GRASS? If yes, how can I?
>
> Not a reply, just a comment: this is a very nice question! If this gets a
> smart solution consider posting it at GIS@SE. GRASS GIS deserves/needs to
> attract more points over there.
>
> Nikos
> _______________________________________________
> grass-user mailing list
> [email protected]
> http://lists.osgeo.org/mailman/listinfo/grass-user
# -*- coding: utf-8 -*-
"""
Created on Mon Jun 11 17:45:47 2012

Code to create canopy boxes around stem locations
works for field data from Reserva Duke - Maria Hunter
as exported by Grass db.select!!

Altered on 8/8/2012 to work with directly with the csv data

@author: Daniel Victoria <[email protected]>
Sustainable Landscapes Brazil project
"""

class boxCanopy():
    """Class to create and export canopy box around stem"""
    
    def __init__(self, trans, stem, utm_e, utm_n, cnpN, cnpS, cnpE, cnpW):
        """
        Initialize canopy box.
        stem: Stem number
        utm_e & n: stem position
        cnpyN,S,E,W: crown dimensions
        """
        
        # cleaning extra "s
        stem = stem.replace('"','')
        self.stem = 'stem_'+stem
        self.cat = stem
        self.trans = trans            
        self.utm_e = float(utm_e)
        self.utm_n = float(utm_n)
        self.cnpN = float(cnpN)
        self.cnpS = float(cnpS)
        self.cnpE = float(cnpE)
        self.cnpW = float(cnpW)
        
        # calculating N, S, E, W extremes
        self.N = self.utm_n + self.cnpN
        self.S = self.utm_n - self.cnpS
        self.E = self.utm_e + self.cnpE
        self.W = self.utm_e - self.cnpW
    
    def spitAscii(self):  
        """Return canopy box in Grass ascii format"""
        if self.skip == 1:
            #print("No data for this point")
            return
        linha = []
        linha.append("B 5\n")
        linha.append("  %.8f %.8f\n" % (self.W, self.S))
        linha.append("  %.8f %.8f\n" % (self.E, self.S))
        linha.append("  %.8f %.8f\n" % (self.E, self.N))
        linha.append("  %.8f %.8f\n" % (self.W, self.N))
        linha.append("  %.8f %.8f\n" % (self.W, self.S))
        linha.append("C 1 1\n")
        linha.append("  %.8f %.8f\n" % (self.utm_e, self.utm_n))
        linha.append("  1 %s\n" % (self.cat))
        return linha
        
    def lasSubset(self):
        """Return csv file appropriate for subseting las files
        using the QRSC Lidar program.
        CSV must have: id, easting, northing, major, minor axis, azimuth
        When subsetting rectangles, major axis is N-S dir
        azimuth is clockwise from north"""
#        if self.skip == 1:
#            #print("No data for this point")
#            return
        l = (self.stem, (self.E + self.W)/2, (self.N + self.S)/2, (self.cnpN + self.cnpS)/2, (self.cnpE + self.cnpW)/2)
        linha = ("%s, %.8f, %.8f, %.4f, %.4f, 0\n" % l)
        return linha
    
        
        
def processGrass(infile, outfile):
    """Process file with the data for GRASS output"""
    entrada = open(infile,'r')
    #skipping header
    entrada.readline()
    dados = entrada.readlines()
    entrada.close()
    saida = open(outfile, 'w')
    saida.write("VERTI:\n")
    for i in dados:
        p = i.split('|')
        if [p[2], p[4], p[5], p[13], p[14], p[15], p[16]].count('NA') == 0:
            # there is data to process...
            ponto = boxCanopy(p[1], p[2], p[4], p[5], p[13], p[14], p[15], p[16])
            linhas = ponto.spitAscii()
            for i in linhas:
                saida.write(i)
    saida.close()
    return

def processLasSubset(infile, outfile, point_filter = ''):
    """Process file with the data for las subset
    point_filter is a filter to export only data for a specific plot
    sometimes it has to be surrounded by quotes
    c is a list with the col numners of the needed data eg.
    c = [trans, stem, utmE, utmN, crownN, crownS, crownE, crownW]"""
    
    # c for duke_v4, 2009, 2011 and average, tapajos
    #c = [0, 1, 2, 3, 15, 16, 17, 18]
    #c = [0, 1, 2, 3, 28, 29, 30, 31]
    #c = [0, 1, 2, 3, 33, 34, 35, 36]
    c = [0, 1, 2, 3, 4, 5, 6, 7]
    
    entrada = open(infile,'r')
    #skipping header
    entrada.readline()
    dados = entrada.readlines()
    entrada.close()
    saida = open(outfile, 'w')
    for i in dados:
        p = i.split(',')
        if [p[c[1]], p[c[2]], p[c[3]], p[c[4]], p[c[5]], p[c[6]], p[c[7]]].count('NA') == 0:
            # there is data to process...
            ponto = boxCanopy(p[c[0]], p[c[1]], p[c[2]], p[c[3]], p[c[4]], p[c[5]], p[c[6]], p[c[7]])
            if point_filter == '':
                #process all data
                linha = ponto.lasSubset()
                saida.write(linha)
            else:
                #process only selected data
                if ponto.trans == point_filter:
                    linha = ponto.lasSubset()
                    saida.write(linha)
    saida.close()
    return
_______________________________________________
grass-user mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/grass-user

Reply via email to