Latest Contributions by XanderBakker

What then does the documentation refer to ? You are completely right. I does not seem to do what is described in the help. Even the slightly change sample script does not return what one should expect based on the description. import arcpy fc = r'C:\Project\_Forums\ZonStatPol\test.gdb\Intersect' feature_count = int(arcpy.GetCount_management(fc).getOutput(0)) if feature_count == 0:     arcpy.AddError("{0} has no features.".format(fc)) else:     arcpy.AddMessage("{0} has {1} features.".format(fc, feature_count)) print arcpy.GetMessages() Results in: Executing: GetCount C:\Project\_Forums\ZonStatPol\test.gdb\Intersect Start Time: Tue Feb 25 15:46:45 2014 Row Count = 10 Succeeded at Tue Feb 25 15:46:45 2014 (Elapsed Time: 0,19 seconds) ... which is the output of the geoprocessing tool GetCount_management Said, but true (or I simply don't get it). Kind regards, Xander ... View more

'Challenge' 1: import arcpy arcpy.env.workspace = r'C:\Project\999999\myFGDB.gdb'  # example FileGeodatabase as workspace arcpy.env.workspace = r'C:\Project\999999\mySHPfolder'  # example folder (with shapefiles) as workspace fclasses = arcpy.ListFeatureClasses() for fclass in fclasses:     fctype = arcpy.Describe(fclass).shapeType     print "{0} is a {1} feature class".format(fclass, fctype) 'Challenge' 2: import arcpy, os arcpy.env.workspace = r'C:\Project\999999\mySHPfolder' fclasses = arcpy.ListFeatureClasses() out_ws_all = r'C:\Project\999999\test_all.gdb' out_ws_pol = r'C:\Project\999999\test_pol.gdb' for fclass in fclasses:     fctype = arcpy.Describe(fclass).shapeType     newname = fclass     if newname[-4:].upper() == ".SHP":         newname = newname[:-4] # trim the .shp extension from the name     out_fc = os.path.join(out_ws_all, fclass)     arcpy.CopyFeatures_management(fclass, out_fc)     if fctype == 'Polygon':         out_fc = os.path.join(out_ws_pol, fclass)         arcpy.CopyFeatures_management(fclass, out_fc) Kind regards, Xander ... View more

Hi Jollong, In case of going for an intersect, you should start with converting the raster to polygons. Make sure that the "Simplify polygons" option is switched off. Intersect this polygon featureclass with your zones (circles). The resulting polygon featureclass will have areas containing information for each (part of a) pixel within each zone (circle). Since you know the area of the circles you can calculate for each feature (combination of zone and pixel) the fraction of the circle area. If the area of the circle varies, you can join it from the zones to the intersect result. See the image below: [ATTACH=CONFIG]31721[/ATTACH] The ZoneArea is obtained by joining the original area of the zones and the percentage is calculated by dividing the Shape_Area by the ZoneArea (* 100). Please note that some combinations of zones and pixels occur twice in the table. These is caused by the overlap between the zones (circles). You should aggregate the results to obtain a more useful result. Kind regards, Xander ... View more

Hi Jollong, In that case you should divide your zones featureclass into multiple featureclass without any overlap and perform the zonal statistics for each featureclass. The alternative is to convert your raster to polygons and perform an intersect with your circles. If I'm correct the overlapping circles will not be a problem in the vector analysis. Kind regards, Xander ... View more

Hi Tony, The arcpy.AddMessage() will add a message that is visible if you run your script as a tool from a toolbox. The arcpy.GetMessages() will return the messages from a geoprocessing tool used. In your example script you don't run any tool, so no message is returned. If I create a simple script: import arcpy arcpy.AddMessage("Something...") arcpy.AddMessage("Blah...") cnt = 0 for message in arcpy.GetMessages():     cnt += 1     arcpy.AddMessage("Message {0} is '{1}'".format(cnt, message)) arcpy.AddMessage("arcpy.GetMessageCount() results in: {0}".format(arcpy.GetMessageCount())) The output in the tool window will be: Executing: Script Start Time: Tue Feb 25 08:20:02 2014 Running script Script... Something... Blah... arcpy.GetMessageCount() results in: 0 Completed script Script... Succeeded at Tue Feb 25 08:20:02 2014 (Elapsed Time: 0,00 seconds) Kind regards, Xander ... View more

Hi Nick, If the address is all you have, it will not be possible. If your address location has information on population density (X pop/km2), then you could create a buffer around the point: radius = SQRT((30.000 pop / (X pop/km2)) / PI) * 1000 (in meters) If you have building blocks (or center points of building blocks), you could use some code to sort the distance of each block (only use block that are within a certain distance) to the address location and sum the population per building block until you reach the 30K. Merge the polygons together or create a convex hull to obtain the polygon you're looking for. ... or maybe Grouping Analysis works for you. Kind regards, Xander ... View more

If you look at the help topic on "Zonal Statistics (Spatial Analyst)", you will find the following text: If the zone input is a feature dataset, a vector-to-raster conversion will be internally applied to it. To ensure that the results of the conversion will align properly with the value raster, it is recommended that you check that the extent and snap raster are set appropriately in the environment settings and the raster settings. It is possible that when you manually convert the blue circles to raster (using the existing raster as snap raster), that the number of pixels per circle may be slightly different due to their relative position to the pixels of the other raster. Kind regards, Xander ... View more

Hi Mark, Although this type of analysis is normally done using the raster format, it can be done using a featureclass (shapefile). Assuming that you are interested in an angle between the start point and the end point of each feature you could do this: open the attribute table of your shapefile add a new field that will hold the cardinal flow direction, make sure it is a text field (let's call this field "CardFD") right click on your new field "CardFD" and select "Field Calculator..." to open the Field Calculator Choose Python as parser Switch the "Show Codeblock" on Paste the code below in the "Pre-Logic Script Code": import arcpy from math import atan2, pi def getCardinal(angle):     lstCardinal = [[22.5,"E"], [67.5, "NE"], [112.5, "N"],                    [157.5, "NW"], [202.5, "W"], [247.5, "SW"],                    [292.5, "S"], [337.5, "SE"], [360, "E"]]     for item in lstCardinal:         value = item[0]         if angle < value:             cardinal = item[1]             break     return cardinal def calcGeomCardinality(polyline):     pnt1 = polyline.firstPoint     pnt2 = polyline.lastPoint     angle_deg = (atan2(pnt2.Y - pnt1.Y, pnt2.X - pnt1.X)) * 180.0 / pi     if angle_deg < 0:         angle_deg = 360 + angle_deg     return getCardinal(angle_deg) Paste the line below in as formula (just below "CardFD =") calcGeomCardinality( !SHAPE!) Click "OK" This will take each polyline feature and extract the from and to node from it. Next it will calculate the angle and parse it to cardinal flow direction. You can also write the angle itself to a new field: add a new field (double or  float) that can hold the angle between the start and end point . Let's call this field "FlowAngle" right click on your new field "FlowAngle" and select "Field Calculator..." to open the Field Calculator it will probably show the code from cardinal flow direction. Replace the code in the "Pre-Logic Script Code" by import arcpy from math import atan2, pi def getAngle(polyline):     pnt1 = polyline.firstPoint     pnt2 = polyline.lastPoint     angle_deg = (atan2(pnt2.Y - pnt1.Y, pnt2.X - pnt1.X)) * 180.0 / pi     if angle_deg < 0:         angle_deg = 360 + angle_deg     return angle_deg Paste the line below in as formula (just below "FlowAngle=") getAngle( !SHAPE!) Click "OK" If you are going to perform some analysis to determine the predominant angle, please be aware that the resulting angle does not have to be representative for the area. Kind regards, Xander ... View more

Hi Holly, I recently did an open source project with Python and found the website below very helpful: http://www.lfd.uci.edu/~gohlke/pythonlibs/#gdal I installed the executable "GDAL-1.10.1.win32-py2.7.exe" and was able to use "gdal, ogr and osr" without any further pip or easy_install. Please pay attention to the note on that page: This distribution includes all required files. Do not use together with OSGeo4W or gdalwin32. Kind regards, Xander ... View more

Hi Robert, I tried the code and resulted with the same problems as you encountered; some polygons have a zero area and all that draw are overlapping exactly. I thought that this might be caused by a lacking coordinate system and even tried to project the polygon to a projected coordinate system. After this I notices that both seem to work. See code below: import arcpy arcpy.env.overwriteOutput = True infile = r'C:\Project\_Forums\CSV2Polygon\csv\Metadata_RP.csv' outshp_gcs = r'C:\Project\_Forums\CSV2Polygon\shp\poly_gcs.shp' outshp_pcs = r'C:\Project\_Forums\CSV2Polygon\shp\poly_pcs.shp' sr_in = arcpy.SpatialReference(4283) # GCS_GDA_1994 sr_out = arcpy.SpatialReference(28355) # GDA_1994_MGA_Zone_55 lstPolsGCS = [] lstPolsPCS = [] with arcpy.da.SearchCursor(infile, ["o_minx","o_miny","o_maxx","o_maxy"]) as curs:     for row in curs:         xvalues = [row[0], row[2]]         yvalues = [row[1], row[3]]         print "{0}, {1}, {2}, {3}".format(row[0], row[1], row[2], row[3])         pnt1 = arcpy.Point(min(xvalues),min(yvalues))         pnt2 = arcpy.Point(min(xvalues),max(yvalues))         pnt3 = arcpy.Point(max(xvalues),max(yvalues))         pnt4 = arcpy.Point(max(xvalues),min(yvalues))         polygon_gcs = arcpy.Polygon(arcpy.Array([pnt1, pnt2, pnt3, pnt4, pnt1]), sr_in)         lstPolsGCS.append(polygon_gcs)         polygon_pcs = polygon_gcs.projectAs(sr_out)         lstPolsPCS.append(polygon_pcs) arcpy.CopyFeatures_management(lstPolsGCS, outshp_gcs) arcpy.CopyFeatures_management(lstPolsPCS, outshp_pcs) In my case I create two outputs; both polygon shapefiles, one for the geographic coordinates and one for the projected coordinates. Essential is providing the proper spatial reference when creating the geometry (see sr_in). I create a new polygon with the projected copy of the input polygon. Both are added to a list and written to their corresponding shapefiles. BTW; the reason I first add the x and y coordinates of a row to lists, it to assure to get the minimum and maximum values and create polygons that are clockwise. There is one record (polygon # 10) where the xmin > xmax. o_minx    o_miny    o_maxx    o_maxy 144,8203  -37,7621  144,8201  -37,761 Kind regards, Xander ... View more

Received by email: Just I inform you that I tried to combine in a unique raster named Constraints: - The water raster obtained from rivers and lakes layers (having the value 1). - Roads raster which has the value 100. - NonForest raster which has the values 1 and 0 (1 for NonForest pixels and 0 for Forest pixels). So I used Raster Calculator: [Constraints]=[water] + [Roads]+[NonForest] The result was incomprehensible raster representing points scattered in space and having the values 101 and 102. I found the same result when I tried to create the following raster (Raster Calculator): [Soil Loss]=[Factors]-[Constraints]. First of all I would like to ask you not to send me data to my work email. However, I did take a short look at it and I notice the following things: Assuming that the water raster is derived from "riviere" and "lac", if you zoom in close you will notice that there is a shift in pixel location. Please read the Help topic on "How the Snap Raster environment setting works" and make sure your output extent and snap raster is always the same. Make sure that you use one pixel size along your analysis (I see 82m for the water raster, and 50m for the forest raster). The water raster contains 1 for water and NoData for no water. If you sum a value (known value) with a NoData value, the result will be NoData. NoData is not the same as 0. Please read the topic on "NoData and how it affects analysis". So, if NoData means 0, make sure you assign a 0 value to the NoData pixels before combining it with other rasters or handle NoData appropriately. Kind regards, Xander ... View more

Nice!  I pulled from one of our impmentations with lots of pandas DataFrame processing, but I like this suggestion for just numpy work. Hi James, it's posts like yours (thinking outside the box) that lets me learn something new every day... (+1 for that). Kind regards, Xander ... View more

Hello, I have two layers : Rivers (as a polyline shapefile) and lakes (as a polygon shapefile). I want to create a new layer of proximity factor to the streams (Rivers and lakes) which varies by intervals from 5 (near the streams) to 1 (the most distant zones). I used Merge and Union tools. Also I transformed these layers to Raster and I used Euclidean Distance in ArcGIS 10.2... but these solutions didn't work. Can you show me which solution is the most appropriate to this study case ? Thank you in advance Hi Jallouli, There are different solutions to this problem. In case you do this the raster way (Spatial Analyst), you will have to convert the rivers and lakes to rasters (separately). Before doing that it maybe wise to add a new field and fill it with the value 1. Next you will have to create a new raster combining the "rivers" and "lakes" rasters using a Con statement that looks something like: water = Con(IsNull([rivers]), [lakes], [rivers]) After this you can use the water raster in the Euclidean Distance tool and finish with a reclassify to create the intervals. It is easier though to do this in vector format. Since it is not possible to merge different geometry types into the same featureclass, you should choose a geometry type that best represents both types of data (lakes and rivers). Since a river has a surface, you can buffer the rivers with a small value (best to choose half the average river width). This results in a polygon featureclass. Merge the buffered rivers with the lakes to create a single polygon featureclass (let's call this "water"). Now buffer the "water" polygons with the "Multiple Ring Buffer" tool and specify the distances you wish. This method is more accurate and less work. Kind regards, Xander ... View more

While we're at it, using part of the idea James suggested (the numpy part) and mixing it with some list comprehensions I came up with this: import arcpy, numpy intbl = r'C:\Project\_Forums\CommaDelimitField\test.gdb\intable' outtbl = r'C:\Project\_Forums\CommaDelimitField\test.gdb\tst01' flds = ('number_','Code','Year_') # list comprehensions lst_in = [row for row in arcpy.da.SearchCursor(intbl, flds)] lst_out = [(num, row[1], row[2]) for row in lst_in for num in row[0].split(', ')] # use numpy to store the table npa = numpy.array(lst_out, numpy.dtype([('number', '|S10'), ('code', '|S10'), ('year', numpy.int32)])) arcpy.da.NumPyArrayToTable(npa, outtbl, ("number", "code", "year")) Kind regards, Xander ... View more