Latest Contributions by XanderBakker

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

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 ... View more

Although I don't think this has anything to do with arcpy / ArcGIS, there is a module you could load called calendar: import calendar print calendar.month(2014,2) results in:    February 2014 Mo Tu We Th Fr Sa Su                                 1  2 3   4   5   6  7   8   9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 sorry, but the formatting changes the output. So I added an image: [ATTACH=CONFIG]31572[/ATTACH] See for more functionality: http://docs.python.org/2.7/library/calendar.html?highlight=calendar#calendar Kind regards, Xander ... View more