POST
|
What I meant by using ArcMap tools, was not actually using the raster calculator, but the spatial analys "Con" tool. What I would do is open ArcMap, open the catalog window, right-click on Toolbox (or My Toolboxes) folder and choose "New" Model. Perfect, this worked!!! I separated the process into two different models. I performed a Minus first, and used that as the input raster which had the conditional evaluated on. Thanks so much for your help and patience!
... View more
07-01-2013
05:14 PM
|
0
|
0
|
321
|
POST
|
This works: arcpy.gp.RasterCalculator_sa("""elevation = Con("Surface_RAS"-"Ground_RAS">=0,"Surface_RAS"-"Ground_RAS",0)""","C:/Summer_13/Geo.gdb/elevation") But trying to edit it so that it takes the variables doesn't work. DEM = arcpy.Raster(raster_list[0])
DSM = arcpy.Raster(raster_list[1])
expression = """elevation = Con({}-{}>=0,{}-{},0)""".format(DSM,DEM,DSM,DEM)
raster_out = os.path.join(FD_addr,"elevation")
arcpy.gp.RasterCalculator_sa(expression,raster_out) Probably something is wrong with the expression variable. I have to look into figuring out how to convert the """....""" to make my expression
... View more
07-01-2013
03:56 PM
|
0
|
0
|
321
|
POST
|
Yes, I can create the elevation raster from the two rasters produced by the code by using the Raster Calculator. I tried copying the code as a python snippet, but I read on the documentation for the Raster Calculator that it is only supported inside arcmap and will not work in a python script. The equivalent is supposedly how I'm trying to implement it. In Raster Calculator, I type in the following: elevation = Con("Surface_RAS"-"Ground_RAS">=0,"Surface_RAS"-"Ground_RAS",0) The python snippet is: arcpy.gp.RasterCalculator_sa("""elevation = Con("Surface_RAS"-"Ground_RAS">=0,"Surface_RAS"-"Ground_RAS",0)""","C:/Summer_13/Geo.gdb/elevation") (Surface_RAS and Ground_Ras are the two raster inputs that the code creates.) Works perfectly. I'm creating a new Python toolbox, and running the script as a tool. DEM format? I'm not sure I know what you're referencing. Not really sure how to see Arcmessages as I do everything in the command line, and not really familiar with the Arc command window/python tools. How are you running this? In the command window, or making a script tool or what? If running in console window, after it runs, if you type:
print elevation
print raster0
print raster1
do you get any results? also, can you create the elevation raster that you are after within arcmap using the tools? I am wondering if it is creating an empty output, so is not saving it. I would make sure that it is possible to do what you are trying using arcmap tools/model and get the desired output. then work on the python code (export from the model builder and modify). Might also look into the DEM format. I have tried to use this in the past for certain tasks, and found that it does not support state plane coordinates. Perhaps one of the input rasters is not supported or valid for this operation? If you can do the raster calculations in ArcMap, then that would rule this out as well. R_
... View more
07-01-2013
03:23 PM
|
0
|
0
|
321
|
POST
|
So this is what the part of the code that isn't creating the raster looks like (libraries and extensions included)
import arcpy
import os
from arcpy import env
from arcpy.sa import *
arcpy.CheckOutExtension("Spatial") # for arcpy.sa
env.workspace = FD_addr # DEM/DSM raster locations
raster_list = arcpy.ListRasters("*", "ALL") # adds DEM/DSM raster
size = len(raster_list) #### Make some exception here if != 2
raster0 = arcpy.Raster(raster_list[0])
raster1 = arcpy.Raster(raster_list[1])
elevation = arcpy.sa.Con(raster1-raster0 >= 0, raster1-raster0,0)
elevation.save("elevation") # saved in FD_addr
... View more
07-01-2013
02:04 PM
|
0
|
0
|
321
|
POST
|
It doesn't report anything. I'm not sure how I see command line messages (I'm very new to this.) The only messages I can see are the ones I code in by arcpy.AddMessage("Message") So the first option didn't give me any message (though I tried it with the arcpy.AddMessage() as well) and I tried elevation.save() I have since changed it to elevation.save("elevation") because I figured it needed the name to save the raster to since the workspace was already set and knew where to save. I assume I have writing privileges because before this I have a terrain to raster method, and the terrains are being saved to the same location as I'm trying to save this raster (FD_addr), but this is before I set the workspace. I only set the workspace because I wanted a list created with the rasters that are saved in the FD_addr location. So unless setting the workspace alters the writing privilege, I assume I can write files to that location. I'm still at a loss. I feel like I've tried so many options and still don't know what's wrong or how to fix this. It seems as if it should work! Interesting. Are you sure you hard-coded it to an "existing" FGDB? If so, what do you mean it didn't work? Didn't create the raster with your hard-coded name, created a blank raster? Error messages? What does it report if make the modifications in red?
import arcpy, os
arcpy.CheckOutExtension('Spatial')
arcpy.AddMessage("Creating Elevation Raster")
arcpy.env.workspace = FD_addr
list = arcpy.ListRasters("*", "ALL")
raster0 = arcpy.Raster(list[0])
raster1 = arcpy.Raster(list[1])
elevation = arcpy.sa.Con(raster1 - raster0 >= 0, raster1 - rasters0, 0)
print "temporary: ",elevation.isTemporary
elevation.save(FD_addr)
print "temporary: ",elevation.isTemporary
R_ Also, it looks like elevation.save(FD_addr) is telling it to save with the filename of the workspace. Workspace is normally a folder or FGDB, not a file. It's my understanding that elevation.save() would save it in the current workspace as raster dataset "elevation". Now that I think about it, before trying the above code, what happens if you try this (of course, make sure the FD_addr is set to a FGDB that the current user has write access to):
import arcpy, os
arcpy.CheckOutExtension('Spatial')
arcpy.AddMessage("Creating Elevation Raster")
arcpy.env.workspace = FD_addr
list = arcpy.ListRasters("*", "ALL")
raster0 = arcpy.Raster(list[0])
raster1 = arcpy.Raster(list[1])
elevation = arcpy.sa.Con(raster1 - raster0 >= 0, raster1 - rasters0, 0)
elevation.save()
... View more
07-01-2013
02:01 PM
|
0
|
0
|
759
|
POST
|
I don't see in the code anywhere that you are setting the variable 'FD_addr' (maybe it is above the imports somewhere?) to actually set the env.workspace. Could try somehting like:
elevation.save("C:/output/file_gdb.gdb/elevation") # or whatever the path/filename you want for your out raster
If this works, at least then you know it's an issue with the workspace/path being set properly. R_ Still didn't work 😕 I'm setting the workspace here: env.workspace = FD_addr list = arcpy.ListRasters("*", "ALL") (not sure how you make it so that code appears in gray?) I feel that has to be working because the rasters being inserted into list are the rasters in that location (I double checked.) FD_addr is set way in the beginning of the code.
... View more
06-27-2013
11:58 AM
|
0
|
0
|
759
|
POST
|
Yep - on the last line try: elevation.save("my_output") The raster "my_output" should then be saved to the 'FD_addr' workspace. Question though: Does 'FD_addr' point to a feature dataset within a geodatabase? Last I checked, I thought that you can't write rasters to a feature dataset (has to be directly to the GDB). Tried both elevation.save("my_output") my_output.save(FD_addr) "my_output".save(FD_addr) But again, neither worked/saved. FD_addr is really just bad naming on my behalf. This is part of something much larger, and when I first began this, I had intended for FD_addr to be the location inside the (file) geodatabase where the feature dataset would be created.. So FD_addr points to within the geodatabase and not within the feature dataset. Bad naming, I know... But I have other rasters being created with FD_addr as their output address and they do indeed go to the correct spot.
... View more
06-27-2013
09:12 AM
|
0
|
0
|
759
|
POST
|
The inputs need to be "raster objects" to work in teh new map algebra. Try this: import arcpy, os
arcpy.CheckOutExtension('Spatial')
arcpy.AddMessage("Creating Elevation Raster")
arcpy.env.workspace = FD_addr
list = arcpy.ListRasters("*", "ALL")
raster0 = arcpy.Raster(list[0])
raster1 = arcpy.Raster(list[1])
elevation = arcpy.sa.Con(raster1 - raster0 >= 0, raster1 - rasters0, 0)
elevation.save(FD_addr) While this went through without giving any error message, no raster was created (or saved?).
... View more
06-27-2013
08:37 AM
|
0
|
0
|
759
|
POST
|
I cannot get a conditional evaluation calculation to be done on two rasters. I want an ouput raster to be created from the positive difference of two rasters (so no negative numbers) The part of the code that isn't working import arcpy import os from arcpy import env from arcpy.sa import * import arcgisscripting arcpy.CheckOutExtension('Spatial') I've tried two different ways, both failing Method one arcpy.AddMessage("Creating Elevation Raster") env.workspace = FD_addr list = arcpy.ListRasters("*", "ALL") elevation = Con(list[1]-list[0]>=0, list[1]-list[0],0) elevation.save(FD_addr) Method two arcpy.AddMessage("Creating Elevation Raster") expression = "Con({}-{}>=0, {}-{},0)".format(list[1], list[0], list[1], list[0]) OR = "elevation3" OutRaster = os.join.path(FD_addr,OR) gp = arcgisscripting.create() gp.SingleOutputMapAlgebra_sa(expression, OutRaster) my list has two rasters (I've checked to make sure they are there): a DSM and a DEM, and I want to find the positive elevation difference between the two. Any help would be appreciated! I've been stuck on this for ages.
... View more
06-27-2013
07:56 AM
|
0
|
15
|
8172
|
POST
|
In the terrain to raster tool, you have two options for the sampling distance. ArcMap seems to automatically choose the cell_size number when you select CELLSIZE as the method. What math/logic goes behind choosing this number? How would I go about determining for myself what the number should be instead of relying on the ArcMap supplied number? Thanks!
... View more
06-25-2013
07:42 AM
|
0
|
0
|
548
|
POST
|
Hi, I am trying to create a terrain from an input .las file going from .las to multipoint and then the multipoint feature class to terrain. I can successfully create .las to multipoint, but cannot go from the multipoint to terrain step. I am copying the entire code so it is easier to follow and will bold the section with a problem I can get up to the point where I need to add the feature class(es) to the terrain, but at this point, the code stops. I do not get an error warning, it just states code completed without finishing. I know the error is in how I'm trying to add the feature class, but I can't seem to get it to work. import arcpy import os # ArcGIS Parameter Options Geo_output = arcpy.GetParameterAsText(0) # Input, Folder Geo_name = arcpy.GetParameterAsText(1) # Input, String SR = arcpy.GetParameterAsText(2) # Input, Spatial Reference InLAS = arcpy.GetParameterAsText(3) # Input, File, Multi Value, Filter = las pnt_spacing = arcpy.GetParameterAsText(4) # Input, Double OutFC_1 = arcpy.GetParameterAsText(5) # Input, String class_1 = arcpy.GetParameterAsText(6) # Input, Long, Optional return_1 = arcpy.GetParameterAsText(7) # Input, String, Optional, Default = ANY_RETURNS OutFC_2 = arcpy.GetParameterAsText(8) # Input, String class_2 = arcpy.GetParameterAsText(9) # Input, Long, Optional return_2 = arcpy.GetParameterAsText(10) # Input, String, Optional, Default = ANY_RETURNS clip = arcpy.GetParameterAsText(11) # Input, Feature Class, Optional OutFC = [OutFC_1, OutFC_2] classif = [class_1, class_2] returns = [return_1, return_2] FD_name = Geo_name + "_Feature_Dataset" FD_addr = os.path.join(Geo_output,Geo_name + ".gdb") FC_addr = os.path.join(FD_addr, FD_name) try: arcpy.CreateFileGDB_management(Geo_output, Geo_name) arcpy.CreateFeatureDataset_management(FD_addr, FD_name, SR) for fc,cls,ret in zip(OutFC,classif,returns): result= arcpy.LASToMultipoint_3d(InLAS,os.path.join(FC_addr,fc),pnt_spacing, cls,ret,"#",SR,"las","1","NO_RECURSION") name = fc + "_Terrain" arcpy.CreateTerrain_3d(FC_addr,name,pnt_spacing,"50000","#", "ZTOLERANCE","ZMIN","NONE","1") terrain = os.path.join(FC_addr,name) arcpy.AddTerrainPyramidLevel_3d(terrain,"ZTOLERANCE",".5 1000") embed = fc + "_embed" feature_class = result + ' Shape Mass_Points 1 0 0.5 true false' + embed+ ' <None> false' arcpy.AddFeatureClassToTerrain_3d(terrain, feature_class) # Insert Clip arcpy.BuildTerrain_3d(MP,"NO_UPDATE_EXTENT") except: print arcpy.GetMessages(2)
... View more
06-17-2013
02:31 PM
|
0
|
0
|
482
|
Online Status |
Offline
|
Date Last Visited |
11-11-2020
02:24 AM
|