Python tutorial: Converting a raster dataset to XYZ in Python

One of the most popular posts in my blog is about converting a raster image into XYZ text file.  Converting raster image to XYZ file may be necessary because machine learning algorithms (outside proprietary software) requires input to be a table. I had written a post about converting a raster file into XYZ using ARCGIS some three years ago, which seems to attract many visitors to my blog.

Here I will write a way to do it using Python. For reading geo-referenced raster file, we will use rasterio package which is a wrapper for gdal that provides clean and fast I/O for geospatial raster images. The package is written in cython therefore it is very fast. Reading a raster file with rasterio is a one liner code. You can download binary of rasterio here. Its a binary file so installing rasterio is just a matter of clicking binary executable file. After reading the raster file we will get bounding box of the image and compute XY of each pixel and later write in a csv file together with pixel values. The code works for any number of bands in the image.

I hope this piece of code is useful for you. In near future I will show you to do exact same thing within QGIS. Stay tuned.

Output XYZ file with band values

In addition, the above code can be combine with the code (Data exchange between MATLAB and Python: Reading and writing .mat files with Python)  so that  you can easily export multi-spectral and hyperspectral data as a mat file for MATLAB. Many people seem to be having problem with reading remote sensing image with multiband read function in MATLAB. If you are one them, use above codes and bypass multiband read function to get your remotes images straight in MATLAB.

Data exchange between MATLAB and Python: Reading and writing .mat fileswith Python

I am a MATLAB guy. I love MATLAB. I started with Fortan during my first degree at my Bachelor level. I came across with MATLAB only in year 2010 as a part of my Msc studies. Since then I been using it continuously for various tasks. Just recently, I started  using Python and I LOVE IT. The biggest reason being that Python is open-source and there are lots of Python libraries that can be used. But there are times; I have to shuffle data between MATLAB and python, because for some tasks I prefer MATLAB as I am accustom to it. So if you are in the same boat, then here I show you simple way to transfer data between MATLAB and Python. One has to use scipy.io module in Scipy for that purpose. If you have some old data or got some data online that are saved as MATLAB’s .mat file format, you can simply import it as:

import numpy as np
import scipy.io as sio
mydata = sio.loadmat('mydata.mat')

Now your mydata contains a dictinary with keys corresponding to the varible names saved in the original mydata.mat Saving a python variable to .mat file is also straight forward

# write one variable
x = np.arange(1,10,1)
# file name of the file
fname ='export_from_python.mat'
sio.savemat(fname, {'x':x})

When you read export_from_python.mat, you will get 'x' varialbe into MATLAB. If you want to write more than two variable:

# write two variables
x = np.arange(1,10,1)
y = np.ones((5,5))
fname ='export_from_python_1.mat'
sio.savemat(fname, {'x':x, 'y':y})

I do lot of work with classification of remote sensing images using many different types of machine learning algorithms. Proprietary software like ENVI and ERDAS are not flexible enough for me, as these software donot allow an efficient way to tune hyper-parameters that are algorithm specific. I dont like training models with defualt parameters. I do my model traning and classification in python. Here is the typical workflow.

1.  Import image and training data in python using GDAL
2. Train machine learning model in python using Sci-kit Learn
3. For accuracy assessment, i do it within python but for accuracy assessment visualization, i export  targets and outputs vectors from Python to MATLAB
4. Use 'plotconfusion' function from Neural Network toolbox for a visualization of confusion matrix. Here is the confusion Matrix  information that i got from MATLAB. Its nice, isnt it ?

Confusion matrix

Example of extracting grid nodes in MATLAB

For last few posts, I am writing about morphological image analysis with python with open source package scikit-image. In this post I want to show benefit of morphological image analysis with some real application. I been working with road extraction from satellite images for past few weeks and  one of the common task is to extract the grid corner nodes of the road network. So let’s assume the road grid has been already extracted and we want to extract intersection of grid network. Here I am showing three approaches for that based on synthetic data. Obviously with real data, the task is more complex and might require many slight modifications.

Here is the grid. You can generate such grid with checkerboard function in MATLAB with minor processing (edge detection and hole filling). Intersections of of grid lines have a property that it has 4 white pixels in 4-N neighborhood.  We can exploit this property to extract it. I am going to show to ways by which it can be done and in addition, another method based on extraction of corners will also be shown.

%make a checkerboard
I=checkerboard(20, 10,10)>0.5;
imshow(I)
% detect edge with canny
BW = edge(I,'canny');
figure, imshow(BW)
%make square structural element to fill holes with closing
SE = strel('square', 3);
BW1 = imclose(BW, SE);
figure, imshow(BW1)
% make a skeleton of the to lines 1 pixel thick
BW2 = bwmorph(BW1,'skel',Inf);
figure, imshow(BW2)

 


In mathematics and, in particular, functional analysis, convolution is a mathematical operation on two functions f and g, producing a third function that is typically viewed as a modified version of one of the original function .In simple word, convolution of image is a process where you scan your image from left to right and top to bottom within local neighborhood defined by user, and do mathematical calculation within the local neighborhood. Here we are going to use it for simply sum the number of pixels which are which are white within 3x3 neighborhood. As shown by above figure, if the point is an intersection grid point then at that point sum should be 5. Any points with less than 5 white pixels are not intersection grid points.

% with block processing and inline function with convolution
tic;
kernel = [0 1 0; ... %# Convolution kernel
1 1 1; ...
0 1 0];
sumX = conv2(double(BW2),kernel,'same');
result=sumX;
% only consider pixels which are in grid image
result (BW2==0)= 0
result(sumX<5)=0;
[r,c] = find(result>0);
toc;
figure,imshow((BW2)), hold on, plot(r,c,'r*');title ('with Convolution')
hold off

Way Two: With Mathematical Morphology (MM)

If you are working in the field of Geo-information and you still saying what the heck this ‘MM’ then seriously you should be acquainted with this field of image processing. You can use MM for classification, edge detection, filtering, segmentation, building detection and many other things. Grab a cup of coffee and Google MM in remote sensing; there is tons of stuff to read.  There are courses in MM taught in different universities which are one semester long, so you got an idea how vast the subject is. If you want to seriously know MM, then there is no better book than ‘ Morphological Image Analysis’ by P Soillie. Here we use ‘Erosion’ method of MM.  Time: 0.011184 seconds.

% method two with morphological analysis
tic;
SE = strel ('disk', 1);
result2 = imerode(BW, SE);
[r,c] = find(result>0);
toc;
figure,imshow((BW2)), hold on, plot(r,c,'b*'), title ('With MM Erosion')
hold off

Way Two:  Mathematical Morphology (MM)

There are many corner detection techniques such as Moravec, Harris, SIFT etc. Here, I am going to use Harris corner detector for intersection grid points. Time: 0.094848 seconds.

% detect corner by harris corner detector
tic;
C = corner(BW2,500 );
timet2 = toc;
toc;
figure,imshow(BW2);
hold on
plot(C(:,1), C(:,2), 'g*'); title('with harris corner detector')
hold off

Way Two:  With Harris Corner Detector

So you have seen there are number of ways to solve a particular problem. All above methods detected intersection points successfully. Among three methods, convolution required longest time for processing and MM required shortest time. The gain in computation cost between MM and Corner detection is 9 folds.  I did all this in MATLAB but it can easily coded in python with scikit-image.

A little teaser for using scikit learn for classification of remotesensing images

I work with classification of remote sensing images a lot both with supervised as well as unsupervised classification. Unsupervised classifications don’t need any external input where as supervised classifications need samples or training areas for an algorithm to learn. For processing remote sensing images, there are many proprietary software like ENVI, ERDAS, PCI Geomatica, Global Mapper, eCognition and many more. Even after paying thousands of Euros, classifications algorithm available in Costs-off-the-self (COTS) software is far from satisfactory.

Software	Available Algorithms
eCognition 8.7 (cost >1000 Euros/year)	KNN Decision Tree DT SVM (no Y and C parameters available for tunning, no way to perform grid search for optimal determination of Y and C, Only linear and rbf kernel available) Random Forest  RF (only available in 8.8
ENVI (cost >1000 Euros/year)	Maximum Likelihood (ML) SVM (no way to perform grid search for optima, unless you get your hands dirty with IDL programming Neural Network NN ( Numbers of hidden nodes cannot be assigned) Some otheralgorithms more suitable for Hyperspectral imagesy SAM, SID, Spectral Unmixing.
Scikit learn Free Opensoure	Neighrest Neigbour (NN) Decision Tree (DT) SVM ( Grid search, cross validation flexibility) Random Forest RF AdaBoost Naives Bayes Linear Discriminant Analysis (LDA) Quadratic Discriminant Analysis (QDA)

Here is the little teaser of classification accuracy with many algorithms that are available in scikit-learn for a remote sensing imagery. In near future, I will blog with more illustration and with code. Till then go and make your hands dirty with Python and Scikit-Learn. Make that your new year resolution and trust me, you will thank me for that.



Here, algorithms hyperparameters were not optimally tuned hence superior machine learning algorithms like SVM has very low accuracy for test samples which are not seen by trained model.

Image processing in python with scikit-learn and pymorph

I have been using MATLAB quite a bit for some time. But being a proprietary language it has some drawbacks. One has to fork out lot of money to buy it and renew it yearly license. On top of that one hast to busy many toolboxes.Check this blog for finding various advantages of switching over from MATLAB to Python. You will find very helpful insights into key difference  between MATLAB and python numpy here.

Recently for couples of weeks now, I have to experimenting with python which is an open source software for some of my tasks that I do mainly related to classification, segmentation of remote sensing image etc during my weekend free time. So after sacrificing two of my weekends time for scouring python and different libraries which otherwise I would spend for my leisure activity, I am now pretty comfortable with python language as well for image processing tasks. Sweet :-) .

In this post I will show some morphological image analysis with python. I can see this will be series of post related to python. I plan to cover classification, clustering, edge detection in my subsequent posts.

For installing python, I have used pythonXY which comes with many libraries that might be needed for the scientific data analysis such as scipy, numpy, matplotlib and many others.PythoXY also comes with an Interactive development environment (IDE) for python which is called Spyder. With Spyder, you don’t have run python from command prompt. Other benefit includes integrated help, variable explorer, history console etc. In addition its all FREE.

Here is the interface how it looks.

I will be using library called “scikit-image” which comes in a bundle with pythonXY and in addition I would be using a library called “PyMorph” for morphological image analysis. Scikit-image also has a morphologicl module but there you would only find basic morphological operators like opening, closing, erosion and dialation. In PyMorph there are many advance morphological operator such opening by reconstruction and closing by reconstruction, Auto sequential filter(ASF), ASF by reconstruction etc. So after installing pythonXY, download the PYMoprh library and install it.

So, in this post I would show you how to use scikit-image for creating some basic shapes, calculating i area for each shapes or blobs and writing areas in the figure, finding bounding box of each blobs and plotting it in the figures. After finding areas of each figures, I will show how to remove certain blob  which has less area t than specified area.

% -*- coding: utf-8 -*-
Created on Mon Nov 18 15:42:21 2013
@author: shailesh
"""
%%%
% call necessary libraries
import math
import matplotlib.pyplot as plt
import numpy as np
from skimage.draw import ellipse
from skimage.draw import polygon
from skimage.draw import circle
from skimage.morphology import label
from skimage.measure import regionprops
from skimage.transform import rotate
import matplotlib.patches as mpatches
import pymorph as MM

%%%
% draw some arbitary shapes
image = np.zeros((1000, 1000),dtype=uint8)
% create an ellipse with centre (350,350) with minor axis = 100 and major = 220
rr, cc = ellipse(300, 350, 100, 220)
image[rr,cc] = 1

% create a polygon
x = np.array([1, 7, 4, 1])
x = np.array([1, 70, 40, 1])
y = np.array([1, 20, 80, 1])
rr, cc = polygon(y, x)
image[rr, cc] = 1

% create a polygon
rr, cc = circle(200, 200, 50)
image[rr, cc] = 1
fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(6, 6))
ax.imshow(image)

%lable connected regions
label_img = label(image)
%find properties of connected regions
regions = regionprops(label_img,['Centroid', 'BoundingBox', 'Area'])

%loop through connected regions and find
% Centroid', 'BoundingBox', 'Area' for each blobs
for props in regions:
y0, x0 = props['Centroid']
minr, minc, maxr, maxc = props['BoundingBox']
area =props['Area']
rect = mpatches.Rectangle((minc, minr), maxc - minc, maxr - minr,
fill=False, edgecolor='red', linewidth=2)
ax.add_patch(rect)
ax.text(x0, y0, str(area), fontsize=10, color="blue")
%show figure
plt.gray()
plt.axis((0, 1000, 1000, 0))
plt.show()

%%% remove small blobs
image = image.astype(np.uint8)
%remove blobs with area less than 10000 '
% MM.areaopen is coming from PyMorph
b=MM.areaopen(image,10000)
fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(6, 6))
ax.imshow(b)
ax.title('blobs with areas greater than 10000 ')
plt.gray()
plt.axis((0, 1000, 1000, 0))
plt.show()
% to close the pixels
plt.close("all")

Here are the figures produced by above code.

Initial figure after some blobs were made. The area of each blobs and bounding box were plotted as well


Figure where blobs with area less than 10000 pixels are remove with pymorph.


So this is a hypothetical simple illustration. But the real application of this kind of analysis could be many in remote sensing field. For example, after detection buildings you can easily remove buildings with size less than minimum mapping unit (say 100 m2). Similarly it can be used for removing some artifacts which are small in area that you might not be interested in. And the list can go on.

In next few blogs, i will be writing more about other morphological image processing operators as well as classfiers such Support Vector Machine (SVM), Random Forest (RF) etc.

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.