Help with creating geostat layer in python

3676
4
Jump to solution
04-24-2014 08:21 AM
BenSciance
New Contributor II
I'm trying to loop through feature classes and create a geostatistical layer and then convert the geo layer to points.

I've attempted to debug this problem myself, but would appreciate extra eyes on my code in case someone else can spot an error.

Here is a resource that I've referenced trying to fix my code:
http://forums.arcgis.com/threads/92508-Need-Help-Automation-of-Kriging-using-Model-builder-or-python...

import arcpy from arcpy import env from arcpy.sa import* from arcpy.ga import*   arcpy.env.workspace = "C:\Users\mbs7038\Documents\New_Landsat_Imagery\For_Area_Calc\Gauges_For_AccuracyAssessment\S_Gauges\S_Gauges.gdb" arcpy.env.overwriteOutput=True  lyrpath = "C:\Users\mbs7038\Documents\New_Landsat_Imagery\For_Area_Calc\Gauges_For_AccuracyAssessment\S_Gauges\GA_lyrs"  ptpath="C:\Users\mbs7038\Documents\New_Landsat_Imagery\For_Area_Calc\Gauges_For_AccuracyAssessment\S_Gauges\GA_Krig.gdb"   # Check out the ArcGIS Geostatistical Analyst extension license arcpy.CheckOutExtension("GeoStats")   #modelTemplate="C:\Users\mbs7038\Documents\New_Landsat_Imagery\For_Area_Calc\Gauges_For_AccuracyAssessment\S_Gauges\GA_Model.lyr"  modelxml=("C:\Users\mbs7038\Documents\New_Landsat_Imagery\For_Area_Calc\Gauges_For_AccuracyAssessment\S_Gauges\EmpiricalBayesianKriging.xml")   print 'listing feature classes' fcs=arcpy.ListFeatureClasses("a734184")  for fc in fcs:     print fc     try:         print 'listing fields'         fields = arcpy.ListFields(fc,"Wlevel")         for field in fields:             print ("{0} is a type of {1} with a length of {2}".format(field.name, field.type, field.length))          zfield="Wlevel"          outlyr=lyrpath+"/"+fc+"_GA.lyr"         outpts=ptpath+"/"+fc+"_GA_pt"         #cellsize=0.0008333333333          print'creating geostatistical layer'         infc=fc + " Wlevel"         arcpy.GACreateGeostatisticalLayer_ga(modelxml,infc,outlyr)          print 'converting GA layer to grid'         arcpy.GALayerToPoints_ga(outlyr,fc,zfield,outpts)      except:         print arcpy.GetMessages()  print 'script complete'


I receive these error messages: a734184_GA.lyr does not exist or is not supported 

I've checked my output folder that I've defined to store the new geostat layers, and there is nothing there, even though I received no error messages when running the CreateGeostatisticalLayer tool......

>>>  listing feature classes a734184 listing fields Wlevel is a type of Double with a length of 8 creating geostatistical layer converting GA layer to grid Executing: GALayerToPoints C:\Users\mbs7038\Documents\New_Landsat_Imagery\For_Area_Calc\Gauges_For_AccuracyAssessment\S_Gauges\GA_lyrs/a734184_GA.lyr C:\Users\mbs7038\Documents\New_Landsat_Imagery\For_Area_Calc\Gauges_For_AccuracyAssessment\S_Gauges\S_Gauges.gdb\a734184 Wlevel C:\Users\mbs7038\Documents\New_Landsat_Imagery\For_Area_Calc\Gauges_For_AccuracyAssessment\S_Gauges\GA_Krig.gdb\a734184_GA_pt ALL Start Time: Thu Apr 24 12:12:05 2014 Failed to execute. Parameters are not valid. ERROR 000732: Input geostatistical layer: Dataset C:\Users\mbs7038\Documents\New_Landsat_Imagery\For_Area_Calc\Gauges_For_AccuracyAssessment\S_Gauges\GA_lyrs/a734184_GA.lyr does not exist or is not supported Failed to execute (GALayerToPoints). Failed at Thu Apr 24 12:12:05 2014 (Elapsed Time: 0.00 seconds) script complete >>> 
0 Kudos
1 Solution

Accepted Solutions
SteveLynch
Esri Regular Contributor
Ben

The geostatistical layer that you create is an in memory layer. Therefore
outlyr=lyrpath+"/"+fc+"_GA.lyr" should not have a path and an extension. Just a name, viz., 'mylyr'

You could, for debugging, persist this in memory layer using arcpy.SaveToLayerFile_management('mylyr', r'c:\temp\galyr.lyr')
and then drag this .lyr into ArcMap.

-Steve

View solution in original post

0 Kudos
4 Replies
SteveLynch
Esri Regular Contributor
Ben

The geostatistical layer that you create is an in memory layer. Therefore
outlyr=lyrpath+"/"+fc+"_GA.lyr" should not have a path and an extension. Just a name, viz., 'mylyr'

You could, for debugging, persist this in memory layer using arcpy.SaveToLayerFile_management('mylyr', r'c:\temp\galyr.lyr')
and then drag this .lyr into ArcMap.

-Steve
0 Kudos
BenSciance
New Contributor II
Interesting....Thanks for the info.  Must have missed that when I was reading up on all things geostatistical analyst.

I'll give it a go.
0 Kudos
SteveLynch
Esri Regular Contributor
Ben

Some additional comments regarding your script.

- EBK is a geoprocessing tool and therefore you don't need to 'create' it via the CreateGeostatisticalLayer tool.
- if you are creating raster output you can use the environmental MASK which will speed up execution a lot
- it appears that your data are in GCS (geographic coordinate system), not a good idea, rather project the data.


-Steve
0 Kudos
BenSciance
New Contributor II
Thanks for the info Steve.  I'll definitely keep this is mind.
Ben
0 Kudos