I'm trying to identify the Location (lat/lon) and upstream area of every raster point that flows into the Arctic Ocean (there are a lot so it has to be fully automated). I'm starting with the following rasters:

  • Flow Direction
  • Accumlated Flow Area

My approach so far has been to use the above rasters to create a delineated basin shapefile. From there I was hoping I could use the shapefile to look at the Flow Accumulation Raster and identify the location of the maximum value in each basin. In terms of output, I just need a table that shows [Basin#, Longitude, Latitude, UpstreamArea] for every basin.

I've already looked at How to find the locations of the highest values in a raster? but I don't know how to restrict the raster to only the area bounded by the polygons in the shapefile]

**Using: ArcGIS10, Python, etc.*

  • For those familiar with the hydrology tools, I'm essentially trying to find exactly 1 pour point for every basin, regardless of its size. If I use the basin shapefile the input to Snap Pour Point, I get a pour point for every raster grid cell. Hence why I'm apraching this outside of the Hydrology tools.
    – sea_hydro
    Commented Nov 21, 2012 at 16:01

2 Answers 2


I assume that there is only 1 pixel per polygon that has a maximum value.

1. Find Maximum Value

... use the shapefile to look at the Flow Accumulation Raster and identify the location of the maximum value in each basin.

To solve this task (find the maximum pixel value in each polygon) you can use Zonal Statistics (Spatial Analyst).

in_zone_data = your polygon file

zone_field = basin number

in_value_raster = flow accumulation raster

statistics_type = MAXIMUM

2. Convert basin polygon to raster

pixel value = basin number

3. Find location of the maximum value

Use Con (Spatial Analyst) to find the location of the maximum in each basin.

con ( <raster result from step 1> == <raster flow accumulation>, <raster from step 2>, 0)

The result raster contains one pixel for each basin. The pixel value shows the basin number.

4. Convert raster to point

convert result raster from step 3 to point file (only for pixel values not 0)

5. calculate lat/lon for result from step 4

  • Thanks for the input. I also need to find the raster grid location of that maximum which is why I haven't been able Zonal Statistics. Also, the shapefile is quite large so it would be best if I could use a scripting aproach where I can loop over the basin polygons.
    – sea_hydro
    Commented Nov 21, 2012 at 15:57
  • You're about one-eighth of the way there: you have found the value of the maximum. To obtain the location, you need to go further. There are several approaches. One way is to create a zonal maximum raster--not a table--and compare that to the original DEM. This gives a 0-1 grid showing where the maxima occur. Use this to pick out the values in a grid of x coordinates, another grid of y coordinates, the flow accumulation grid, and a grid of basin IDs. Output these values as a table.
    – whuber
    Commented Nov 21, 2012 at 16:38
  • @sea_hydro: I edited my first answer while whuber makes a comment with an other aproach.
    – Jens
    Commented Nov 21, 2012 at 16:55

Here's what I came up with after trying a number approaches using ArcMap tools. My python skills are in their infancy but it works (although it took about 36 hours to loop over the 110,000 basins.

  1. Create lat/lon arrays that correspond to the basin and flow area files.
  2. Loop over each basin, finding the maximum value and location indices in each
  3. Find the lat/lon at each index location.
  4. Write out.

    import numpy as np
    from osgeo import gdal
    #Load input Rasters (Basin Mask and Accumulated Upstream Area)
    #Input rasters need to be the same size
    f1 = gdal.Open("/pour_points/basins.asc")
    f2 = gdal.Open("/pour_points/source_area.asc")
    basins = np.array(f1.GetRasterBand(1).ReadAsArray())
    areas = np.array(f2.GetRasterBand(1).ReadAsArray())
    basins = np.ma.masked_values(basins,-9999)
    areas = np.ma.masked_values(areas,-9999)
    #Extract raster information (dimensions, min,max, etc)
    width = f1.RasterXSize
    height = f1.RasterYSize
    gt = f1.GetGeoTransform()
    res = gt[1]
    minx = gt[0]
    miny = gt[3] + width*gt[4] + height*gt[5]
    maxx = gt[0] + width*gt[1] + height*gt[2]
    maxy = gt[3]
    #Make arrays of lats and lons
    lons = np.linspace(minx,maxx-res,num=width)
    lats = np.linspace(maxy-res,miny,num=height)
    #Make arrays of same dimensions as input arrays of lat/lon values
    x,y = np.meshgrid(lons,lats)
    #Setup basin in/out arrays
    basin_ids = np.arange(np.min(basins),np.max(basins))
    #Test basins = basin_ids = np.array([57087,53728,55400,37250,60966,71339])
    basin = np.zeros(len(basin_ids),dtype='i')
    max_area =  np.zeros(len(basin_ids),dtype='i')
    x_outlet = np.zeros(len(basin_ids),dtype='f')
    y_outlet = np.zeros(len(basin_ids),dtype='f')
    min_x = np.zeros(len(basin_ids),dtype='f')
    max_x = np.zeros(len(basin_ids),dtype='f')
    min_y = np.zeros(len(basin_ids),dtype='f')
    max_y = np.zeros(len(basin_ids),dtype='f')
    #Loop over every basin id, finding the maximum upstream area [and location]
    #and record the basin#,longitude,latitude,area
    for i,j in enumerate(basin_ids):
        #print 'i = ',i
        x_basin = x[np.nonzero(basins==j)]
        y_basin = y[np.nonzero(basins==j)]
        max_area[i] = np.int(max(areas[np.nonzero(basins==j)]))
        max_ind = np.argmax(areas[np.nonzero(basins==j)])
        x_outlet[i] = x_basin[max_ind]
        y_outlet[i] = y_basin[max_ind]
        min_x[i] = min(x_basin)
        max_x[i] = max(x_basin)+res
        min_y[i] = min(y_basin)
        max_y[i] = max(y_basin)+res
    #save the list of pour points as a comma seperated text file
    # This file is directly importable into arcgis for validation purposes
    out_file = 'arctic_pour_points.txt'
    with file(out_file,'w') as outfile:

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.