I have a 1x1meter DEM and the original LiDAR point cloud. I need to create a vector feature with roads. Are there any tools in ArcGIS to do this? Or any ideas for the algorithms to try to do this in R or Matlab?
3 Answers
You can use a segmentation process like meanshift in Orfeo Toolbox in QGIS or Monteverdi standalone application. Also, meanshift it's available in ArcGIS.
-
1You have the LiDAR Data? Is the area forested? If so, use con functions to isolate places that have the same elevation in both the DEM and the DTM. Call the other areas NoDATA. This will generate a new surface with a lot of the areas not roads removed.– GBGSep 22, 2017 at 22:57
can you provide a link to the DEM? Considering that roads have an altitude not so different from that of their vicinity, and that they are not lowest points as to apply hydrological functions, i think texture is a feature that may help you to pop out man-made surfaces, in r (rgrass7) you can use grass's r.texture command; you may also consider using other kind of imagery, i used your image (the jpg) for calculating texture (above) and below there's an ortho photo i had where roads are really highlighted, then you can select pixels by digital number and extract the roads
r.import input=C:\Users\Elio\Downloads\b5dhD.jpg output=b5dhD -o
g.region raster=b5dhD.1
r.texture input=b5dhD.1_size5 -a output=text size=5 distance=2 --overwrite
-
1here is a link to the hillshade dropbox.com/s/cuvkth0b9o3f3te/hillshade.7z?dl=0– LizaSep 22, 2017 at 17:54
-
-
1
-
I added the code above; I did it straight in GRASS 7.2.1 command line, in R you have to add execGRASS("r.texture" , input="b5dhD.1", ...) Sep 22, 2017 at 18:49
I doubt very much that full automation of this task is possible at all. Let's test semi-automated procedure suggested in my comments. I placed 2 points at the ends of the longest roads 'visible' on hillshade (note that points are labelled by their FID):
arcpy.gp.Slope_sa("dem", "../SLOPE", "PERCENT_RISE")
arcpy.SelectLayerByAttribute_management("POINTS","NEW_SELECTION", """"FID" = 0""")
arcpy.gp.CostBackLink_sa("POINTS", "SLOPE", "./blink")
arcpy.SelectLayerByAttribute_management("POINTS","SWITCH_SELECTION")
arcpy.gp.CostPath_sa("POINTS", "SLOPE", "blink","./path", "EACH_CELL", "FID")
arcpy.gp.RasterCalculator_sa("""Int(Power(2,"blink"-1))""", "../fdir")
arcpy.gp.StreamToFeature_sa("path", "fdir", "../ROAD.shp", "SIMPLIFY")
Picture below shows newly derived 1800m long road (dashed line) on the top of publicly available centrelines:
Table below is NEAR distance statistics of "ROAD" vertices against existing centrelines:
Not too bad in my opinion and I have reservations about resolution of existing lines.
In terms of efficiency, however it seems too slow compared to methodology described earlier for similar post. It took me under 1 minute to accurately digitise 3500 m of 'visible' roads.