Publishing an ArcGIS Geoprocessing Service which uses Python

The Scenario:

I have a Python script which I wish to use within a geoprocessing tool. I want to publish this tool as a geoprocessing service using ArcGIS Server.


The Problem:

In the past I've written scripts, only to struggle integrating it into a production environment. In the ArcGIS environment, this often occurs when I've finished a Model or Python Script, which works fine in ArcGIS Desktop, but not as a published Geoprocessing Service.


The Solution:

The main thing to understand is where these different technologies interact (eg. Python, ModelBuilder, ArcGIS Server). It is these interfaces between your script, your model, your toolbox and your data that must be understood.

Before I walk through the process of publishing a Python script as a GeoProcessing Service, I will highlight a few important things to consider:

Make sure all required resources are available to the ArcGISSOC user from the ArcGIS Server instance. This includes your data, SDE connection files, your scripts, network resources, etc. Anything that your script and/or model depends on must be accessible by the ArcGISSOC user from the ArcGIS Server. This can be achieved by either:
A) Sharing the resources on your network, eg. by creating a single directory for all resources; or
B) Copying the required resources to the server and updating links within the scripts and models (annoying at best).


1) Create the Python Script
Make sure you set the return parameter using the SetParameterAsText(index,value) method.

Heres an example Python script:

import arcgisscripting
import sys
import random

gp=arcgisscripting.create()

num_points=gp.GetParameterAsText(0)
ws=r"%in_memory%"
gp.Workspace=ws
myFCName=r"random_points"

gp.CreateFeatureClass_management(gp.Workspace,myFCName,"POINT")
gp.AddMessage("Feature class created")

fcDesc=gp.Describe(gp.Workspace+"/"+myFCName)

cur=gp.InsertCursor(myFCName)
pnt=gp.CreateObject("Point")

for i in range(num_points):

feature=cur.NewRow()
    pnt.y=random.random()*180 - 90
    pnt.x=random.random()*360 - 180
    feature.shape=pnt
    cur.InsertRow(feature)

del cur

myFeatureLayer=gp.MakeFeatureLayer_management(myFCFullPath,myFLName)
gp.SetParameterAsText(1,myFeatureLayer)

2) Import the script into an ArcToolbox
As far as I can tell, this process just registers a Python script with a toolbox and wraps it in a ModelBuilder Tool template. It allows the script to behave like a Tool in any other toolbox. As part of this process, you will need to specify the location of the script. This is why your script must be accessible by the ArcGISSOC user as described above. You must also, specify the input and output parameters that your script is expecting and their data types.

Example values for an output parameter:
Display Name: BufferedFeatures
Data Type: Feature Layer
Type: Derived
Direction: Output

3) Create a Geoprocessing Model using ModelBuilder.
The main thing to remember here is to use the correct input and output parameter data types. There are limitations on which data types can be exposed via the REST interface. These restrictions only apply to the input and output values of the MODEL, not the Script! Details can be found here

4) Add the registered script into the model
Drag the new Script from your Toolbox (in ArcCatalog) into the new Model. This will add the Python script into your Model and will behave like any other tool. Link the input and output parameters to the script, by opening the script within your model and selecting the parameters you created earlier.

5) Publish the Toolbox
Your toolbox can now be published to ArcGIS Server. Simply right click on the Toolbox in ArcCatalog and select "Publish to ArcGIS Server...". Select next to the prompt.


James.

Evolution of the Human Processor

More and more frequently, I hear people talking about the importance of processing huge datasets. These grand ideas often push the limitations of technology. I can understand the urge to pursue these ideas, but I feel that the mentality is misguided. Mother nature worked out a long time ago that processing data in this way wastes precious neurons. Like us, nature faces a limitation of resources, both in storage capacity and processing power. I've heard (admittedly, from a dubious source) that the human brain has a capacity of about 4 terabytes. This may sound like a lot, but it is not when you consider the storage requirement for all a video and audio over your whole lifetime. So how does the human mind deal with this flood of information?

By prioritising knowledge. Think about the structure and operation of the human eye. We have central vision, with which we continuously change focus between things that we consciously or subconsciously decide is important. We also have peripheral vision which allows us to detect important information within our current view. Between these two extremities, we have a gradual transition between very high resolution and very low resolution.

What we don't have is: a uniform grid of retinal receptors to process data from all directions equally. Unlike a bug, we have evolved to assign attention to very small subset of the information that surrounds us, and then to seek out further detail.

Over millions of years of evolution, nature has determined that all information should not be processed equally. It does this by extracting a large amount of information from a small field of view (think hunting), and only processing small amounts from our peripheral vision (sudden movements).

So how can we apply these observations to real-world problems?

This philosophy should be applied to any high-volume data processing scenario. Just remember a few concepts:

1) Devise rules for deciding what characterises valuable information
2) Formulate efficient methods to extract the valuable information.

In my next post I will discuss techniques for improving algorithm efficiency.

James.

Displaced REST Map Service in Flex Application

The Scenario:
Recently, I've been doing some development with the ArcGIS API for Flex and ESRI's open-source "Flex Viewer" application. The application includes several widgets which work out-of-the-box. After configuring an application with the basic widgets, I got to the point where I needed extra functionality. In this particular situation, I wanted to be able to load external services (eg. ArcGIS REST or WMS) at runtime. I decided that the most appropriate approach would be to extend the "LayerList" widget.

The Problem:
The modification appeared to be working fine on the development environment. I could load several layers, both REST and WMS, which displayed beautifully.
However, when the flex app was migrated to the production system, I started to notice some issues. The application was fired up, and I loaded an external REST service. It seemed to be working, but when I maximised the window, the data no longer aligned with the basemap. At first I thought that my code had modified the spatial reference for the layer.

The Solution:
After some investigation it occurred to me that the top left corner of the image was in the correct location. I also realised that the returned image was over the correct extent. It turned out that the image which I was requesting was larger than the maximum image size configured for this service. (The default for ArcGIS REST services is 2048x2048 pixels.)

The solution is to modify the service configuration file. On windows servers: C:\arcgis\server\user\cfg\service_directory\service_name.cfg
Simply update the values for MaxImageHeight and MaxImageWidth. I set both of mine to 3072, like so:

<ServerObjectConfiguration>
<Properties>
...
<MaxImageHeight>3072</MaxImageHeight>
<MaxImageWidth>3072</MaxImageWidth>
...
</Properties>
</ServerObjectConfiguration>

I encountered a similar situation with WMS Services, but unlike the REST services, the limited image size (2048x2048) was positioned correctly.
I am assuming that the WMS services are processed more easily because the maximum file size information is visible to the client in the GetCapabilities file. ArcGIS REST services do not expose this information.

More information about ArcGIS Service configuration files can be found here:
http://webhelp.esri.com/arcgisserver/9.3/java/index.htm#cfg_file_service.htm

Augmented Reality in Commercial Applications

ANZ Bank have recently released an iPhone app for real estate analysis

I have to say, the application looks amazing, one of the best I've used. The user can search for houses on the market by location, price range, etc. Then the results can be displayed in three different ways:
1) A list (boring)
2) As points on a map (good)
3) In an augmented reality view (great)

The map view is cool, but is becoming expected within these kinds of applications. The thing really impressed me was the AR view. Augmented reality uses phone location (GPS + compass + accelerometer) to add information to the camera display. In this implementation, ANZ have presented the real-estate information to the user in a format that they can relate to (more so than a map). It is a bit of a gimmick, mainly because it heavily relies on the compass sensor at this stage, but impressive none the less.

What would be really cool, would be to see some live GPS data streaming to the display. For example, someone is waiting for their bus, so they open the AR display and can see that they have missed it by a few minutes. Or to interact with nearby friends, Twitter would become more like a phone call and less like a messaging service.

It's great stuff and am looking forward to seeing other implementations. I've already used the Layar app to overlay Wikipedia articles and other sources on the camera display. Let me know if you know of any other good ones!

Flex Sandbox Security Error #2048

I came across an annoying limitation of Flex today. After deploying my ArcGIS Flex application to an IIS server, I was having trouble accessing a Web Map Service (WMS). I could run the application locally and access the remote services, but deploying the application resulted in a Sandbox Security Error #2048:

"SecurityError: Error #2048: Security sandbox violation: http://mymachine/myflexapp/index.swf cannot load data from http://myserver.mydomain.com/arcgis/rest/services/myMapService/1?f=json..."

After perusing various AtionScript boards, I was finally directed to the solution on the ArcGIS Server 9.3 online help:

http://resources.esri.com/help/9.3/arcgisserver/apis/flex/help/index.html

"To access data from a different server than the one hosting your Flex application, the remote server needs to have a cross-domain file in the root directory."

As I understand it, this helps to prevent cross-site request forgery. Fair enough, but it's still annoying that Flex has this requirement, while many other client apps do not. My Flex application accesses the WMS service just like any other. The difference is this:

Web applications operate in a common environment, your browser. This environment shares common resources, such as cookies, which is a security concern.

Short of publishing the application as a standalone AIR app, I must rely on configuring the WMS server to allow access. It seems that this goes against some of the objectives of the OGC (seamless interoperability). I suppose javascript applications manage this is a similar way...?

Initially I was concerned about opening security holes, so I checked some well-known websites for examples on how to configure the file:
http://www.youtube.com/crossdomain.xml
http://www.google.com/crossdomain.xml

I was able to resolve the problem by asking the WMS Server administrators to set up the crossdomain.xml file at the root level of their website as described on the ESRI site. An example crossdomain.xml might be:

<?xml version="1.0" ?>
<cross-domain-policy>
<allow-access-from domain="*" />
<site-control permitted-cross-domain-policies="all" />
<allow-http-request-headers-from domain="*" headers="*" />
</cross-domain-policy>

This may not always be possible so an alternative approach is to configure a reflector on my own web server. Basically this will forward requests to the real destination, but to the web app, they appear to be coming from my own server. Obviously, copyright issues need to be considered when setting up something like this. The benefit of having a reflector is that you only need to configure a single crossdomain.xml file on the web server.

Amateur Radio and Python GDAL

I have installed GDAL python bindings to my Karmic Ubuntu box and my MacBook. I intend to test this Python library on the Landsat imagery which I have downloaded (Canberra and surrounds of course). There were a number of sample gdal-python scripts included with the MacOSX download. I still need to check these out.

After a number of Amateur Radio websites have caught my attention, I came across the linux software "GPredict". This software provides a mapping interface which displays the locations of hundreds of satellites. Also, it includes a window to configure a radio device. I shall investigate further.

Processing ESRI Shapefiles with Python

For the last few days, I've been working on translating data between an ESRI ArcSDE instance and KML for display within Google Earth. For 'point' and 'line' feature classes this was a fairly easy task, but translating 'polygon' feature classes turned out to be more hassle than I was expecting, so I thought I'd share my experiences. I was working with an ArcSDE 9.2 instance, but I am quite confident that very little is different, in terms of the geometry schema, in ArcGIS 9.3 software.

The ArcGIS geometry storage schema is quite complex, so I will try to restrict this post to the fundamental concepts. Basically, accessing geometry proceeds as follows:

Each feature class is made up of many rows.
Each row has an associated geometry.
A geometry may consist of several parts.
Each part is composed of either a single point or an array of points
Each point has X, Y and Z values.

Additionally, the geometry type can be identified by the 'Type' attribute of the geometry class.

The full geoprocessor model for ArcGIS 9.3 can be found here

Now, single point feature classes are easily processed as follows:

rows = gp.SearchCursor(myFC)
thisRow = rows.next()

while (thisRow):
    rowGeom=thisRow.getvalue("SHAPE")
    for part in range(rowGeom.partcount):
        geomPart=rowGeom.getPart(part-1)
        thisPoint=geomPart.get(0)
    thisRow=rows.next()

Extracting single lines is a bit trickier, as you must process the array of points:

rows = gp.SearchCursor(myFC)
thisRow = rows.next()

while (thisRow):
    rowGeom=thisRow.getvalue("SHAPE")
    for part in range(rowGeom.partcount):
        geomPart=rowGeom.getPart(part-1)
        pointArray=geomPart.get(0)
        thisPoint=pointArray.next()
        while (thisPoint):
            #Do some processing on thisPoint.
            thisPoint=pointArray.next()
    thisRow=rows.next()

Extracting polygons is where things get complicated, due to the way in which 'inner rings' are stored. Inner rings are linear rings that contain an area which is not to be included in the original polygon. Visually, these appear as 'holes' inside the polygon. An example of this is the polygon representing the political region of South Africa:



The empty area in the middle of South Africa is Lesotho. In this example the external boundary of South Africa forms the 'outer ring' and the border of Lesotho forms an 'inner ring'. The region bounded by the outer ring minus the region bounded by the inner ring(s) defines the polygon.

The ArcGIS schema stores polygon information in the geometry's array of points. As described above, the point array is made up of a list of points which can accessed using the iterator function 'next()'. For polygons, the list stores an initial 'outer ring' and optionally, several 'inner rings'. The outer ring is always the first contiguous list of points in the array. An inner ring can then be identified by a null point followed by the list of points forming a second linear ring. Additional inner rings may be added in the same way. This data structure dictates that processing polygons with inner rings can be done using the following logic:

rows = gp.SearchCursor(myFC)
thisRow = rows.next()

while (thisRow):
    rowGeom=thisRow.getvalue("SHAPE")
    for part in range(rowGeom.partcount):
        geomPart=rowGeom.getPart(part-1)
        pointArray=geomPart.get(0)
        thisPoint=pointArray.next()
        while (thisPoint):
            #Do some processing on thisPoint.
            thisPoint=pointArray.next()

            #Checking if this polygon contains an inner ring.
            if not thisPoint:
                thisPoint=pointArray.next()
                if thisPoint:
                    innerRing=True
                    #Begin point list for inner ring
    thisRow=rows.next()

 

I will extend this logic code to create KML in my next post.

Cheers,
JamesO