Encontrar ubicaciones de trama máximo dentro de los polígonos: Delimitación de la Cuenca

Question

Estoy tratando de identificar la Ubicación (lat/lon) y el área de upstream de cada punto de trama que desemboca en el Océano Ártico (hay un montón, así que tiene que ser totalmente automatizado). Estoy empezando con el siguiente rásteres:

• La Dirección Del Flujo
• Accumlated Área De Flujo

Mi enfoque ha sido hasta ahora el uso de la por encima de los rásteres para crear un delineado de la cuenca del shapefile. Desde allí tenía la esperanza de que podría utilizar el archivo de forma que buscar en el Flujo de Acumulación de Trama y de identificar la ubicación de el valor máximo en cada cuenca. En términos de producción, sólo tengo una tabla que muestra [Cuenca#, Longitud, Latitud, UpstreamArea] para cada cuenca.

Todas las ideas o sugerencias sería apreciado. [Tenga en cuenta que ya he mirado aquí , pero no sé cómo restringir la trama solo a la zona delimitada por los polígonos en el shapefile]

*Uso: ArcGIS10, Python, etc.

Gracias,

Joe

Answer 1

Supongo que sólo hay 1 píxel en cada polígono que tiene un valor máximo.

1. Encontrar El Valor Máximo

... utilizar el archivo de forma que buscar en el Flujo de Acumulación de Trama y identificar la ubicación de el valor máximo en cada cuenca.

Para resolver esta tarea (de encontrar el máximo valor de píxel en cada polígono) puede utilizar Zonal de Estadísticas (Spatial Analyst).

in_zone_data = polígono archivo

zone_field = cuenca número de

in_value_raster = flujo de acumulación de trama

statistics_type = MÁXIMO

2. Convertir la cuenca del polígono a trama

valor de píxel = cuenca número de

3. Encuentre la ubicación de el valor máximo

Uso Con (Spatial Analyst) para encontrar la ubicación del máximo en cada cuenca.

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

El resultado de trama contiene un píxel para cada cuenca. El valor del píxel muestra la cuenca del número.

4. Convertir la trama a punto de

convertir el resultado de la trama del paso 3 al archivo de punto (sólo para los valores de los píxeles no 0)

5. calcular lat/lon para el resultado del paso 4

Answer 2

He aquí lo que se me ocurrió después de probar un número de enfoques, utilizando herramientas de ArcMap. Mi python habilidades están en su infancia, pero funciona (a pesar de que tomó cerca de 36 horas para recorrer los 110.000 cuencas.

1. Crear lat/lon matrices que corresponden a la cuenca y el área de flujo de archivos.
2. Un bucle en cada cuenca, encontrar el valor máximo y la ubicación de los índices en cada uno de los
3. Encontrar la lat/lon en cada ubicación del índice.

4. Escribir.

   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[i]=np.int(j)
       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:
       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=',')