|
POST
|
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
02-24-2014
10:12 PM
|
0
|
0
|
1670
|
|
POST
|
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
02-24-2014
09:21 PM
|
0
|
0
|
1232
|
|
POST
|
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
02-24-2014
05:16 AM
|
0
|
0
|
810
|
|
POST
|
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
02-24-2014
12:15 AM
|
0
|
0
|
1670
|
|
POST
|
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
02-24-2014
12:05 AM
|
0
|
0
|
3196
|
|
POST
|
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
02-23-2014
09:57 PM
|
0
|
0
|
1259
|
|
POST
|
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
02-20-2014
10:21 PM
|
0
|
0
|
1272
|
|
POST
|
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
02-20-2014
08:48 PM
|
0
|
0
|
3539
|
|
POST
|
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
02-20-2014
03:00 AM
|
0
|
0
|
2232
|
|
POST
|
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
02-19-2014
09:56 PM
|
0
|
0
|
3539
|
|
POST
|
Hi Rafael, Is the error message pointing to a specific line in the code? What is the exact content of the parameters? If I understand this correctly you create a feature layer and for each row you create a selection of that row, using the FID. I suppose there is more code which is not included in your post. What do you want to achieve with your code (what happens after the selection is created)? A few pointers: In your code you concatenate the path with a name "territorios_lyr". If your path is a true path, the result will be missing the trailing slash at the end of the path. It is better to use os.path.join(path, "territorios_lyr") , this will require import os However, when you create a feature layer, it resides in memory. You only specify a name (not a path) when creating the feature layer, so change the code to: territory_lyr = "territorios_lyr" It is better to use "format" when combining text and numbers (for building the where clause): where = "FID={0}".format(currentID) Since the where clause depends on the type of data (with respect to field delimiters) and since the name of the FID column can be derived using arcpy.Describe, you could change it to: fldFID = arcpy.Describe(territory).OIDFieldName fields= (fldFID, "avg_mindis", "I_conc", "Max_lenght","TamanoTerr") with arcpy.da.UpdateCursor(territory, fields) as cursor: for row in cursor: currentID = row[0] where = "{0}={1}".format(arcpy.AddFieldDelimiters(territory, fldFID), currentID) arcpy.SelectLayerByAttribute_management(territory_lyr, "NEW_SELECTION", where) # do something with the selection... Kind regards, Xander
... View more
02-19-2014
09:21 PM
|
0
|
0
|
1243
|
|
POST
|
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
02-19-2014
01:32 PM
|
1
|
0
|
2232
|
|
POST
|
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
02-19-2014
12:27 PM
|
0
|
0
|
1948
|
|
POST
|
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
02-19-2014
02:29 AM
|
0
|
0
|
1318
|
|
POST
|
Hi Amy, I think this will do the job: import arcpy, os # access your table(s) intbl = r'C:\path\to\your\FileGeodatabase.gdb\inputtable' outtbl = r'C:\path\to\your\FileGeodatabase.gdb\outputtable' # optionally create output table, using the input table as schema # CreateTable_management (out_path, out_name, {template}, {config_keyword}) outpath, outname = os.path.split(outtbl) arcpy.CreateTable_management(outpath, outname, intbl) # have a close look at the fieldnames, some names may be restricted # and are changed by ArcGIS when imported from for instance Excel flds = ['number','Code','Year'] # do a search cursor on the input table and an insert cursor on the output table outcurs = arcpy.da.InsertCursor(outtbl, flds) with arcpy.da.SearchCursor(intbl, flds) as incurs: for row in incurs: numbers = row[0].split(',') for number in numbers: # for each 'number' a record is written outcurs.insertRow((number.strip(), row[1], row[2])) del outcurs Kind regards, Xander
... View more
02-19-2014
02:11 AM
|
0
|
0
|
2232
|
| Title | Kudos | Posted |
|---|---|---|
| 6 | 12-20-2019 08:41 AM | |
| 1 | 01-21-2020 07:21 AM | |
| 2 | 01-30-2020 12:46 PM | |
| 1 | 05-30-2019 08:24 AM | |
| 1 | 05-29-2019 02:45 PM |
| Online Status |
Offline
|
| Date Last Visited |
11-26-2025
02:43 PM
|