1

I have computed Hydrologic Sensitivity Index (HSI) as a raster for a region using slope, Hydraulic conductivity, Depth to the Restrictive Layer and Flow accumulation. All the files are in rasters of same pixel size. The rasters are float type rasters.

Now, I would like to plot Hydraulic conductivity vs HSI, Depth vs HSI for every pixel etc.. In other words, for each pixel, I need the computed HSI value, Depth, Hydraulic conductivity, and Flow accumulation in a column or dataframe.

How can I accomplish this in ArcMap?

What I tried:

  1. First I tried to use the Combine tool. But Combine tool doesn't seem to work with float raster even if I convert it into integers there will be too many unique values to make this combination.

What if two different pixels had same values?

In that case, I do not think this tool will show those two pixel separately (right?).

  1. Convert each Raster to Int Raster than to Polygon than use Union tool. I couldn't accomplish this because the HSI raster file is very detailed and takes too long to convert to polygon and seems to crash in the middle.

What I think could work:

  1. I know there is a way to convert Raster data into a Numpytable. This way I will have four separate Numpy table for each of the four data but how will I map the value in such a way that i ensure that they are for the same pixel?

I know I can use arcpy.Describe.polygon.OIDFieldName to create id in Polygon but is there similar tool for Raster. But i am also not sure when I am using this tool for four different raster.

Will the pixel in the same location get the same oid values?

Improve this question
8
  • If they all have the same dimensions to begin with then try arcpy.NumpyArrayToRaster() 4 times, np.dstack() to create a (y_size, x_size, 4)-shaped array, then stack.reshape((-1, 4)) to achieve your 4 columns
    – mikewatt
    Apr 19, 2021 at 23:51
  • @mikewatt Thank you for this suggestion. They all have the same pixel size. when i do the stacking using np.dstack() how does the function know that it is maping the values for the same pixel without us providing any id to compare to (like we do while doing a polygon join)?
    – Samrat
    Apr 20, 2021 at 0:02
  • They all have to have the same width, height, pixel size, and geographic location for this to work. But if they do, then the arrays that you load will all have the same shape and the order of pixels will be consistent across all of them. Any numpy operations like dstack() will always preserve the order as appropriate (see the examples)
    – mikewatt
    Apr 20, 2021 at 0:16
  • If you really after entire raster (not just streams) compute and combine rows and columns using gis.stackexchange.com/questions/204037/… You can use Sample to collect as many rasters values as you want.
    – FelixIP
    Apr 20, 2021 at 3:44
  • @mikewatt When i try your suggestion i run into two problems: 1. when i use arcpy.RasterToNumpyArray() i get extremely large negative values in my array. I checked my Raster file and it seems to be fine. Also, all three raster shows exactly the same value and all large negative number. 2.i ignored the large values and went to the next step. I was able to stack the values from three raster together but i couldnt convert the array into gistable. Is there a easy way to export the final array out of arcpy environment? Any suggestions?
    – Samrat
    Apr 21, 2021 at 23:09

2 Answers 2

Reset to default
1

Set environment extent and cell size to your raster. Make sure workspace is also specified. In Python window type:

arcpy.gp.SingleOutputMapAlgebra_sa("$$ROWMAP","nRow")

This will produce integer raster of rows:

Similarly compute raster of columns:

arcpy.gp.SingleOutputMapAlgebra_sa("$$COLMAP","nCol")

enter image description here

In Arcmap Options increase number of unique numbers to render to nr*nc, where nr is number of columns in you raster, nc -columns. Use Combine tool from spatial analyst to combine nRow and nCol into single raster where each cell is unique:

enter image description here

Use Sample to sample your rasters, e.g.

arcpy.gp.Sample_sa("FACC;FDIR", "Combine", "in_memory/abc", "NEAREST", "Value")

Output table will have all you need, even coordinates of cell centres:

enter image description here

Improve this answer
11
  • This worked so well. Thank you!
    – Samrat
    Apr 24, 2021 at 1:05
  • Perhaps mark it as solved, please.
    – FelixIP
    Apr 24, 2021 at 1:31
  • i tried to upvote it last time but it said i dont have enought "reputation" to upvote a post.. I didnt know there was a "Solved" button. I checked a mark in your answer. Thanks again!
    – Samrat
    Apr 25, 2021 at 21:36
  • Ok, cheers mate
    – FelixIP
    Apr 26, 2021 at 0:25
  • sorry to bug you again. Previously i wrote a code for a small subarea and the process worked fine. But when i do the same thng on whole area, I am getting following error after creating nrows and ncols, and increasing render limit (78403059) while using Combine tools: Runtime error Traceback (most recent call last): File "<string>", line 1, in <module> File "c:\program files (x86)\arcgis\desktop10.7\arcpy\arcpy\geoprocessing_base.py", line 510, in <lambda> return lambda *args: val(*gp_fixargs(args, True)) ExecuteError: ERROR 001143: Background server threw an exception.
    – Samrat
    Apr 26, 2021 at 18:26
1

A variant on @FelixIP's approach is first take one of your rasters and run it through the raster to point tool. This creates a point for every pixel. You now have a point dataset which you could apply a selection to, e.g. a particular catchment, then run it through the Sample tool as he describes above.

Improve this answer
2
  • Thank you. Your explanation helped.
    – Samrat
    Apr 24, 2021 at 1:06
  • It is an option, but much slower, than one with rasters only, I guess it is so, because no shapes is created in the process.
    – FelixIP
    Apr 26, 2021 at 3:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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