[pdal] Fwd: calculate CHM

[Luca Delucchi]

Hi Andrew,

On 22 June 2015 at 14:42, Andrew Bell <andrew.bell.ia at gmail.com> wrote:
>
>
>
> Filters can feed one another:
>

ok

> <?xml version="1.0"?>
> <Pipeline version="1.0">
> <Writer type="writers.las")
>   <Option name="filename">out.las</Option>
>   <Filter type="filters.programmable">
>   <Option name="source">
> import numpy as np
> def myfunc(ins,outs)
>   val = ins['Y']
>   outs['Y'] = val * 1.05
>   return True
>   </Option>
>   <Filter type="filters.crop>
>     <Option name="bounds">
>         ([1,1.1],[2,2.2])
>     </Option>
>     <Reader name="readers.las"
>       <Option="filename">in.las</Option>
>     </Reader>
>   </Filter>
>   </Filter>
> </Writer>
> </Pipeline>
>

I'm using the following xml file

<Pipeline version="1.0">
<Writer type="writers.las">
<Option name="filename">/my/output/chm.las</Option>
<Filter type="filters.programmable">
<Option name="function">chm</Option>
<Option name="module">anything</Option><Option name="source">
import numpy as np
import struct
import osgeo.gdal as gdal

def get_value(x,y, band, band_type, geomt):
    px = int((x - geomt[0]) / geomt[1])  # x pixel
    py = int((y - geomt[3]) / geomt[5])  # y pixel
    try:
        structval = band.ReadAsArray(px, py, 1, 1)
        val = structval[0][0]
        if cmp(val,0) == -1:
            val = 0
        return val
    except:
        return None

def chm(ins,outs):
    inrast = '/my/input/dtm.tif'
    rast = gdal.Open(inrast)
    band = rast.GetRasterBand(1)
    geomtransf = rast.GetGeoTransform()
    band_type = band.DataType
    Zs = ins['Z']
    Xs = ins['X']
    Ys = ins['Y']
    newZ = []
    import pdb; pdb.set_trace()
    for i in range(len(Xs)):
        try:
            z = get_value(Xs[i], Ys[i], band, band_type, geomtransf)
        except:
            z = None
        if z:
            nz = Zs[i] - z
            newZ.append(nz)
        else:
            newZ.append(None)
    pdb.set_trace()
    outs['Z'] = np.array(newZ)
    return True
</Option>
<Filter type="filters.crop">
<Option name="polygon">POLYGON ((711348.5 5129215.5,711348.5
5132163.5,713522.5 5132163.5,713522.5 5129215.5,711348.5
5129215.5))</Option>
<Reader type="readers.las">
<Option name="filename">/my/input/file.las</Option>
</Reader>
</Filter>
</Filter>
</Writer>
</Pipeline>

 and it seems to work correctly, using pdb I can see the values of newZ variable

min(newZ), max(newZ)
(None, 69.119999999999891)

but when I check the output file I got only 0 for the Z values

{
          "average": 0,
          "count": 8227397,
          "maximum": 0,
          "minimum": 0,
          "name": "Z",
          "position": 2
},


> We're working on getting more examples written and posted in the
> documentation.  Write back if you have questions.
>

I can contribute with mine when they will work


-- 
ciao
Luca

http://gis.cri.fmach.it/delucchi/
www.lucadelu.org