|
POST
|
Hello, I am trying create a raster that shows road density (m/ha) within 3km radius of the cell. The line density tool in Spatial Analyst appears to be the tool to do this as it allows for input of a line file (roads), a search radius (3000m), output raster cell size (10m), and area unit (ha). However, upon running the tool I getting results that I don't believe to be correct. For example the maximum value I am getting is 0.27 m/ha, which is way to low of a value considering to the high road density in the area. To double check I converted a portion of the raster to points and buffered those points by 3 km and found the area of road within this buffer. Using this information the road density values come out to around 9 or 10 m/ha which is much more reasonable for the density of roads in the area. The dataframe and road layer is in UTM Zone 12N so the units should be in meters. Any ideas why these results are not correct? I wonder if it has to do with the population field being set to none as in the tool help it indicates the population field is a: Numeric field denoting population values (the number of times the line should be counted) for each polyline. Values in the population field may be integer or floating point. The options and default behaviours for the field are listed below. Use None if no item or special value will be used and each feature will be counted once. Therefore I'm wondering if it is only counting each road line feature the first time it intersects a cell search radius instead of counting it within the search radius of each cell. Any advice or suggestions would be greatly appreciated. Below is my inputs if that helps. Thanks in advance. Hi Karl, If you look at the explanation in the Help you will see that the density is calculated as follows: Density = ((L1 * V1) + (L2 * V2)) / (area_of_circle) Lines L1 and L2 represent the length of the portion of each line that falls within the circle. The corresponding population field values are V1 and V2. Source: http://resources.arcgis.com/en/help/main/10.2/index.html#//009z00000012000000 In your case, you are not applying a population field, so it is reduced to: Density = (L1 + L2) / (area_of_circle) Since you are using a search radius of 3000m, the area of the circle is ± 28274333m² (or 2827ha or 28.27km²). This causes the density to result in very low values. In Help you will also find this: A default area unit is selected based on the linear unit of the projection definition of the input polyline feature data or as otherwise specified in the output coordinate system environment setting. For line density, when an area unit factor is specified, it converts the units of both length and area. For example, if the linear unit is meters, the output area units will default to SQUARE_KILOMETERS and the resulting line density units will convert to kilometers per square kilometer. The end result, comparing an area scale factor of meters to kilometers, will be the density values being different by a multiplier of 1,000. (... so not 1,000,000) Hope this helps you to interpret the results, although it may not be what you are looking for... Kind regards, Xander
... View more
09-05-2013
01:31 AM
|
0
|
0
|
784
|
|
POST
|
Awesome idea Xander. Thank you. Although I only have 8 polygons to extrude over the whole USA. So mine won't look as cool as yours. Plus some of my polygons with small values are behind some with much larger values. It would cover the stack. I could create a new polygon shapefile with just squares at the center of each large region, then just join the region $ data to each of these new polygons. Then extrude just this layer so as not to have anything overlapping. And maybe extrude the regional polygons negatively so that it doesn't cover up the bar. Or I could export the tilted map with the regions shaded to a 2D png/jpg and insert it into a PowerPoint slide and insert each 3D bar graph over the regions with the original excel file. Its easy for my purposes and looks pretty cool for not having an expensive graphics program. I'll try both and see what goes over better with the boss. Thanks again. Hi Jessica, You have enough creative ideas to make this work. Maybe you could share the result on this forum to inspire others ... Kind regards, Xander
... View more
09-04-2013
10:19 PM
|
0
|
1
|
2435
|
|
POST
|
Is this in the process of being fixed? It would be nice to make maps like this http://www.researchgate.net/post/How_do_you_create_this_visual_pseudo-3d_map_representation without having to buy and learn a graphics program. Hi Jessica (and Jonathan), Not sure if this issue was fixed, but with the current functionality you can do some interesting visualizations too. For a project I wanted to show the spatial distribution of ammoniac emission by farm animals in the Netherlands and I choose to do this in ArcScene (10.1 at the time). This was the result (extracted from a presentation): [ATTACH=CONFIG]27190[/ATTACH] What I had was 4 attribute columns with the emission values per square (2.5x2.5 km²). I added this layer 4 times to ArcScene (or copy and paste the layer in the TOC). For the first emission (species) set the Base Height to 0 and the Extrusion to the attribute value itself. For the second emission (species), set the Base Height to the first emission value and the Extrusion to the attribute value itself and so on. The last one starts (Base Height) at the sum of the other emissions... BTW; I just tried to use the stack graph in 10.2 and the data still seems to disappear. Kind regards, Xander
... View more
09-03-2013
11:22 PM
|
0
|
0
|
2435
|
|
POST
|
Hi, I�??d like to normalize map data in raster format to the data range using Raster Calculator. I'm not so much concerned with the accuracy of the method. e.g. v�??(i) = (v(i) - minA) / (maxA - minA) What is the syntax I need to do this? I am getting problems using the following script: ("raster" - min("raster")) / (max("raster") - min("raster")) Thanks! Hi, If I could convince you to do this in the Python Window, it would be a lot easier. In that case you can make use of the Raster object in arcpy. Just have a look at the code below: import arcpy
my_raster = arcpy.Raster(r'obj\objecten')
result = (my_raster - my_raster.minimum) / (my_raster.maximum - my_raster.minimum)
Basically you load the arcpy module, you define a variable "my_raster" to which you will assign the raster you want to process. In this case I used a raster called "objecten" which resides in a grouplayer called "obj", but you can just specify the name of the raster as it is called in your table of content. The last line performs the calculation (make sure you have the Spatial Analyst extension switched on). As you can see you can access the minimum and maximum values of a raster, since they are properties of the raster object. See also: http://resources.arcgis.com/en/help/main/10.1/index.html#//002100000017000000 http://resources.arcgis.com/en/help/main/10.1/index.html#//018z00000051000000 To save the result, you can add another line:
result.save('yourNameGoesHere')
It will save it in your current workspace. Hope this helps you. Kind regards, Xander
... View more
09-03-2013
04:22 AM
|
0
|
0
|
956
|
|
POST
|
Hi Yang, I suppose you have experience with ArcGIS 9.x and are used to find the functionality of 3D analyst in the 3D Analyst toolbar. In 9.x you would see something like this: [ATTACH=CONFIG]27167[/ATTACH] All of these tools (and more) can be accessed now through the ArcToolbox: [ATTACH=CONFIG]27168[/ATTACH] This screenshot is from 10.2, but you will see something similar at 10.1. Hope this is what you are looking for... Kind regards, Xander
... View more
09-03-2013
03:41 AM
|
0
|
0
|
2804
|
|
POST
|
Hi Zuoqi Chen, Having a dataset locked can be a difficult thing to solve. IMHO I thing the locking mechanism in ArcGIS has been implemented to drastic, causing more problems than it avoids... I assume you are writing your raster datasets to a file geodatabase. It is possible to test the scema locking before deleting and place it in a "try except" (as shown below), but this doesn't solve anything if you have a lock on your dataset. Always be sure to delete any references to your dataset and delete cursors, rows, etc...
if arcpy.TestSchemaLock(data):
try:
arcpy.Delete(data)
except:
# catch error here
else:
print "Unable to acquire schema lock"
source: http://resources.arcgis.com/en/help/main/10.2/index.html#/TestSchemaLock/018v0000002m000000/ What I have also come across is using the garbage collector. This may sometimes solve the problem. See snippet below:
# at beginning of script
import gc # Garbage Collector
gc.collect()
# at end of script
collected = gc.collect()
print "Garbage collector: collected %d objects." % (collected)
# try to delete now...
sources: http://gis.stackexchange.com/questions/19408/arcgis-10-0-python-searchcursor-file-locking http://docs.python.org/2/library/gc.html http://pymotw.com/2/gc/ If this doesn't work, you could go back to the gold old Esri grid data format (write raster to folder), which is probably one of the last formats for which the locking mechanism hasn't been implemented. 😄 This did solve some nasty locking problems I had...
... View more
09-03-2013
01:17 AM
|
0
|
0
|
1037
|
|
POST
|
Hi yang li, Judging from your first screenshot, the 3D extension is already switched on (otherwise the ArcScene and ArcGlobe icons would have been grayed out). To enable the icons (functionality) on the 3D Analyst toolbar, you need to add data to your table of content. When you add a valid raster you will see this: [ATTACH=CONFIG]27159[/ATTACH] Most options are enabled, only the Profile Graph option remains disabled, since it requires a selection of 3D lines. Kind regards, Xander
... View more
09-02-2013
11:22 PM
|
0
|
0
|
2804
|
|
POST
|
Hi Peter, I'm glad that worked out for you. Could you click on the green check mark? See also the note below: Did you post and someone replied with a good answer? Help the community to find these great posts with answers by clicking the green check mark to indicate that it has been answered for you. Kind regards, Xander
... View more
09-02-2013
10:12 PM
|
0
|
0
|
853
|
|
POST
|
Hi Peter, In fact there is (at least in 10.1 and 10.2, I didn't check older versions). The arcpy.mapping.MapDocument has a property "isDDPEnabled": import arcpy myMXD = r'C:\Project\_LearnPython\DDPwelaan.mxd' mxd = arcpy.mapping.MapDocument(myMXD) if mxd.isDDPEnabled: # do something with Data Driven Pages pass else: # do something else here pass Links: http://resources.arcgis.com/en/help/main/10.1/index.html#//00s30000000n000000 http://resources.arcgis.com/en/help/main/10.2/index.html#//00s30000000n000000 http://resources.arcgis.com/en/help/main/10.2/index.html#/DataDrivenPages/00s300000030000000/
... View more
09-02-2013
05:19 AM
|
0
|
0
|
853
|
|
POST
|
Hi Stuart, Below you will find a snippet of Python code that can be used to calculate the maximum M for each feature in a featureclass. It is a standalone script that can be executed from your Python IDE of choice. Before you do, you should change the location/name of the featureclass used and the name of the (existing) field where the result will be written. In short it will: access your featureclass loop through the features calculate per feature the maximum M, using GetMmax() update the feature (row) def main():
import arcpy
# settings, replace with your own dataset and existing column name:
FClassOrTable = r'C:\Project\_tmp\test.gdb\Mawares'
outputField = 'MMax'
# .. or make variable
# FClassOrTable = arcpy.GetParameterAsText(0)
# outputField = arcpy.GetParameterAsText(1)
# update featureclass
fields = ['SHAPE@',outputField]
with arcpy.da.UpdateCursor(FClassOrTable, fields) as cursor:
for row in cursor:
row[1] = GetMmax(row[0])
cursor.updateRow(row)
del row
del cursor
def GetMmax(geom):
# get the max(M) from a geometry
mmax = -9999
for part in geom.getPart():
for point in part:
m = point.M
if m > mmax:
mmax = m
return mmax
if __name__ == '__main__':
main()
You can also add the script to a toolbox. The tool can eventually be used in modelbuilder in case you want to combine it with additional tools. In case you want to create a tool from it, you should uncomment the "arcpy.GetParameterAsText()" lines and define the parameters properly. first: the location/name of the featureclass (Type: Feature Layer) second: the field that will store the max M (Type: Field, obtained from the featureclass, and apply a fields filter to only show the double or floating fields). Links to additional information: http://resources.arcgis.com/en/help/main/10.1/index.html#/UpdateCursor/018w00000014000000/ Cheers, Xander
... View more
09-02-2013
04:48 AM
|
0
|
1
|
987
|
|
POST
|
Hi Ulises, I notice that you have a friction-surface that contains the estimated time to cross a pixel: 2) a friction-surface raster where each pixel contains the estimated time (in seconds) required to cross that pixel (based on different travel speeds according to different land classes that include roads). The resolution of my rasters is 30m by 30m. You should however express the friction in time per unit (sec/m) in your case. See the Help section: http://resources.arcgis.com/en/help/main/10.1/index.html#/How_the_cost_distance_tools_work/009z00000025000000/ The cost assigned to each cell represents the cost per unit distance for moving through the cell. The final value per cell is the cell size multiplied by the cost value. For example, if the cost raster has a cell size of 30, and a particular cell has a cost value of 10, the final cost of that cell is 300 units. To correct this, simply divide your friction raster by the pixelsize (30) and try again. Good luck!
... View more
08-30-2012
12:03 AM
|
0
|
0
|
644
|
|
POST
|
Hi Ryan, Thanks for the explanation, but I still don't subscribe to this point of view... In my stubborn intent to combine the Hillshade command with the save command, something unexpected happened. For the record; I didn�??t have my hopes up high, but wanted to know how ArcMap would respond. I loaded a DEM called myDEM to ArcMap, created a new file geodatabase and set it as my default geodatabase. After activating the Spatial Analyst extension, I opened the Python window and started with the line: from arcpy.sa import * My Hillshade command would be as simple as: myHS = Hillshade("myDEM",315,45) To force the temporary raster to be saved as �??myHS�?�, I started to type �??save(�?� at the end of the command. The moment the left bracket is typed, the command is executed. Every time you type another character the command is invoked again and again. In some cases an error about the background geoprocessing was thrown, and in another case ArcMap threw an error and closed�?� Now, the nice thing is, that if you type the command in (say Notepad) and copy and paste it in the python window, the command is executed correctly and stores the raster with the given name: myHS = Hillshade("myDEM",315,45).save("myHS") The only strange behavior is that the resulting raster is not added to the TOC (although this is indicated in the Geoprocessing options) OK, so I shouldn't do this, but is this the way ArcMap should react to such a situation? In the old days when I used to work with ILWIS, they implemented two types of command lines. Output = command and Output := command The first resulting in a virtual object which would be (re)calculated in the moment it was requested and would change if the input data would change and the latter creating the physical output raster with the name provided. I wonder if it would be possible to implement something like that (meaning := creates the raster with the name indicated, using only the = sign results in a temporary dataset) . Why is there a difference between arcpy commands and arcpy.sa commands? If I use the arcpy.Clip_management command I don�??t use an output = command construction, but I am able to define the output name: arcpy.Clip_management("myDEM","1952602 294196 1953546 296176","myClip", "#", "#", "NONE") In the command above I am able to define the output name of the raster (it will be saved as "myClip")�?� Kind regards,
... View more
11-07-2011
04:00 AM
|
0
|
0
|
976
|
|
POST
|
I tried reproducing this but could not. nbrCircle, nbrWeight, nbrIrregular, etc. produced no results other than to open that particular help topic. Thanks for the effort! Thanks for testing. It does sound like it is most likely specific to Windows XP. Hi Rhond, any luck reproducing the bug, or is it just me? Kind regards,
... View more
10-28-2011
12:14 AM
|
0
|
0
|
309
|
|
POST
|
Workspace and scratchWorkspace doing strange things when using arcpy.sa? An odd thing just happened. I was working on some analysis using arcpy.sa when a simple Focal Statistics: >>> AHN2CFILL2FM2 = FocalStatistics("AHN2CFILL2",NbrRectangle(3,3,"CELL"),"MEAN","DATA") â?¦ resulted in the following error: Runtime error File <string>, line 1, in <module> File c:\program files\arcgis\desktop10.1\arcpy\arcpy\sa\Functions.py, line 4788, in FocalStatistics File c:\program files\arcgis\desktop10.1\arcpy\arcpy\sa\Utils.py, line 47, in swapper File c:\program files\arcgis\desktop10.1\arcpy\arcpy\sa\Functions.py, line 4782, in wrapper File c:\program files\arcgis\desktop10.1\arcpy\arcpy\geoprocessing\_base.py, line 490, in <lambda> <class 'arcgisscripting.ExecuteError'>: ERROR 000875: Output raster: C:\Project\_Wetterskip\data\Wetterskip.gdb\C:\Project\_Wetterskip\data\py\py\Wetterskip.gdb\FocalSt_AHN21's workspace is an invalid output workspace. ERROR 000581: Invalid parameters. Failed to execute (FocalStatistics). Apparently, for some reason my workspace was set to a non valid location. When I listed the workspace and scratchworkspace I noticed this: >>> print env.workspace C:\Project\_Wetterskip\data\Wetterskip.gdb >>> print env.scratchWorkspace C:\Project\_Wetterskip\data\py\py\Wetterskip.gdb My workspace was still defined correctly, but for some strange reason my scratch workspace (normally the same location as my workspace) had been set to a non existing subfolder \py\py, where obviously it wasnâ??t going to find any geodatabaseâ?¦ But why would my output workspace been set to a combination of my workspace and my scratchWorkspace? The only thing that I can imagine being the reason for this, is that shortly before executing the Focal Statistics command I had used the option "Save Asâ?¦" from the context sensitive menu, shown when clicking the right mouse button in the python window. On saving I had created a new folder "py" (but only 1 "\py" subfolder and not "\py\py"). Also, when opening the Focal Statistics tool in the toolbox, the default output location was pointing to the non existing location (see Screenshot_FocalStatisticsTool). Investigation the Environment setting (through Geoprocessing menu, option Environmentsâ?¦) also showed the wrong location for scratch workspace. I closed the tool and the environment setting, went back to my python window and set my scratchWorkspace to my current workspace: >>> env.scratchWorkspace = env.workspace >>> print env.scratchWorkspace C:\Project\_Wetterskip\data\Wetterskip.gdb Next step was to perform the Focal Statistics again (through Python), which this time resulted in: Runtime error File <string>, line 1, in <module> File c:\program files\arcgis\desktop10.1\arcpy\arcpy\sa\Functions.py, line 4788, in FocalStatistics File c:\program files\arcgis\desktop10.1\arcpy\arcpy\sa\Utils.py, line 47, in swapper File c:\program files\arcgis\desktop10.1\arcpy\arcpy\sa\Functions.py, line 4782, in wrapper File c:\program files\arcgis\desktop10.1\arcpy\arcpy\geoprocessing\_base.py, line 490, in <lambda> <class 'arcgisscripting.ExecuteError'>: ERROR 000622: Failed to execute (Focal Statistics). Parameters are not valid. ERROR 000628: Cannot set input into parameter in_raster. What in the worlds name was happening with my "in_raster"â?¦ :confused: I took a closer look at my TOC and noticed that for a really strange reason all the datasources had become invalidâ?¦ (see Screenshot_TOC). Fortunately, using the Repair Datasource option fixed all the datasource errors. 🙂 I Think Iâ??m going to restart my system, and see if it happens againâ?¦ Anyone with similar experiences?
... View more
10-26-2011
02:34 AM
|
0
|
3
|
5544
|
|
POST
|
There is a thing that I simply don???t understand or what I am not willing to accept. When using ArcPy.sa to create raster output in the Python window, the name the raster is given in the workspace (in my case a fgdb) is based on the function applied (e.g. ifthe_ras* when using a Con statement, or FocalSt_* when using focal statistics, etc). What I would really like (and I know I don???t stand alone in this) is that the name I have assigned when executing the command is used for the name of the raster in my output workspace. Use the env.overwriteOutput setting and when the output raster exists and should not be overwritten throw an error. But please don???t make my script twice as long since after every command: myRaster = Some ArcPy.sa command I have to say myRaster.save("myRaster").:mad: I???m glad that ArcGIS 10(.1) is a lot more stable than the 9.4 beta. When I was using 9.4, I regularly ended up with a crash and a fgdb filled with a lot of rasters for which it was impossible to derive what they were based on??? Should I submit this in the ideas section, or can we just consider this a bug, that can be resolved in the next release?
... View more
10-26-2011
01:47 AM
|
0
|
2
|
3310
|
| 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 |
2 weeks ago
|