Latest Contributions by XanderBakker

hello Xander. thanks for you reply. The only data I have is the medical records, so I know the region the patient came from. I have a map of the regions and I want to put all the records there, in oder to have them classified per region, or cancer type etc. I don't have any summarized records, this is what I want to create PS: sorry about the cancer talk, but this is the subject If the tabular data is already in ArcGIS table format you can read in the help topic "Summarizing data in a table" how you can summarize it. If your data is in Excel you can use Pivot Table to summarize the data. Or use the options of you DBMS of choice. You summarize on the field with region "key". This key is also used for joining the summarized data to the regions shapes. Hope this helps, Xander ... View more

You can create a Fishnet. See Help topic Create Fishnet (Data Management) Kind regards, Xander ... View more

Greetings to everyone. My problem is this: I want to analyze data which refer to cancer occurrences in my country and examine their spatial correlations. This means that my database consist of patient data, that have a reference to the region the patient comes from. My question is how can I relate all the data to the regions shapefile. That is, I have a polygon shapefile with all 52 regions, but, as expected, each region has more than one occurrence. One solution is to have 52 separate shapefiles, one for each region, and then produce spatial statistics for each. What I want to do, is have them unified in one database, in order to make comparisons, for example, between region A & B etc. Any help is much appreciated. What you are looking for is normally done with a join. It's a good start to read this help topic: About joining and relating tables If you don't have any spatial component in your database with cancer occurrences (coordinates, addresses, medical centers) then the best thing to do is summarize the data in the database. I suppose you will want to know the number of cases. The resulting table can be joined to the shape. Further analyses like occurrences in relation to population can be performed in the joined featureclass. If there is a spatial component in your database, make it spatial (geocode addresses or make XY Event layer). Next, spatial join the points to the featureclass with the regions. Since it is a 1:many relationship you can summarize the attributes. Kind regards, Xander ... View more

In addition the points pointed out correctly by Jake, there are a few other things: 1) The "env.workspace" is pointing to a featureclass (shapefile) instead of a workspace. Better change this to env.workspace = "K:/Projects/Project Number" 2) it is OK to check for the existence of the featureclass "shapefile.shp", but it won't work when the workspace is not correctly set change: verification = arcpy.Exists(fc) print verification del verification to: if arcpy.Exists(fc):     # continue here with indented code 3) You import "fileinput, os, string", but since you don't use them, it may be better not to import those. 4) It is not necessary to delete the cursor since you are using a "with" statement. 5) it may be wise to check if the column exists before adding the field to avoid errors outFld = "Field3" if len(arcpy.ListFields(fc, outFld)) == 0:     arcpy.Addfield_management (fc, outFld, "DOUBLE") Kind regards, Xander ... View more

Hi Fred, You're right that .TXT is not a supported output format in the Table to Table conversion tool. If you work with a featureclass instead of a table the "Export Feature Attribute to ASCII (Spatial Statistics)" would be an option. If not than you can use the snippet below I found on stack exchange: http://gis.stackexchange.com/questions/17933/export-table-to-x-y-z-ascii-file-via-arcpy import arcpy,csv table =r'c:\path\to\table' outfile = r'c:\path\to\output\ascii\text\file'      #--first lets make a list of all of the fields in the table fields = arcpy.ListFields(table) field_names = [field.name for field in fields] with open(outfile,'wb') as f:     w = csv.writer(f)     #--write all field names to the output file     w.writerow(field_names)     #--now we make the search cursor that will iterate through the rows of the table     for row in arcpy.SearchCursor(table):         field_vals = [row.getValue(field.name) for field in fields]         w.writerow(field_vals)     del row Please note that things may go wrong with binary (blob) fields and perhaps with Null values as well. Kind regards, Xander ... View more

Hi Teresa, Sounds OK. Could you share a small part of your data to see if it can be reproduced? Kind regards, Xander ... View more

The code returns me an error Runtime error <type 'exceptions.IOError'>: UNKNOWN layerName=r'M:\UTILISATEURS\Pierre\SERAIL\MISE_A_JOUR.gdb\TEST' I work from a class annotation. I don�??t want to change the data source,  but simply change the label expression or display rather like my data is  a class annotation. Hi Pierre, I notice your layername. I used the name of the layer as it appears in the table of contents (TOC) of the mxd. You are using a reference to the data source itself and that's probably not how the layer is named in the TOC. So change this: layerName=r'M:\UTILISATEURS\Pierre\SERAIL\MISE_A_JOUR.gdb\TEST' to this: layerName='TEST' ... if your layer is called 'TEST' in the TOC. You could also try with manually entering the label expression below in the layer properties: replace(replace([TEXTSTRING],"Ec.","Ecole"),"Egl.","Eglise") The error Runtime error <type 'exceptions.IOError'>: UNKNOWN you obtained refers to: Raised when an I/O operation (such as a print statement, the built-in open() function or a method of a file object) fails for an I/O-related reason, e.g., "file not found" or "disk full". Do you have write permissions to save the mxd? But the biggest difference is that you are not labeling based on a field, but you have an (feature-linked?) annotation class. This works completely different. In stead of a label expression there is the "Display Expression". This display expression is not (yet) supported by arcpy. In this case the best thing you can do is manually enter the properties of the annotation class and activate the Display TAB and press the Expression button: [ATTACH=CONFIG]28710[/ATTACH] In the Display Expression dialog enter the expression: replace(replace([TEXTSTRING],"Ec.","Ecole"),"Egl.","Eglise") [ATTACH=CONFIG]28711[/ATTACH] When there are more parts of text you want to change in the feature, it may be better to use a multiline expression like this: Function FindLabel ( [TextString] )   Label = replace([TextString], "Ec.", "Ecole")   Label = replace(Label, "Egl.", "Eglise")   Label = replace(Label, "Text to search for", "Text to replace it with")   'etc   FindLabel = Label End Function I tried this with a dummy dataset and if I hit the "Verify" button it shows the replaced text correctly. If I look at the map however it doesn't display the changed text. Please note that I switched on the option "Show MapTips using the display expression"... See steps below. [ATTACH=CONFIG]28714[/ATTACH] What is even stranger, when I hover over the text in the map it shows a kind of tool tip with the corrected text... This is odd.. maybe you're better of not using annotation... Kind regards, Xander ... View more

Hi Xander, Thanks for your response. Actually I was looking about "Feature Attachments", which involves a relationship table to store the attachments for features. Please go through "Query attachments" section in below link, which talks about performing the same task in .Net. http://help.arcgis.com/en/sdk/10.0/arcobjects_net/conceptualhelp/index.html#//0001000001qr000000 Python is more limited than ArcObjects (although you could implement COM types in Python and use ArcObjects). If you only want to save the attachments and not use any related data from the features then the code remains mostly the same: from arcpy import da import os # fc = r"C:\Project\_LearnPython\Attachments\test.gdb\yourFCname" tbl = r"C:\Project\_LearnPython\Attachments\test.gdb\yourFCname__ATTACH" fldBLOB = 'DATA' fldAttName = 'ATT_NAME' outFolder = r"C:\Project\_LearnPython\Attachments" with da.SearchCursor(tbl,[fldBLOB,fldAttName]) as cursor:    for row in cursor:       binaryRep = row[0]       fileName = row[1]       # save to disk       open(outFolder + os.sep + fileName, 'wb').write(binaryRep.tobytes())       del row       del binaryRep       del fileName In this case you access the attachment tabel ("yourFCname__ATTACH"). The name of the attachment is stored in the field "ATT_NAME" and the binary data in the field "DATA". If you want to access the feature itself, you will have to use the "REL_OBJECTID" field from the attachment table and query the featureclass on OBJECTID = REL_OBJECTID... Kind regards, Xander ... View more

It ran trough, but drop-value problem remains. (noticed a tiling effect, so there might be a projection problem here too -> looking into that currently) best regards, Henri Hi Henri, You're right the Blockstatistics causes this effect. Change it to FocalStatistics and that problem is most likely solved (parameters remain the same): meanElev = arcpy.sa.BlockStatistics(filledDEM, arcpy.sa.NbrRectangle(3, 3, "CELL"),"MINIMUM", "DATA") becomes: meanElev = arcpy.sa.FocalStatistics(filledDEM, arcpy.sa.NbrRectangle(3, 3, "CELL"),"MINIMUM", "DATA") Kind regards, Xander ... View more

Hello, I work on a table of about 1200 objects. I wonder if it is possible to write a standalone Python script to replace one word with another without changing the values �??�??of my table. As building label expressions. Exemple: I would like to replace Ec. for Ecole other exemple I would like to replace Egl. for Eglise Thanks for your help Hi Pierre, I guess the expression that you're looking for is something along the lines of: replace(replace([yourLabelFieldName],"Ec.","Ecole"),"Egl.","Eglise") I don't know why you want a standalone Python script to change the expression, but it could be something like: import arcpy mxdPath = r"C:\Project\_LearnPython\Ohio2.mxd" layerName = "OH30" fldName = "yourLabelFieldName" # in case of shapefile of personal geodatabase labelExpression = 'replace(replace([{0}],"Ec.","Ecole"),"Egl.","Eglise")'.format(fldName) # in case of file geodatabase # labelExpression = 'replace(replace("{0}","Ec.","Ecole"),"Egl.","Eglise")'.format(fldName) mxd = arcpy.mapping.MapDocument(mxdPath) for lyr in arcpy.mapping.ListLayers(mxd):     if lyr.name == layerName:         if lyr.supports("LABELCLASSES"):             lyr.showLabels == True             for lblClass in lyr.labelClasses:                 lblClass.expression = labelExpression mxd.save() del mxd Change the bold text to represent your situation. Activate the italic code in case your layer is in a file geodatabase. The code will: access your mxd loop through the layer until it finds the layerName provided change the label expression of the layer saves the mxd Kind regards, Xander ... View more

Hi Gludo and Elisabeth, I think you can use some Python to achieve this. See example below: def main():     import arcpy     # assume data is available as point featureclass     # and fc has projected coordinate system     fc = "c:/folder/yourGeodatabase.gdb/yourFCname" # change this (path to featureclass GPS points)     fldDistance = "DistGPSpoints"                   # change this (output name distance field)     # add distance field if it doesn't exist     if FieldExist(fc,fldDistance) == False:         arcpy.AddField_management(fc, fldDistance, "DOUBLE")     # Create update cursor for feature class     fields = ['SHAPE@XY', fldDistance]     with arcpy.da.UpdateCursor(fc, fields) as cursor:         X = 0         Y = 0         cnt = 0         for row in cursor:             cnt+=1             prevX = X             prevY = Y             X = row[0][0]             Y = row[0][1]             if cnt == 1:                 # for first point distance is 0                 distance = 0             else:                 distance = getDist(X,Y,prevX,prevY)             # Update row             row[1] = distance             cursor.updateRow(row)     del row def getDist(X1,Y1,X2,Y2):     import math     return math.hypot(X2 - X1, Y2 - Y1) def FieldExist(featureclass, fieldname):     import arcpy     fieldList = arcpy.ListFields(featureclass, fieldname)     fieldCount = len(fieldList)     if (fieldCount == 1):         return True     else:         return False if __name__ == '__main__':     main() You will need to change the path to the featureclass and the name of the output field. It assumes the data is available as point featureclass and the point should be projected (not in a geographic coordinate system). The script will: add the output distance field (now named "DistGPSpoints") if it doesn't exist yet loop through the points and calculate the distance to the previous point Kind regards, Xander ... View more

Is there are a way to perform OLS only on selected records/years from my feature layer? I am trying to analyze species occurrence over time. I have a "Year" field that denotes different years worth of data. I am using ArcGIS 10.1 SP1 on Win7 x64 Thanks! Hi Robert, If you have a lot of years in your data you may want to consider using Python: import arcpy, os from arcpy import env ws = "c:/Folder/yourFileGeodatabase.gdb" # change this fcName = "yourFC"                        # change this # fields in input FC fldYear = "Year" fldID = "yourUniqueIDfieldName"          # change this fldDepVar = "Rainfall"                   # change this: Dependent_Variable fldExpVar = "Temperature;..."            # change this: Explanatory_Variables Coefficient_Output_Table = "olsCoefTab" Diagnostic_Output_Table = "olsDiagTab" olsReport = "c:/Folder/OLSreport" # set workspace env.workspace = ws # create unique set of values in Year field fc = ws + os.sep + fcName values = [row[0] for row in arcpy.da.SearchCursor(fc, (fldYear))] uniqueYears = set(values) # loop through the values for year in uniqueYears:     # create feature layer for year     qDef = "{0} = {1}".format(fldYear,year)     fcFL = "Year_{0}".format(year)     arcpy.MakeFeatureLayer_management(fc,fcFL,qDef,"#")     # perform OLS on feasturelayer     fcOLS = "OLS_{0}".format(year)     ols = arcpy.OrdinaryLeastSquares_stats(fcFL, fldID, fcOLS, fldDepVar, fldExpVar, "{0}{1}".format(Coefficient_Output_Table,year),"{0}{1}".format(Diagnostic_Output_Table,year),"{0}{1}.pdf".format(olsReport,year)) This code will make a unique list of years, and then loop through the years, create a featurelayer based on that year and perform an OLS on it. The result for the year 2000 will be: OLS_2000 (OLS featureclass) olsCoefTab2000 (Coefficient Output Table) olsDiagTab2000 (Diagnostic Output Table) OLSreport2000.pdf (OLS report as PDF) I haven't tested the code. See for more explanation on how to define the required variables for OLS: Ordinary Least Squares (OLS) (Spatial Statistics) Kind regards, Xander ... View more

So maybe I should just drop the arcscripting? Hi Henri, I think it is a good idea to drop the arcgisscripting, since it requires loading the toolboxes and the location of those depend on the installation location, which makes it less interchangeable: # Load toolbox gp.AddToolbox("C:/ESRI/Desktop/Desktop10.1/ArcToolbox/Toolboxes/Spatial Analyst Tools.tbx") I think the code I attached in the previous post will work with 10.1 too (haven't tried though). Just change the lines: # Get the name of folder env.workspace = arcpy.GetParameterAsText(0) #Input DEM InDEM = arcpy.GetParameterAsText(1) #Output TWI raster NAME OutTWI = arcpy.GetParameterAsText(2) into: env.workspace = "G:\\xxxx\\LiDAR\\TWI\\" InDEM = 'dem_ldm.tif' OutTWI = "TWI01" # or a different name and see if it runs... Kind regards, Xander ... View more

I'll answer to myself: It seems that all the values are rounded to the nearest 0,5. So 0,401 m / 2.8284 m = 0,14177, would instead be 0,5/3 = 16,6667. Or if distance is 2 metres, then 0,5/2 = 25 %. In the manual it says the values are rounded in FLAT areas, which seems to be a little off considering it does that where slope is 25 %. Another remark is that in flat areas (lakes) the drop percentages are indeed continuous! [ATTACH=CONFIG]28643[/ATTACH] Here's a sample showing what I mean. Darker gray equals flow direcetion to lower left and percentage 16,666 % drop while lighter gray here has 25 % drop and flow direction to the left. Hi Henri, I had a look at the original script and your data noticed a few things: the script assumes a cell size of 5m you are working with a DEM with 2m cell size it uses the old arcgisscripting library The assumption of 5m cell size is based on the following lines in the script: gp.times_sa(Output_Plus_raster, "25", Output_Times_raster) # calculates area of cell 5x5=25m Input_true_raster_or_constant_value = "5" # distance between cell centers horizontally and vertically false_raster_or_constant_value = "7.07" # distance between cell centers diagonally I adopted the code to arcpy (10.2) and added some code to retrieve the cell size from the input DEM. No additional checks are added though. I used it on a DSM (including buildings and flat terrain giving high and low values) and couldn't reproduce the effects you describe. Give it a try and let me know how it went. Example usage: [ATTACH=CONFIG]28648[/ATTACH] #   Created By: Prasad A Pathak # Purpose: Calculate Topographic Wetness Index: natural logarithm of area divided by slope # It indicates the probable water saturation level of the ground # By default, ArcGIS calculates slope by considering just 3 X 3 neighborhood # For this index a slope is calculated from the topmost pixel to any particular pixel is needed # This is acheived by calculating the minimum elevation raster # small adaptions to ArcGIS 10.2 by Xander Bakker # - adapted to arcpy # - cellsize is determined from input DEM # - output TWI is name (string) not raster dataset # Import system modules import arcpy import os import math from arcpy import env env.overwriteOutput = True # Get the name of folder env.workspace = arcpy.GetParameterAsText(0) #Input DEM InDEM = arcpy.GetParameterAsText(1) #Output TWI raster NAME OutTWI = arcpy.GetParameterAsText(2) # note this was a raster, is now changed to a name # Check out license for Spatial Analyst arcpy.CheckOutExtension("Spatial") # determine cellsize of input raster desc = arcpy.Describe(InDEM) cellsize = desc.meanCellHeight # new # Variables used Output_surface_raster = env.workspace + os.sep + "eldodem_fill" Output_flow_direction_raster = env.workspace + os.sep  + "eldodem_flw" Output_drop_raster = env.workspace + os.sep + "eldodem_drp" Output_accumulation_raster = env.workspace + os.sep  + "eldodem_acc" Output_Plus_raster = env.workspace + os.sep  + "eldodempls" Output_Times_raster = env.workspace + os.sep  + "eldodemtim" Output_Mean_elev_raster = env.workspace + os.sep  + "mean_elev" Output_edrop = env.workspace + os.sep  + "edrop" Out_change = env.workspace + os.sep  + "change_elev" Out_distance = env.workspace + os.sep  + "distance" Out_slope = env.workspace + os.sep  + "slope" Out_preatan = env.workspace + os.sep  + "preatan" Out_TWI_raster = env.workspace + os.sep  + OutTWI # new # Fill DEM to make it depressionless filledDEM = arcpy.sa.Fill(InDEM) filledDEM.save(Output_surface_raster) print "DEM filled" arcpy.AddMessage("DEM filled") # Flow Direction flowDir = arcpy.sa.FlowDirection(filledDEM,"NORMAL",Output_drop_raster) flowDir.save(Output_flow_direction_raster) print "Flow direction raster created successfully" arcpy.AddMessage("Flow direction raster created successfully") # Flow Accumulation flowAcc = arcpy.sa.FlowAccumulation(flowDir,"","FLOAT") flowAcc.save(Output_accumulation_raster) print "Flow accumulation raster created successfully" arcpy.AddMessage("Flow accumulation raster created successfully") # One is added to each pixel to get an count of how many pixel including the current are contributing the flow flowAccPlus1 = flowAcc + 1 flowAccPlus1.save(Output_Plus_raster) print "addition raster created successfully" arcpy.AddMessage("addition raster created successfully") # Contributing Area using the number of pixels and flowAccPlus1Times = flowAccPlus1 * cellsize**2 flowAccPlus1Times.save(Output_Times_raster) print "multiplication raster created successfully" arcpy.AddMessage("multiplication raster created successfully") # Process: Block Statistics... meanElev = arcpy.sa.BlockStatistics(filledDEM, arcpy.sa.NbrRectangle(3, 3, "CELL"),"MINIMUM", "DATA") meanElev.save(Output_Mean_elev_raster) print "MIN elevation raster created successfully" arcpy.AddMessage("MIN elevation raster created successfully") eDrop = filledDEM - meanElev eDrop.save(Output_edrop) print "EDROP raster created successfully" arcpy.AddMessage("EDROP raster created successfully") Input_false_raster_or_constant_value = 0.005 outChange = arcpy.sa.Con(eDrop >= 0.005, eDrop, Input_false_raster_or_constant_value) outChange.save(Out_change) print "Change elevation raster created successfully" arcpy.AddMessage("Change elevation raster created successfully") inListras = arcpy.sa.InList(flowDir,[1,4,16,64]) Input_true_raster_or_constant_value = cellsize false_raster_or_constant_value = cellsize * math.sqrt(2) outDist = arcpy.sa.Con(arcpy.sa.IsNull(inListras),false_raster_or_constant_value,Input_true_raster_or_constant_value) outDist.save(Out_distance) print "Distance raster created successfully" arcpy.AddMessage("Distance raster created successfully") outSlope = outChange / outDist outSlope.save(Out_slope) print "Slope raster created successfully" arcpy.AddMessage("Slope raster created successfully") outPreatan = flowAccPlus1Times / outSlope outPreatan.save(Out_preatan) print "Pre-atan raster created successfully" arcpy.AddMessage("Pre-atan raster created successfully") outTWI = arcpy.sa.Ln(outPreatan) outTWI.save(Out_TWI_raster) print "TWI raster created successfully" arcpy.AddMessage("TWI raster created successfully") Kind regards, Xander ... View more

Hi, I'm trying to compute the least cost distance between a source (point A) and several destinations (points B) using as cost surface the slope of a DEM. I've reclassified the slope in 9 categories, and I do not have any value equal to 0 or missing. To compute the cost distance raster I used point A as the source layer, and I limit the maximum distance to 10km. I've computed the cost distance raster directly using spatial analyst toolbox and using python, in both cases when I display the cost distance raster in ArcmMap it appears displaced with respect the source: Point A is not in the center of the raster, but it appears in the border with a high cost associates to its location. Boths layers: source and reclassified slope have the same projection WGS1984_UTM16. I've tried to do it converting the point to a raster, then the pixel appears again displaced far away from the center of the cost distance raster. I would like to know if anyone has any clue of what could be wrong. I attach a image of the cost distance raster, the red flag represents the point A or source and the red circle is a 10km buffer around point A. Thanks for your comments.[ATTACH=CONFIG]28629[/ATTACH] Hi Teresa, This type of error is normally caused by different projections and the lack of a proper transformation. Since you indicate that all your data use the same projection, maybe there are other data with different projections in your TOC (check the dataframe properties) or your environment settings contain a different setting for output coordinate system. Could you check these settings? Kind regards, Xander ... View more