Issues with asc raster math

[moovida]

Hi, 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.
• 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 tere a better way to create a raster based on an existing in order to have exact the same grid?)

[jericks]

I am not exactly sure what is going on, but I added two units tests to make sure ArcGrid Rasters work with the setValue and getValue commands.

2e55d6a

The only thing I noticed was:

def outChmRaster = Format.getFormat(dsmPath).read()

Should dsmPath be chmPath?

[moovida]

def outChmRaster = Format.getFormat(dsmPath).read()

Should dsmPath be chmPath?

Actually not. This is the trick I am using in order to not have to create the Raster manually (part of point 2). Usually one wants to have a new raster of the exact grid of the used rasters, so the simplest is to re-read one of the rasters and use it to overwrite. Is there a better way using the data of an existing raster?

I will look into your commit and see if it works for me.

Mind that the first point (0.5 resolution) is a bug in my opinion. Your testcases do not consider this, since the rasters you create are not aware of cell resolution:

For example in this case:

        Bounds bounds = new Bounds(0, 0, 7, 5, "EPSG:4326")
        List data = [
                [0,0,0,0,0,0,0],
                [0,1,1,1,1,1,0],
                [0,1,2,3,2,1,0],
                [0,1,1,1,1,1,0],
                [0,0,0,0,0,0,0]
        ]
        Raster raster = new Raster(data, bounds)

what is resolution? 1?

[moovida]

I created two mini ASC files to reproduce in a testcase the cell resolution problem, which I think is a bug:

CHM.zip

Script used:

import geoscript.layer.*

def base = "/your/path/to/data/"
def dsmPath = base + "test_dsm.asc"
def dtmPath = base + "test_dtm.asc"
def chmPath = base + "test_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)

Instead with those two it seem to not have the strange diff problem I was experiencing. Will need to check that on my side.

[moovida]

Last note about the problem 2 in the original issue. The script still doesn't work on large rasters. I have no idea why. But it does if I use the data list approach.

import geoscript.layer.*
import geoscript.geom.*

def base = "path"
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 outChmRasterData = []

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

for ( row in 0..(rows-1)) {
    def rowData = []
    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
            if(chmValue < 1.0) {
                chmValue = nv
            }
            rowData << chmValue
        } else {
            rowData << nv
        }
    }
    outChmRasterData << rowData
}

def outChmRaster = new Raster(outChmRasterData, dsmRaster.getBounds())
Format.getFormat(chmPath).write(outChmRaster)

This is just for the record, at the moment I am not able to understand why.

[jericks]

I think I may have discovered the problem. If you try to set values on a raster that exceed the maximum value of the existing raster, the values seem to wrap when you write to disk. If you save the same raster to a new file using a new Format it doesn't wrap and it seems to work.

This script works on a larger raster (the natural earth terrain raster).

https://gist.github.com/jericks/bb84e294b25afc89d6fda5a74dbe6237

I think we should look into the GeoTools GridCoverageBuilder as a way to creating new Raster.

http://docs.geotools.org/stable/javadocs/org/geotools/coverage/grid/GridCoverageBuilder.html

[moovida]

If you try to set values on a raster that exceed the maximum value of the existing raster, the values seem to wrap when you write to disk.

By maximum value you mean the row or col? I am not doing that, so I am not sure I understand what you mean.

Where you able to guess why the 0.5 resolution breaks the getValue in the right border positions?

[jericks]

Bands have minimum and maximum values. I think when you try to write a Raster back to an existing Raster with these values are enforced (not by GeoScript but by GeoTools).

I am not sure about the 0.5 resolution issue yet.

[moovida]

Ah, now I understand, the min/max are on the value itself. That makes sense then. Thank you!
Then the only right way to do it is by creating a new raster with the list of data.

[jericks]

I would like to support a better way to dealing with this. Does this help?

7cf901b