2

I'm using a script on python "format_snow_param.py". This script takes geotiff data for snow band parameterization and converts it to the needed text format for the VIC model

Usage: python format_snow_params.py <template raster> <elevation raster>
<output snow band file> <interval for snow bands>

on the terminal I type the following command:

Format elevation band file: $ python format_snow_params.py
../data/input/gis/Nyando_grid.tif
../data/input/gis/Nyando_basin_SRTMSnap.tif ../data/input/snow.param 100

The results give this error:

> Traceback (most recent call last):   File "format_snow_params.py",
> line 203, in <module>
>     main()     File "format_snow_params.py", line 197, in main
>     format_snow_params(sys.argv[1],sys.argv[2],sys.argv[3],sys.argv[4])  
> File "format_snow_params.py", line 112, in format_snow_params
>     bands = np.arange(minelv, maxelv+interval, interval) ValueError: arange: cannot compute length

Here is the modification of my question, following the comments. I used Luke's script, in which he added the print, and the output results is the following:

Min: --
Max: --
Int: 100
Traceback (most recent call last):
  File "format_snow_params1.py", line 194, in <module>
    main() 
  File "format_snow_params1.py", line 188, in main
    format_snow_params(sys.argv[1],sys.argv[2],sys.argv[3],sys.argv[4])
  File "format_snow_params1.py", line 100, in format_snow_params
    bands = np.arange(minelv, maxelv+interval, interval)
ValueError: arange: cannot compute length

the output results of the command gdalinfo -stats Nyando_basin_SRTMSnap.tif :

gdalinfo: Relink `/home/lgee/miniconda3/bin/../lib/libpng16.so.16' with `/lib/x86_64-linux-gnu/libpthread.so.0' for IFUNC symbol `longjmp'
Driver: GTiff/GeoTIFF
Files: Nyando_basin_SRTMSnap.tif
Size is 3300, 1500
Coordinate System is `'
Origin = (34.750000000000000,0.050000000000000)
Pixel Size = (0.000303030303030,-0.000333333333333)
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  34.7500000,   0.0500000) 
Lower Left  (  34.7500000,  -0.4500000) 
Upper Right (  35.7500000,   0.0500000) 
Lower Right (  35.7500000,  -0.4500000) 
Center      (  35.2500000,  -0.2000000) 
Band 1 Block=3300x1 Type=Float32, ColorInterp=Gray
  Minimum=1134.000, Maximum=2862.000, Mean=1735.713, StdDev=473.560
  NoData Value=-3.4028234663852886e+38
  Metadata:
    STATISTICS_MAXIMUM=2862
    STATISTICS_MEAN=1735.7128876922
    STATISTICS_MINIMUM=1134
    STATISTICS_STDDEV=473.55999424241

the output results of the command gdalinfo -stats Nyando_grid.tif :

gdalinfo: Relink `/home/lgee/miniconda3/bin/../lib/libpng16.so.16' with `/lib/x86_64-linux-gnu/libpthread.so.0' for IFUNC symbol `longjmp'
Driver: GTiff/GeoTIFF
Files: Nyando_grid.tif
Size is 20, 10
Coordinate System is `'
Origin = (34.750000000000000,0.050000000000000)
Pixel Size = (0.050000000000000,-0.050000000000000)
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  34.7500000,   0.0500000) 
Lower Left  (  34.7500000,  -0.4500000) 
Upper Right (  35.7500000,   0.0500000) 
Lower Right (  35.7500000,  -0.4500000) 
Center      (  35.2500000,  -0.2000000) 
Band 1 Block=20x10 Type=Byte, ColorInterp=Gray
  Minimum=1.000, Maximum=1.000, Mean=1.000, StdDev=0.000
  NoData Value=0
  Metadata:
    STATISTICS_MAXIMUM=1
    STATISTICS_MEAN=1
    STATISTICS_MINIMUM=1
    STATISTICS_STDDEV=0

attached script:

# import dependencies
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
warnings.simplefilter("ignore")

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
    try:

        # 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.ma.masked_where(elvhires<0,elvhires)

    # 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:
        os.remove(outSnow)    

    nbands = [] # blank list

    # try to write to output snow parameter file
    try:

        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]

                        # find min and max values for interval
                        minelv = np.min(tmp.astype(int)) - (np.min(tmp.astype(int))%interval)
                        maxelv = np.max(tmp.astype(int)) + (np.max(tmp.astype(int))%interval)

                        # create an array of band limits  
                        try:
                            bands = np.arange(minelv, maxelv+interval, interval)
                        except:
                            print('Min: %s\nMax: %s\nInt: %s'%(minelv, maxelv, interval))
                            raise

                        # create blank array for number of bands calculation...
                        if mask[i,j] == 1: # ...if it is not a masked pixel

                            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(cnt)) # write grid cell id

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

                                    f.write('{0:.4f}\t'.format(frac))

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

                                    except IndexError:
                                        muelv = 0

                                    f.write('{0:.4f}\t'.format(muelv))

                                # calculate the precipitation fractions and write to file
                                for c in range(maxbands):
                                    try:
                                        idx = np.where(bcls==uniqcnt[c])
                                        frac = np.float(idx[0].size) / np.float(bcls[np.where(bcls>=0)].size)
                                    except IndexError:
                                        frac = 0

                                    f.write('{0:.4f}\t'.format(frac))

                                f.write('\n') # write return value for new line

                        if pass_cnt == 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')

    return

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...")
        sys.exit()

    else:
        # Pass command line arguments into function
        format_snow_params(sys.argv[1],sys.argv[2],sys.argv[3],sys.argv[4])

    return  

# Execute the main level program if run as standalone
if __name__ == "__main__":
    main() 
  • 1
    Just as to debug the possible problem, could you add print(minelv, maxelv, interval, maxelv+interval) before bands = np.arange(minelv, maxelv+interval, interval) and see if they are properly created? – RJJoling Apr 9 '18 at 17:11
  • Thanks for your comment, I tried that but it's still the same problem – Bouramtane Tarik Apr 9 '18 at 17:20
  • 1
    That's to be expected, but what did the last line read (the print statement should have worked before the error occurred). – RJJoling Apr 9 '18 at 17:26
  • 1
    @BouramtaneTarik No one can help you if you can't simply read a line from your computer screen. – Jon Apr 9 '18 at 17:52
  • 1
    @BouramtaneTarik Adding the "print" statement will output the values of minelv, maxelv, interval, and maxelv+interval so you know what they are. It will do nothing to "fix" your code. It is a diagnostic tool. So you need to tell RJJoling what the values of those variables are. – Jon Apr 9 '18 at 19:42
1

Your basin mask has a different pixel size and a different shape (rows and cols) to your elevation grid. From what I can see, the basin mask should be the same size for the code to work.

Your Answer

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

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