Find locations of raster maximum within polygons: Basin Delineation

[JoeHamman]

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


•     Accumulated 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.

Any ideas or tips would be appreciated. [Note that I've already looked here but I don't know how to restrict the raster to only the area bounded by a polygons in the shapefile]

*Using: ArcGIS10, Python, etc.

Thanks,

Joe

[MarkBoucher]

Joe,

What do you mean by, "the maximum value in each basin"? Do you mean the max value for the flow accumulation grid (fac)?

The max value fac for each basin will be at the coast. Not all of the grids will flow to a stream that flows to the coast. I would think that you will have many very small watersheds, possibly several that are one grid in size, that drain to the ocean. So, in theory, every grid that boarders the ocean will be a point of interest for you.

If you can isolate the grid cells in the fac that boarder the ocean, you can convert them to points using tools. Then you could filter them based on some minimum fac value. Then you can add lat & long fields to them and calculate geometry for the lat & long.

[JoeHamman · Accepted solution]

MBoucher21:  Thanks for the input.  There are alot of basins that only have 1 cell.  Thats ok, I need to know their location too.  We can't assume that the maximum value happens to be at the coast.  Infact, it often is not as their are deltas and islands that often show up as the "coast".

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             basin=np.int(j)             x_basin = x[np.nonzero(basins==j)]             y_basin = y[np.nonzero(basins==j)]             max_area = np.int(max(areas[np.nonzero(basins==j)]))             max_ind = np.argmax(areas[np.nonzero(basins==j)])             x_outlet = x_basin[max_ind]             y_outlet = y_basin[max_ind]             min_x = min(x_basin)             max_x = max(x_basin)+res             min_y = min(y_basin)             max_y = 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:             outfile.write('OID,longitude,latitude,basin_area,min_lon,min_lat,max_lon,max_lat\n')             np.savetxt(outfile,(np.array([basin,x_outlet,y_outlet,max_area,min_x,min_y,max_x,max_y]).T),delimiter=',')