Raster calculation issues

[andrea antonello]

Hi all,
I am trying to do a very simple raster diff (manually to teach
students). The script is the following:

import geoscript.layer.*

def base = "path/to/base/folder/"
def dsmPath = base + "DSM_SolarTirol_small.asc"
def dtmPath = base + "DTM_SolarTirol_small.asc"
def chmPath = base + "CHM.asc"

def dtmRaster = Format.getFormat(dtmPath).read()
def dsmRaster = Format.getFormat(dsmPath).read()
def outChmRaster = Format.getFormat(dsmPath).read()

def cols = dtmRaster.cols
def rows = dtmRaster.rows
def nv = -9999.0

for ( row in 0..(rows-1)) {
for ( col in 0..(cols-1)) {
def position = [col, row]
def dtmValue = dtmRaster.getValue(position)
def dsmValue = dsmRaster.getValue(position)
if(dtmValue != nv && dsmValue != nv){
def chmValue = dsmValue - dtmValue
outChmRaster.setValue(position, chmValue)
}
}
}

Format.getFormat(chmPath).write(outChmRaster)

Issues with this:

- if the raster has a cell size of 0.5, I get an error due to the fact
that the point of the last column is outside the raster. It doesn't
happen if I cheat and change it to 1.0. I am guessing there is a
rounding problem. (opened an issue here:
https://github.com/geoscript/geoscript-groovy/issues/40 )
- in general the script doesn't work. I am not sure what I am not
seeing, but the result is the same as the dsm raster. I am indeed
using the dsmraster as a template to fill in the new data, but the
setValue should overwrite it, which it is not doing.

is there a better way to create a raster based on an existing in order
to have exact the same grid?


Thanks,
Andrea