POST
|
I use GIS everyday, I have an ArcMap licence 10.1 on my PC, but the PC with Spatial Analyst and 3D licence has remained at 9.3 due to cuts of the budget of the office where I work. We use Spatial Analysis not on regular basis, that's why. Since I should run the script with several maps, and our DTM is very heavy, I would prefer to have the script adapted to 9.3, but I understand that might be very diffcult bye Ialina
... View more
10-20-2015
05:54 AM
|
0
|
1
|
1234
|
POST
|
Sorry to have answered so late. I confirm that the gdb was created with ArcGis 9.3. I had already seen that the da commands were not available in 9.3 and was hoping that someone could give me the clue how to bypass the problem, as I tried to bypass the arcpy function, also not available, with something I found in the net. But this is too much for my little knowledge and I have to give up. 3,600 € for the software updating seems to be too much for the needs of the work and for the budget of the company. Thank you anyway Bye Ialina
... View more
10-20-2015
01:42 AM
|
0
|
3
|
1234
|
POST
|
Hi, unfortunately I have the spatial analyst linked to an ArcGis 9.3 licence, so I I tried the script adapted to python 2.6 and ArcGIS9.3, so far I've been able to do searching in the net (I'm a beginner). I run the script and I've got his message: Traceback (most recent call last): File "C:\Python26\Lib\site-packages\pythonwin\pywin\framework\scriptutils.py", line 326, in RunScript exec codeObject in __main__.__dict__ File "C:\Documents and Settings\lfantinato\Documenti\zone_per2.py", line 34, in <module> gp.env.Workspace = r"G:\ialina2\dtm\zon_st.gdb" File "C:\Python26\lib\site-packages\win32com\client\dynamic.py", line 522, in __getattr__ raise AttributeError("%s.%s" % (self._username_, attr)) AttributeError: esriGeoprocessing.GpDispatch.1.env what are the problems? is it possible to fix them? If someone is able to help me, I would be very grateful, thank you Ialina this is the script adjusted by me: def GetPixelValueForPercentile(dctper, percentile): ##"""will return last pixel value where percentile LE searched percentile.""" for k in sorted(dctper.keys()): perc = dctper if perc <= percentile: pix_val = k else: break return pix_val def FieldExist(featureclass, fieldname): ## """Check if field exists""" import win32com.client gp = win32com.client.Dispatch("esriGeoprocessing.GpDispatch.1") fieldList = gp.ListFields(featureclass, fieldname) return len(fieldList) == 1 import win32com.client gp = win32com.client.Dispatch("esriGeoprocessing.GpDispatch.1") ## settings ras = r"G:\ialina2\dtm\dtmx10clipgrd" fc = r"G:\ialina2\dtm\polig_250clip" ## ras = r"C:\Path\To\Folder\Or\Geodatabase\With\Raster" ## fc = r"C:\Path\To\Folder\Or\Geodatabase\With\PolygonFeatureclass" lst_perc = [2, 5, 10, 90, 95, 98] ## list of percentiles to be calculated fld_prefix = "Perc_" gp.env.Workspace = r"G:\ialina2\dtm\zon_st.gdb" ##arcpy.env.workspace = "IN_MEMORY" gp.env.overwriteOutput = True ## add fields if they do not exist in fc ## be aware, field that exist will be overwritten! print "filling field list and adding fields" flds = [] for perc in lst_perc: fld_perc = "{0}{1}".format(fld_prefix, perc) flds.append(fld_perc) if not FieldExist(fc, fld_perc): gp.AddField_management(fc, fld_perc, "LONG") print "flds={0}".format(flds) ## Enable Spatial analyst gp.CheckOutExtension("Spatial") ## loop through polygons print "loop through polygons" flds.append("SHAPE@") i = 0 with gp.da.UpdateCursor(fc, flds) as curs: for row in curs: i += 1 print "Processing polygon: {0}".format(i) polygon = row[flds.index("SHAPE@")] lst_parts = [] if polygon.partCount == 1: for part in polygon: for pnt in part: x, y = pnt.X, pnt.Y lst_parts.append(gp.Point(x, y)) else: for part in polygon: lst_crds = [] for pnt in part: x, y = pnt.X, pnt.Y lst_crds.append(gp.Point(x, y)) lst_parts.append(lst_crds) print " - lst_parts={0}".format(lst_parts) ## Execute ExtractByPolygon (you can't send the polygon object) print " - ExtractByPolygon..." ras_pol = gp.sa.ExtractByPolygon(ras, lst_parts, "INSIDE") outname = "ras{0}".format(i) ras_pol.save(outname) print " - saved raster as {0}".format(outname) ## create dictionary with value vs count print " - fill dict with Value x Count" flds_ras = ("VALUE", "COUNT") dct = dict((row[0],row[1]) for row in gp.da.SearchCursor(outname, flds_ras)) ## calculate number of pixels in raster print " - determine sum" cnt_sum = sum(dct.values()) print " - sum={0}".format(cnt_sum) ## loop through dictionary and create new dictionary with val vs percentile print " - create percentile dict" dct_per = {} cnt_i = 0 for val in sorted(dct.keys()): cnt_i += dct[val] dct_per[val] = cnt_i / cnt_sum ## loop through list of percentiles print " - iterate percentiles" for perc in lst_perc: ## use dct_per to determine percentiles perc_dec = float(perc) / 100 print " - Perc_dec for is {0}".format(perc_dec) pixval = GetPixelValueForPercentile(dct_per, perc_dec) print " - Perc for {0}% is {1}".format(perc, pixval) ## write pixel value to percentile field print " - Store value" fld_perc = "{0}{1}".format(fld_prefix, perc) row[flds.index(fld_perc)] = pixval ## update row print " - update row" curs.updateRow(row) ## return SA license gp.CheckInExtension("Spatial")
... View more
10-01-2015
06:01 AM
|
0
|
6
|
856
|
Online Status |
Offline
|
Date Last Visited |
11-11-2020
02:24 AM
|