Python Final Project Assistance XY points from table to final IDW raster

[annesanta_maria2]

Hi,
I am taking an independent study python class and the professor is just learning python as well. For my final project I want to do the following:

1. add xy data from several tables-Make XY events layer-using WGS84 spatial reference
2. convert that data to a permanent shapefile- copy features
3. take those points and list fields-arcpy.listfields
4. run IDW across several fields for all the points

I can run that process for one set of data from table to IDW raster but I cannot seem to figure out how to turn that into a code that automates twenty tables and many fields in those points layers.

Any suggestions on how to start? Essentially I need a script that I can turn into a toolbox for a final project in April. I thought I had this under control until I realized how hard it is for me to grasp python. I am just learning and I DON'T WANT SOMEONE JUST TO TELL ME HOW to do it. I want to really understand.

import arcpy
from arcpy import env

# set workspace
arcpy.env. workspace ="whatever the workspace would be C:/user/etc"
tablelist = arcpy.Listtables()
.....# this is where I get stuck how do take this table list and turn into into several events layers. Do I need to list all the fields in the table. 

for tables table list:
     arcpy.MakeXYEvents layers(enter parameters)


arcpy.CopyFeatures

I guess what I am looking for is the best next step to take. I appreciate any assistance.
Anne

[XanderBakker]

Hi Anne,

Good for you that want to understand the Python code that solves your problem. In order to give you some pointers, it is necessary to know a little more:

• what version of ArcGIS are you using?


• do you have any programming knowledge?

About Python
In case you don't have (much) experience with Python, consider the next resources:

The thread below provides some useful links to resources:
http://forums.arcgis.com/threads/96424-Learning-Python

A good way to start can be using the Python window in ArcGIS, since you obtain context sensitive Help. Also after you manually execute a ArcGIS tool, you can use the option Copy as Python Snippet in the Results window.

When working with looping through feature classes and tables, have a look at cursors (especially the data access cursors): Accessing data using cursors

And you can follow free courses online at Esri:
http://training.esri.com/gateway/index.cfm?fa=search.results&searchterm=python

Depending on your knowledge and needs, there are many options. Also the Forums are a very useful resource. Always look at the posts by jscheirer (Jason Scheirer). He is on the Python team at Esri and also posts on: http://arcpy.wordpress.com/

I heard about the book "Programming ArcGIS 10.1 with Python" which appears to be a good resource, and there will be a new one "Mastering Python Programming for ArcGIS".


About your question
I assume your post in this threads is related to what you are asking here:
http://forums.arcgis.com/threads/69130-ArcGIS10.0-PythonWindow-Make-XY-Event-Layer-.csv-amp-WGS-1984...

From this thread I extracted that in your workspace you have a number of DBF tables:
AB10_.dbf, AB11_.dbf, AB15_.dbf, CFHMS_.dbf, CFMHN_.dbf, FA_.dbf, HSFHome_.dbf, HSFVAN_.dbf, KF18_.dbf, KF1W_.dbf, LASA_.dbf, MB_HAck_.dbf, PCM10_.dbf, PMB28_.dbf, pmc9_.dbf, SHROCK_.dbf, simo09_.dbf, SIMO94_.dbf, SP150_.dbf, sp42_.dbf, T13East_.dbf, T13West_.dbf, TB53_.dbf, TB63_.dbf, WSFBOW_.dbf, WSFDALLAS_.dbf
I assume they all have the the same coordinate fields: "Lat" and "long". Please note that in your code in the other thread you assign the Lat to x and Long to y (it should be the other way around!).

After converting the DBF's to shapefiles, you want to perform inverse distance weighted interpolation on each attribute. What I don't completely understand is whether you want to merge all the tables together first, or whether for each table/shapefile and attribute the IDW should be performed. Main question here is what you want to do with the rasters that result from the operation.

Looking at the table names, it seems that each table covers a different area. Is that right? In that case merging the features would be an option, but if you want to do that, you will have to look at the schema's of the tables (do they have the same fields?). If the schema's are different, you will probably have to use cursors and more programming to do that or some nasty field mapping. I could advise you on that, but will need to have a look at the data first. I don't know if that's an option for you...

Please note that looping through each field of a shapefile to perform the IDW will result in more rasters then you want, since your lat and long fields will also be used. You will have to include some more intelligence to do this right.

Let me know if I can help you.

Kind regards,

Xander

[annesanta_maria2]

Hi,
Thank you so much for your assistance 🙂 I am really working on this! I am a grad student and it seems all the jobs post school want python which of course the geography department doesn't teach. So the independent study class is the best option. Right now I am working through Learning Python the Hard Way and the ArcGIS python book. I've attached an example of one of the data sheets I am working with. All of the excel/.dbf files look the same.

The end result of the project will be a custom tool so that you can loop through the tables and create event layers and then copy them to a folder and have feature classes created. I think what I will do is leave the looping out of the IDW, just because of complications and because it is probably better if the people I am working with specify what data they want to do IDW on. If that makes sense. So the end result I am working toward is a tool where "they-the users" specify the location of their table, the tool then loops through the tables and creates points layers. Then they specify which table they want to use and the field they want to use for IDW.

So mainly I am working through looping with the tables and XY events layers. I have fixed the code so that it does work better than it did last time. Here is an example:

>>>import arcpy
>>> from arcpy import env
>>> env.workspace = "F:/work/python" # once I turn this into a tool it will be changed to GetParameters as text()
>>> tables = arcpy.ListTables()
>>> xc = "long"
>>> yc = "lat"
>>> out_layer = "F:/work/python" # once I turn this into a tool it will be changed to GetParameters as text()
>>>spref = 4326 #I did this as the hardcode for WGS 84 but there may be a better way. I am still looking.
>>> for tbl in tables:
...     arcpy.MakeXYEventLayer_management(tbl, xc, yc, out_layer, spref)

New error which I checked out. Runtime error  Traceback (most recent call last):   File "<string>", line 2, in <module>   File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\management.py", line 6348, in MakeXYEventLayer     raise e ExecuteError: ERROR 000212: Cannot create XY event source Failed to execute (MakeXYEventLayer)
AND I think it has to do with the coordinate system reference that I am using.

The next step is to do arcpy.CopyFeatures_management(outlayer, savedlayer)

I really appreciate the help and I think I am getting closer to figuring this out and learning python everyday.

[XanderBakker]

Hi Anne,

Indeed you spatial reference is not defined correctly. Change:

spref = 4326

... into:

spref = arcpy.SpatialReference(4326)

This creates an actual spatial reference object from the WKID.

In case the tables contain the same attributes and represent different contiguous areas, it will be better to merge the points together before you perform the interpolation. Interpolation is, as the name suggests, not an extrapolation. This means that the values at the borders of the raster resulting from the interpolation will not be very reliable.

Another thing I notice in your code is the line:

out_layer = "F:/work/python" 

This line of code is pointing to a folder, which in this case is the same folder as the folder that holds the tables. Instead of a folder you should provide a name here, like for instance "measurements". The layer created by this tool is temporary and resides in memory. Within the loop (for tbl in tables) the arcpy.CopyFeatures_management should be executed. This will require a location and name of the output featureclass to be created.

Let's assume you would have an input folder with tables and an output folder to store the featureclass. Your code could like something like this:

import arcpy, os
from arcpy import env

ws_in = "F:/work/python" # could use arcpy.GetParameterAsText(0)
ws_out = "F:/work/python/output" # could use arcpy.GetParameterAsText(1)

env.workspace = ws_in
tables = arcpy.ListTables()

xc = "long"
yc = "lat"

spref = arcpy.SpatialReference(4326)

XYeventname = "measurements"

for tbl in tables:
    arcpy.MakeXYEventLayer_management(tbl, xc, yc, XYeventname, spref)

    # assuming you have DBF tables, strip the extension:
    fclass = os.path.join(ws_out, tbl.strip(".dbf"))

    # make the featureclass
    arcpy.CopyFeatures_management(XYeventname, fclass)

In the code for each table an XY event layer is created in memory. It will be named "measurements" and this layer is overwritten with each table. Using the os module the output workspace "ws_out" is joined with the table name, although in this case the ".dbf" is stripped from the name. The resulting variable "fclass" contains the path and name to the featureclass to be created. Using the copy features tool, the layer is written to a featureclass.

What is important here, is that the variable "fclass" defined what type of output is written. If it references to a name in a folder , it will be a shapefile, but if you are pointing to a name in a file geodatabase, it will create a featureclass in a geodatabase. The reason for not using the same input and output folder, is that a shapefile is created with the same name as the (DBF) table. However, a shapefile consists of a DBF file as well to store the attributes. If you would write the shapefile with the same table name to the same input folder you would get some nasty errors.

Another thing is the attached Excel file. If you would work on this it would change things a bit, since when working with Excel files, the Excel file itself is the workspace and worksheets or named data ranges will be the tables.

Kind regards,

Xander

[annesanta_maria2]

Hi,
Sorry it took so long to say thank you. After working through the tables and the code I finally got it to run properly without any errors. It was amazing. Thank you. I put together this code for the IDW so that someone can select several different fields from the same table and run IDW on the data for that set of points. It's really crude but it works I suppose. Again I want to thank you for all your help. If you have any suggestions on the stuff at the bottom let me know.


import arcpy
from arcpy import env
from arcpy.sa import *

arcpy.CheckOutExtension("Spatial")

env.workspace = "C:/Users/Anne/Desktop/Python/finalproj"

#Set Variables
points = "Export_Output.shp"
zfield1 = "pH"
zfield2 = "Pphate"
zfield3 = "K"
zfield4 = "CA"
zfield5 = "Mg"

#First Run
outidw1 = Idw(points, zfield1)
outidw1.save("C:/Users/Anne/Desktop/Python/finalproj/data/outputs/ab10ph")

#Second Run
outidw2 = Idw(points, zfield2)
outidw2.save("C:/Users/Anne/Desktop/Python/finalproj/data/outputs/ab10pphate")

#Third Run
outidw3 = Idw(points, zfield3)
outidw3.save("C:/Users/Anne/Desktop/Python/finalproj/data/outputs/ab10ll")

#Fourth Run
outidw4 = Idw(points, zfield4)
outidw4.save("C:/Users/Anne/Desktop/Python/finalproj/data/outputs/ab10ca")

#Fifth Run
outidw5 = Idw(points, zfield5)
outidw5.save("C:/Users/Anne/Desktop/Python/finalproj/data/outputs/ab10mg")

[annesanta_maria2]

I'm having a lot of trouble getting the parameters as text. For some resone the arcpy.env.workspace is not defined everytime I run.

import arcpy, os
from arcpy import env
ws_in = arcpy.GetParameter(0)
ws_out = arcpy.GetParameter(1)

arcpy.env.workspace = ws_in
arcpy.env.overwriteOutput = True

tables = arcpy.ListTables()

# Sets local variables
xc = "long"
yc = "lat"

spref = arcpy.SpatialReference(4326)

XYeventname = "measurements"

#Loops through the table and makes points which are saved to the output location
for tbl in tables:
    arcpy.MakeXYEventLayer_management(tbl, xc, yc, XYeventname, spref)
    fclass = os.path.join(ws_out, tbl.strip(".dbf"))
    arcpy.CopyFeatures_management(XYeventname, fclass)

[annesanta_maria2]

Sorry. I am having problems getting the parameters when I create the tool.

So I added an input and output file for the tool. The input is set up to read workspaces and the output is set up as an output workspace. I get this error.

Traceback (most recent call last):
  File "C:\Users\Anne\Desktop\Python\finalproj\tabletopoint.py", line 29, in <module>
    fclass = os.path.join(ws_out, tbl.strip(".dbf"))
  File "C:\Python27\ArcGIS10.1\Lib\ntpath.py", line 96, in join
    assert len(path) > 0
TypeError: object of type 'geoprocessing value object' has no len()

[annesanta_maria2]

So ignore everything. I figured it out!
🙂

[XanderBakker]

Hi,
Sorry it took so long to say thank you. After working through the tables and the code I finally got it to run properly without any errors. It was amazing. Thank you. I put together this code for the IDW so that someone can select several different fields from the same table and run IDW on the data for that set of points. It's really crude but it works I suppose. Again I want to thank you for all your help. If you have any suggestions on the stuff at the bottom let me know.


import arcpy
from arcpy import env
from arcpy.sa import *

arcpy.CheckOutExtension("Spatial")

env.workspace = "C:/Users/Anne/Desktop/Python/finalproj"

#Set Variables
points = "Export_Output.shp"
zfield1 = "pH"
zfield2 = "Pphate"
zfield3 = "K"
zfield4 = "CA"
zfield5 = "Mg"

#First Run
outidw1 = Idw(points, zfield1)
outidw1.save("C:/Users/Anne/Desktop/Python/finalproj/data/outputs/ab10ph")

#Second Run
outidw2 = Idw(points, zfield2)
outidw2.save("C:/Users/Anne/Desktop/Python/finalproj/data/outputs/ab10pphate")

#Third Run
outidw3 = Idw(points, zfield3)
outidw3.save("C:/Users/Anne/Desktop/Python/finalproj/data/outputs/ab10ll")

#Fourth Run
outidw4 = Idw(points, zfield4)
outidw4.save("C:/Users/Anne/Desktop/Python/finalproj/data/outputs/ab10ca")

#Fifth Run
outidw5 = Idw(points, zfield5)
outidw5.save("C:/Users/Anne/Desktop/Python/finalproj/data/outputs/ab10mg")

Hi Anne,

Glad to hear you worked it out. I would like to tell you something on using dictionaries and looping through data. In your sample where you calculate the IDW for 5 fields, you could have used a dictionary containing the input field names and output raster names.

import arcpy, os

arcpy.CheckOutExtension("Spatial")
ws_out = "C:/Users/Anne/Desktop/Python/finalproj/data/outputs"
env.workspace = ws_out # set it to your output workspace
arcpy.env.overwriteOutput = True

#Set Variables
points = "C:/Users/Anne/Desktop/Python/finalproj/Export_Output.shp"

# dictionary with fieldnames and output raster names
d = {'pH': 'ab10ph', 'Pphate': 'ab10pphate', 'K': 'ab10ll',
     'CA': 'ab10ca', 'Mg': 'ab10mg'}

# loop through the dictionary
for fieldname, rasname in d.items():
    # are you sure you want to accept default cellsize, power and search radius?
    outidw = arcpy.sa.Idw(points, fieldname)
    outidw.save(rasname) # no path included, so saved to (output) workspace

# return the sa license
arcpy.CheckInExtension("Spatial")

In case you want to use the field name in the output raster name you could a list as well:

l = ['pH', 'Pphate', 'K', 'CA', 'Mg']
ras_prefix = 'ab10'
for name in l:
    outidw = arcpy.sa.Idw(points, name)
    outidw.save('{0}{1}'.format(ras_prefix, name))

Another thing to remember is that you are not specifying any cellsize, power and search radius settings. If you have your environment settings setup correctly and your points are also located outside your study area, this might give appropriate results. If your area is cutoff (since it will by default take the extent of the points) you may want to define things like arcpy.env.snapRaster and arcpy.env.extent.

Kind regards,

Xander

[annesanta_maria2]

Hmm. Would this be able to be turned into a tool so that the professor could put in any set of points and then specify the fields to be used? (I definitely need to take another class in Python. The ArcGIS book is not entirely helpful, but it isn't bad either). I'm not sure about extents right now and cell size, but I am looking into the data some more to see what the area covered is and how much interpolation they are doing outside points etc. So I am waiting on that, but this is super helpful. I really really appreciate it.

import arcpy, os

arcpy.CheckOutExtension("Spatial")
ws_out = #GetParameterAsText(0)
env.workspace = # GetParameterAsText(1)
arcpy.env.overwriteOutput = True

#Set Variables
points = GetParameterAsText(2)

#Is there a way to set this up so the fields could be selected from a dropdown list in the actual toolbox. Then those would be used to create a list.....
l = ['pH', 'Pphate', 'K', 'CA', 'Mg']
ras_prefix = 'ab10'
for name in l:
    outidw = arcpy.sa.Idw(points, name)
    outidw.save('{0}{1}'.format(ras_prefix, name))


# return the sa license
arcpy.CheckInExtension("Spatial")