0

I have a DEM raster with some negative values. After overlying other shapefiles, I realize that the negative values are for bodies of water (mainly lakes and bays). I want to replace the negative values with the positive values elevation values (lakes). Any suggestion to do this in python? (Con tool will not work: lakes are at different altitude levels).

6
  • Thanks for your anwers. However I was (and still looking) for doing something like this:
    – Delonix R.
    Commented Feb 26, 2016 at 17:10
  • Could you tell us where these replacement values come from? Are they in another DEM? Or are they perhaps (as the phrasing suggests) the absolute values of the negative numbers already present? If this is so, why not just use the absolute value function?
    – whuber
    Commented Feb 26, 2016 at 18:50
  • Hi whuber, the values should come from the surroundings (neighbourhood) of the negative values, in the same DEM. Because some negative values are clustered, usually they are lakes or quarries in the reality (terrain). So my goal was to replace negative values with the positive -or zero- values that surround them. I hope I was able to explain myself. Thanks.
    – Delonix R.
    Commented Feb 26, 2016 at 18:56
  • Finding those surrounding values is the hardest part of this question. It is important that you mention this in the question itself! (I'm not implying it's particularly difficult, though. For instance, you can Regiongroup the negative values, Expand those regions by one pixel, and compute the Zonal minima of the regions after replacing all negative values in the DEM with Nodata.)
    – whuber
    Commented Feb 26, 2016 at 19:10
  • Yes, below that solution was mentioned. And I've already solved that part myself in pretty much the same way... but what I really want to know if you can modify any given cell in the raster, like using python/Numpy or other similar programming tools. Are you aware of this?
    – Delonix R.
    Commented Feb 26, 2016 at 19:25

5 Answers 5

3

I would think you could still use the CON tool if you do a little bit of vector work first. If you have a layer of the lakes, you would want them to have an attribute for the base elevation value (as shown in my embedded graphic). I drew this inside a study area polygon because you might want to avoid problems of NoData values by having the study area (green) be 0 or some value to use in the con statement.

enter image description here

So lets say you convert that vector work into a raster and call it studyArea. It has values of 0, 200, or 300 (in my example).

import arcpy
from arcpy import env
from arcpy.sa import *
env.workspace = "C:/myFolder/data"
outCon = Con(Raster("studyArea") > 0, "studyArea", "dem")
outCon.save("C:/myFolder/output/output.img")
2

Assuming your lakes are polygons:

  1. Assign unique IDs to your lakes and do zonal statistics of DEM, maximum. In theory it should pick elevation of individual lake shores
  2. If output has negatives, create small buffer for lakes polygon and repeat step 1.
  3. Use Con(IsNull("maxdem"),"dem", "maxdem"), or Con("dem"<0),"maxdem", "dem")

Don't forget to set environment extent=dem extent

2

I prefer the min instead of the max, as a lake would otherwise overflow BUT beware of dams where the structure could be invisible on coarse resolution DEM's (if you expect large dams in your study area , the mean is safer)

first, extract your lakes with their bank

  1. Con("DEM"<0, 1, 0) to get each individual lake

  2. region group to assign an unique ID to each lake (and background)

  3. shrink the background by one pixel

Then, get the minimum value of the DEM around the lake

  1. Con("DEM"<0, 9999, "DEM") replace negative DEM values with very large values

  2. Zonal statistics "MIN" gets the minimum pixel values on the lake bank

Finally, a combine min values with original DEM

Con("DEM"<0, "zonalMin", "DEM")
3
  • Shrinking can be tricky for multiple lakes
    – FelixIP
    Commented Feb 18, 2016 at 20:23
  • this is why I shrink the background instead of expanding the lakes. Alternatively you can use a 3 by 3 focal stat (minimum) on the raster where lakes are filled with 9999 and perform the zonal stat on non expanded lakes.
    – radouxju
    Commented Feb 18, 2016 at 20:43
  • Shrinking the background, as described here, guarantees you get no information out of the lakes at all. If you're worried about connecting adjacent lakes, then perform a focal minimum within a 3 X 3 neighborhood and find the zonal minimum of that within each lake region.
    – whuber
    Commented Feb 26, 2016 at 22:05
1
  1. If you don't already have polygons for your lakes, you can pull them out of the DEM by finding clusters of negative values. Identify and extract the negative or 0 cells in the raster using the Extract By Attributes tool.
  2. Convert the lakes raster to a polygon.
  3. Buffer the polygon by whatever radius you want to make up your lake "neighborhood." Create a ring-shaped buffer to leave out the negative cells. Use "OUTSIDE_ONLY" for the line side parameter.
  4. Here's where you have options on what value you want to use to fill in the cells. You can either (a) assign the average (or minimum) value of cells in the buffer to the lake or (b) you can create an interpolated raster to estimate what the landscape would be in the absence of the lakes. FelixIP and radouxju addressed the first option. If you choose to do the second option, you can use your ring-shaped buffer polygons to extract only the cells in your neighborhood using Extract by Mask.
  5. Convert the raster to points. This should produce rings of points around your lakes with the shore elevation values.
  6. Use one of the Interpolation tools to estimate the values inside the lake polygons, using the lake polygons as a mask for output extent. There are many different methods to do interpolation, but if you don't need anything fancy, I prefer using Inverse Distance Weighted because it gives greater weight to points close by and less weight to points far away.
  7. Replace the cells in your original DEM using a Con statement:

    Con("DEM"<0, "InterpolatedRas", "DEM")

This workflow requires Spatial Analyst.

This will produce a raster with varying values in the lake areas. In essence, it estimates what the elevations of the land would be based on the pattern of elevations surrounding the lake. This will be more pronounced if your DEM is full of hills and cliffs. If the surrounding landscape is relatively flat, then this is probably overkill and you can go with assigning a single value for the whole lake.

-1

By using Numpy, you can get something like this:

x    y    value
0    0    0
0    1    0
0    2    1
1    0    0
1    1    1

where x and y are row and columns in the raster, and value is the value of that precise cell in the raster (see code below). This way we can easily spot the unwanted values, make a change to the wanted value and put that value back to the raster. Now, the question is: How can we put that value back to the raster???

Code to obtain rows, columns and values of all the raster's cells.

import arcpy

fOut = open('outputFile.txt', w)                # Open output file

fOut.write('x' + '\t' + 'y' + '\t' + 'value\n') # Write the header

rstArray = arcpy.RasterToNumPyArray(rasterFile) # Change rasterFile to numpy array

rows, cols = rstArray.shape                     # Return the rows, columns

for rowNum in xrange(rows):                     # Loop through the rows

    for colNum in xrange(cols):                 # Loop through the row's columns

        value = rstArray.item(rowNum, colNum)   # Get the value at the cell

        fOut.write(str(rowNum) + '\t')          # Write the row number
        fOut.write(str(colNum) + '\t')          # Write the column number
        fOut.write(str(value) + '\n')           # Write the value and a new line
fOut.close()

This code is in the question: "How to get values of each cell in raster attribute table?" in this same website/forum.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.