Some background first, I am working on a model that works with elevation bands. Elevation bands in sense that user defines the number of elevation bands in the area (just like reclassification of DEM under equal intervals). So, in case of my DEM, the minimum and maximum elevation is 1722.35 meters and 4196.3 meters, and if I want the 5 elevation bands, I would do (maximum elevation - minimum elevation)/5 = 495 meters (interval).Moreover, the fishnet/grid raster are approximately 6 km in resolution. So a particular grid cell can have 0-5 elevation bands in it. The Python script that I am using takes four inputs (the data can be accessed from here (the course resolution grid raster, DEM, .txt file) to write data to and the interval as follows:

python format_snow_params.py ~/VIC_GRID_Rstr1.tif ~/DEM.tif ~/snow_params.txt 495

This script performs two passes. The first one loops through the grid cells to find the maximum number of elevation intervals (ie. the maxbands variable). The second pass then gets the actual elevation band values and fractional coverage for those bands. The number of bands are needed for writing the parameter file because VIC expects each cell to have the same number of band information written even though there may not be that many bands in a grid cell.

maxbands is calculated at the end of the first pass then the script moves to the second pass.

The model requires that the sum of fractional area coverages of the elevation bands in a grid cell is equal to one. a snippet of the output file is as follows:

287909  0.5783  0.3905  0.0000  0.0000  0.0000  2404.7144   2561.4426   0.0000  0.0000  0.0000  0.5783  0.3905  0.0000  0.0000  0.0000  
286980  0.7016  0.2673  0.0000  0.0000  0.0000  2407.4819   2518.9944   0.0000  0.0000  0.0000  0.7016  0.2673  0.0000  0.0000  0.0000  

enter image description here The script (works in py 3) is as follows:

from __future__ import print_function
import os
import sys
import warnings
import numpy as np
from osgeo import gdal
from osgeo.gdalnumeric import *
from osgeo.gdalconst import *

# set system to ignore simple warnings

def format_snow_params(basinMask,elvHiRes,outSnow,interval):
    FUNCTION: format_snow_params
    ARGUMENTS: basinMask - path to template raster to run VIC model at
               elvHiRes - path elevation raster dataset at native resolution
               outsnow - path output snow parameter file
               interval - vertical distance to do equal interval segmentation
    KEYWORDS: n/a
    RETURNS: n/a
    NOTES: Does not return a variable but writes an output file

    band = 1 # constant variable for reading in data

    interval = int(interval) # force equal interval value to be int type

    # make a list of input raster files
    infiles = [basinMask,elvHiRes]

    # try to read in the raster data

        # read basin grid raster
        ds = gdal.Open(infiles[0],GA_ReadOnly)
        b1 = ds.GetRasterBand(band)
        mask = BandReadAsArray(b1)
        maskRes = ds.GetGeoTransform()[1]
        ds = None
        b1 = None

        # read hi res elevation raster
        ds = gdal.Open(infiles[1],GA_ReadOnly)
        b1 = ds.GetRasterBand(band)
        elvhires = BandReadAsArray(b1)
        clsRes = ds.GetGeoTransform()[1]
        ds = None
        b1 = None

    # if not working, give error message
    except AttributeError:
        raise IOError('Raster file input error, check that all paths are correct')

    # mask elevation values less than 0
    elvhires[np.where(elvhires<0)] = np.nan

    # get ratio of high resoltion to low resolution
    clsRatio = int(maskRes/clsRes)

    # check if the output parameter file exists, if so delete it
    if os.path.exists(outSnow)==True:

    nbands = [] # blank list

    # try to write to output snow parameter file

        with open(outSnow, 'w') as f:

            cnt = 1 # set grid cell id counter

            # pass counter
            pass_counter = range(2)

            # perform two passes on the raster data
            # 1) to grab the maximum number of bands for a given pixel
            # 2) to calculate the snow band parameters and write to output file
            for pass_cnt in pass_counter:

                # loop over each pixel in the template raster
                for i in range(mask.shape[0]):
                    cy1 = i*clsRatio
                    cy2 = cy1+clsRatio

                    for j in range(mask.shape[1]):
                        cx1 = j*clsRatio
                        cx2 = cx1+clsRatio

                        # get all hi res pixels in a template pixel
                        tmp = elvhires[cy1:cy2,cx1:cx2]
                        if tmp.size == 0:
                        # create blank array for number of bands calculation...
                        if mask[i,j] > 1:
                             # ...if it is not a masked pixel

                            if np.all(tmp == np.nan) == True:
                                tmp[:,:] = 0
                            # find min and max values for interval
                            minelv = np.nanmin(tmp.astype(int)) - (np.nanmin(tmp.astype(int))%interval)
                            maxelv = np.nanmax(tmp.astype(int)) + (np.nanmax(tmp.astype(int))%interval)
                            # create an array of band limits
                            bands = np.arange(minelv, maxelv+interval,interval)

                            bcls = np.zeros_like(tmp)
                            bcls[:,:] = -1

                            # get the number of bands per pixel
                            for b in range(bands.size-1):
                                bcls[np.where((tmp>=bands[b])&(tmp<bands[b+1]))] = b # band counter

                                # if it's the first pass get number of bands for each pixel
                                if pass_cnt == 0:
                                    uniqcnt = np.unique(bcls[np.where(tmp>0)])
                                    nbands.append(uniqcnt.size) # save to a list for second pass

                            if pass_cnt == 1:
                                uniqcnt = np.unique(bcls[np.where(tmp>0)])
                                #clscnt = np.bincount(tmp.ravel())

                                f.write('{0}\t'.format(mask[i,j])) # write grid cell id

                                # find frational area for each band and write to file
                                for c in range(maxbands):
                                        idx = np.where(bcls==uniqcnt[c])
                                        num = np.float(bcls[np.where(bcls>=0)].size)
                                        if num == 0:
                                            num = np.float(idx[0].size)
                                        frac = np.float(idx[0].size) / num
                                    except IndexError:
                                        frac = 0

                                #calculate the mean elevation for each band and write to file
                                for c in range(maxbands):
                                        idx = np.where(bcls==uniqcnt[c])
                                        muelv = np.nanmean(tmp[idx])

                                    except IndexError:
                                        muelv = 0


                                # calculate the precipitation fractions and write to file
                                for c in range(maxbands):
                                        idx = np.where(bcls==uniqcnt[c])
                                        num = np.float(bcls[np.where(bcls>=0)].size)
                                        if num == 0:
                                            num = np.float(idx[0].size)
                                        frac = np.float(idx[0].size) / num
                                    except IndexError:
                                        frac = 0
                                f.write('\n') # write return value for new line

                        if pass_cnt == 1 & mask[i,j] == 1:
                            cnt += 1 # plus one to the grid cell id counter

                if pass_cnt == 0:
                    maxbands = max(nbands) # maximum number of bands for a pixel

        # print the number of bands for user to input into global parameter file
        print('Number of maximum bands: {0}'.format(maxbands))

    # except raise an error when it doesn't work
    except IOError:
        raise IOError('Cannot write output file, error with output snow parameter file path')


def main():
    n_args = len(sys.argv)

    # Check user inputs
    if n_args != 5:
        print ("Wrong user input")
        print ("Script writes the snow band parameter file for the VIC model")
        print ("usage: python format_snow_params.py <template raster> <elevation raster> <output snow band file> <interval for snow bands>")
        print ("Exiting system...")

        # Pass command line arguments into function


# Execute the main level program if run as standalone
if __name__ == "__main__":

Now the script works fine, but the script does not tell how are the elevation bands are spatially distributed in a grid cell (pixel/fishnet polygon of the grid) as the model is not concerned with spatial information. But part of my analysis is create the same elevation bands in ArcMap, such that it tells me the how the elevation bands are spatially distributed in the grid this would enable me to slice the grid cell polygon according to the fractional areas of the elevation bands in each grid cell.

So far I have used the Contour tool in ArcMap with the same interval as used in the script, and set the processing extent to particular grid polygon, but the distribution does not match although the tool creates the same number of elevation bands (5) just like the script. So I am not using the right way in ArcMap.

How can I perform this task in ArcMap?


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.