Utilizing Numpy to perform complex GIS operation in ARCGIS 10

Today my country Nepal played with Jordan for securing  a place in world cup 2014 qualification. I was hoping for at least a draw but after 20 minutes or so, Nepal was already 3-0 down. I could not follow the live update anymore, therefore decided to grab a beer and  write something for my blog.This is how this post came into existence.

I wanted to write about array manipulation within a arcgis for a long time. The idea came from a problem that was faced by one of my close senior friend who  frequently uses  arcgis in his work. He had a raster (SRTM ) data and he wanted to calculated maximum value in a varying window ( like 5X5, 7X7, 9X9) for the data and marked that maximum value in the window as one. Focal statistics function in arcgis can calculated maximum value in a specified neighborhood but it assigns the maximum value in the centre cell. It does not tell where the maximum value lies within the window. So he asked me for suggestion, and I , being a MATLAB savvy person quickly wrote a function in MATLAB which would do the trick. I send the solution to the friend but the problem was, he was not familiar to MATLAB and he did not have access to MATLAB as well. He wanted to have a solution that works within the arcgis.

So, I have dived into NUMPY which is a fundamental package in Python for large, multi-dimensional array and matrices for solving the problem. It is already integrated into arcgis 10. In arcgis 10, there are two functions namely RasterToNumPyArray and NumPyArrayTORaster. First function converts a raster to matrix and second function converts a matrix to array. These two functions open a new opportunity for performing MATLAB like matrix manipulation within the arcgis.

So, in this post I will show you how problem mentioned before can be solved. First I will show you the code and will explain it concisely.

1. import arcpy


2.  from numpy import *


3. myArray = arcpy.RasterToNumPyArray('basic_raster')


4. [cols,rows]= myArray.shape


5.  winsize= 3;


6.  winhalf = 1


7.  I=zeros((cols, rows))


8.  for n in range(winhalf,rows-winhalf):


9. for m in range (winhalf,cols-winhalf):


10. data = myArray[m-winhalf:m+winhalf+1,n-winhalf+1:n+winhalf+1]


11. maxdata = data.max()


12. indices = where(data ==maxdata)


13.  I[m-winhalf+indices[0],n-winhalf+indices[1]] = 1

Line 1 and Line 2 imports the arcpy and numpy module.

Line 3 stores raster ‘basic_raster’ in a variable myArray

Line 4 stores size of the matrix myArray

Line 5 defines a window size ( three in this case)

Line 6 creates a matrix which is same size of rater ‘basic_raster’ but all filled with zero

Line 8-13 does the main computation. It will iterate through 3X3 window and calculates the maximum value and assigns 1 to the location of cell which contains the maximum value in the window.

Now we could already exports the I matrix using the function NumPyArrayTORaster (I). But doing so has couple of disadvantage because by doing so we would have just exported the matrix as a raster without any information about the cell size, coordinate and projection information. Therefore we need to get all these information from ‘basic_raster’ which is the raster based on which we calculated matrix I. For this we can use a bit of arcpy.

14 descData=arcpy.Describe('basic_raster')
15 cellSize=descData.meanCellHeight
16 extent=descData.Extent
17 spatialReference=descData.spatialReference
18 pnt=arcpy.Point(extent.XMin,extent.YMin)
19 newRaster = arcpy.NumPyArrayToRaster(I,pnt, cellSize,cellSize)

Line 14 uses a function ‘Describe’ to get information about the raster ‘basic_raster’.

Line 15-18 stores different information in various variables

Line 159transforms matrix I into a raster named ‘newraster’ using information from the raster ‘basic_raster’

Now we have a new raster where one represents cells with a maximum values in a sliding window of 3x3 which perfectly aligns with the raster ‘basic_raster’ in the same projection system.

I hope you have gained some insights into enormous advantage of integrating numpy and arcpy in your GIS workflow. If you have some prior programming knowledge, you would find pretty easy to learn python. Even if you don’t have any programming experience, I suggest you to gets your hands dirty with python. It’s not that hard to learn python, believe me :).

Utilizing Numpy to perform complex GIS operation in ARCGIS 10

Today my country Nepal played with Jordan for securing  a place in world cup 2014 qualification. I was hoping for at least a draw but after 20 minutes or so, Nepal was already 3-0 down. I could not follow the live update anymore, therefore decided to grab a beer and  write something for my blog.This is how this post came into existence.

I wanted to write about array manipulation within a arcgis for a long time. The idea came from a problem that was faced by one of my close senior friend who  frequently uses  arcgis in his work. He had a raster (SRTM ) data and he wanted to calculated maximum value in a varying window ( like 5X5, 7X7, 9X9) for the data and marked that maximum value in the window as one. Focal statistics function in arcgis can calculated maximum value in a specified neighborhood but it assigns the maximum value in the centre cell. It does not tell where the maximum value lies within the window. So he asked me for suggestion, and I , being a MATLAB savvy person quickly wrote a function in MATLAB which would do the trick. I send the solution to the friend but the problem was, he was not familiar to MATLAB and he did not have access to MATLAB as well. He wanted to have a solution that works within the arcgis.

So, I have dived into NUMPY which is a fundamental package in Python for large, multi-dimensional array and matrices for solving the problem. It is already integrated into arcgis 10. In arcgis 10, there are two functions namely RasterToNumPyArray and NumPyArrayTORaster. First function converts a raster to matrix and second function converts a matrix to array. These two functions open a new opportunity for performing MATLAB like matrix manipulation within the arcgis.

So, in this post I will show you how problem mentioned before can be solved. First I will show you the code and will explain it concisely.

1. import arcpy
2.  from numpy import *
3. myArray = arcpy.RasterToNumPyArray('basic_raster')
4. [cols,rows]= myArray.shape
5.  winsize= 3;
6.  winhalf = 1
7.  I=zeros((cols, rows))
8.  for n in range(winhalf,rows-winhalf):
9. for m in range (winhalf,cols-winhalf):
10. data = myArray[m-winhalf:m+winhalf+1,n-winhalf+1:n+winhalf+1]
11. maxdata = data.max()
12. indices = where(data ==maxdata)
13.  I[m-winhalf+indices[0],n-winhalf+indices[1]] = 1

Line 1 and Line 2 imports the arcpy and numpy module.

Line 3 stores raster ‘basic_raster’ in a variable myArray

Line 4 stores size of the matrix myArray

Line 5 defines a window size ( three in this case)

Line 6 creates a matrix which is same size of rater ‘basic_raster’ but all filled with zero

Line 8-13 does the main computation. It will iterate through 3X3 window and calculates the maximum value and assigns 1 to the location of cell which contains the maximum value in the window.

Now we could already exports the I matrix using the function NumPyArrayTORaster (I). But doing so has couple of disadvantage because by doing so we would have just exported the matrix as a raster without any information about the cell size, coordinate and projection information. Therefore we need to get all these information from ‘basic_raster’ which is the raster based on which we calculated matrix I. For this we can use a bit of arcpy.

14 descData=arcpy.Describe('basic_raster')
15 cellSize=descData.meanCellHeight
16 extent=descData.Extent
17 spatialReference=descData.spatialReference
18 pnt=arcpy.Point(extent.XMin,extent.YMin)
19 newRaster = arcpy.NumPyArrayToRaster(I,pnt, cellSize,cellSize)

Line 14 uses a function ‘Describe’ to get information about the raster ‘basic_raster’.

Line 15-18 stores different information in various variables

Line 159transforms matrix I into a raster named ‘newraster’ using information from the raster ‘basic_raster’

Now we have a new raster where one represents cells with a maximum values in a sliding window of 3x3 which perfectly aligns with the raster ‘basic_raster’ in the same projection system.

I hope you have gained some insights into enormous advantage of integrating numpy and arcpy in your GIS workflow. If you have some prior programming knowledge, you would find pretty easy to learn python. Even if you don’t have any programming experience, I suggest you to gets your hands dirty with python. It’s not that hard to learn python, believe me :).

ArcPy : Python scripting in ArcGIS 10

My hands were itching to get hold of newest version of ArcGIS for a really long. Finally, couple of days  back, i got ArcGIS in my aresnal. I am quite excited about it. ArcGIS 10 has introduced a lot more functionality for automation and analysis. The layout has now much elegant look compared to previous version where user can dock different window in sides as well as in bottom and its dock able. That means much more space to work.  The new tools such as split raster and raster mosaic are some of the tool that I wished ESRi had incorporated in its products much earlier.

It comes with numerous new and powerful tools such as ArcPY which provides an interface for a developer to program within the ArcGIS ! Till ArcGIS 9.3, writing python scripting in ArcPY required separate Python IDE. I have played around with ArcPY and i am amazed to see numerous operation it can perform and that too in just few lines of python scripting. Moreover it comes with intelligent window which helps you to complete code with instant help.

Today i will show you, how with simple and few lines, you can start doing geo processing with ArcPY.

The first line in picture imports complete modules which gives access to all the geoprocessing function and properties that are within that. This line is generally a starting line in any python scripting in ArcGIS. ArcGIS also includes other modules such as arcpy.sa for spatial analyst, arcpy.ga for geostastical analyst and arypy.mapping for mapping.

The second line specifies the workspace where we are going to store our features. The final line creates a buffer around the feature "schools", of 1000 feet and stored with a named with "Buffer_1000". It's easy like that!

You might say, I can do this thing with tools available within in arc catalogue toolbox easily with few clicks, so what this fuss about arcpy and why this arcpy would even needed! You are right! The true power of arcpy comes to play if you want to do a batch processing i.e. repeating same thing numbers of times.

In above case, if I wanted to carry buffer analysis with 1000 and 2000 feet then I have to carry the buffer analysis two times with same parameter only changing buffer distance. Things like this can easily be carried out with arcpy with a loop. The variable, "i" stores distance values which I want to buffer around the feature "schools". The variable "savenames" will generates string such as Buffer_1000 and Buffer_2000.

Here is the final output as viewed in arcgis 10. Notice the dockable windows in the sides of the picuture and python window at the bottom.

ArcPy : Python scripting in ArcGIS 10

My hands were itching to get hold of newest version of ArcGIS for a really long. Finally, couple of days  back, i got ArcGIS in my aresnal. I am quite excited about it. ArcGIS 10 has introduced a lot more functionality for automation and analysis. The layout has now much elegant look compared to previous version where user can dock different window in sides as well as in bottom and its dock able. That means much more space to work.  The new tools such as split raster and raster mosaic are some of the tool that I wished ESRi had incorporated in its products much earlier.

It comes with numerous new and powerful tools such as ArcPY which provides an interface for a developer to program within the ArcGIS ! Till ArcGIS 9.3, writing python scripting in ArcPY required separate Python IDE. I have played around with ArcPY and i am amazed to see numerous operation it can perform and that too in just few lines of python scripting. Moreover it comes with intelligent window which helps you to complete code with instant help.

Today i will show you, how with simple and few lines, you can start doing geo processing with ArcPY.

The first line in picture imports complete modules which gives access to all the geoprocessing function and properties that are within that. This line is generally a starting line in any python scripting in ArcGIS. ArcGIS also includes other modules such as arcpy.sa for spatial analyst, arcpy.ga for geostastical analyst and arypy.mapping for mapping.

The second line specifies the workspace where we are going to store our features. The final line creates a buffer around the feature "schools", of 1000 feet and stored with a named with "Buffer_1000". It's easy like that!

You might say, I can do this thing with tools available within in arc catalogue toolbox easily with few clicks, so what this fuss about arcpy and why this arcpy would even needed! You are right! The true power of arcpy comes to play if you want to do a batch processing i.e. repeating same thing numbers of times.

In above case, if I wanted to carry buffer analysis with 1000 and 2000 feet then I have to carry the buffer analysis two times with same parameter only changing buffer distance. Things like this can easily be carried out with arcpy with a loop. The variable, "i" stores distance values which I want to buffer around the feature "schools". The variable "savenames" will generates string such as Buffer_1000 and Buffer_2000.

Here is the final output as viewed in arcgis 10. Notice the dockable windows in the sides of the picuture and python window at the bottom.