I have a 10m USGS DEM for my entire map extent (~100sq miles)
Inside of that, I have specific regions of higher quality GPS point data. (From a rover using a OPUS corrected RTK base station)
How do I seamlessly merge my data/data products with the "background" USGS data? If I simply merge/replace the USGS data with mine, the edges don't align and my entire region looks "sunk in" or "raised up" in comparison to the USGS "background"
The gaps that are being filled in by the USGS DEM are up to a full square mile, so I can't just interpolate the blank area from the edges of my survey.
My first though is to "feather" the USGS DEM. Obviously at the edge where the two meet, my data is more accurate, so feathering the USGS data to meet smoothly with my survey seems like it should work, but I am unsure of the technical way to approach this or the accuracy of this method.
The data products that are being produced from this are hillshades, viewsheds, streamflow, contours, and more - and none of them seem to work well with a sharp boundary transition.
Answer
- Create outline of new dem
- convert to polyline Create transition zone buffer around it. to polyline
- Interpolate shape using old dem
- Interpolate shape using new dem
- Triangulate using both. Map calculator will do the rest
UPDATE AND CORRECTIONS:
to simulate distortions I took this DEM:
clipped and flipped part in the middle and merged original with flipped:
Workflow:
arcpy.gp.RasterCalculator_sa("""Con( ~IsNull("flipped"),1)""","C:/FELIX_DATA/SCRARCH/one")
arcpy.RasterToPolygon_conversion("one","C:/FELIX_DATA/SCRARCH/outline.shp")
arcpy.FeatureToLine_management("outline","C:/FELIX_DATA/SCRARCH/boundary.shp")
arcpy.Buffer_analysis("boundary","C:/FELIX_DATA/SCRARCH/buffer.shp","20 Meters","FULL")
arcpy.FeatureToLine_management("buffer","C:/FELIX_DATA/SCRARCH/to_interpolate.shp","#","ATTRIBUTES")
arcpy.gp.RasterCalculator_sa("""Con(IsNull("flipped"),"dem")""","C:/FELIX_DATA/SCRARCH/w_hole")
arcpy.InterpolateShape_3d("w_hole","to_interpolate","C:/FELIX_DATA/SCRARCH/outside.shp")
arcpy.InterpolateShape_3d("flipped","to_interpolate","C:/FELIX_DATA/SCRARCH/inside.shp","#","1","BILINEAR","DENSIFY","0")
arcpy.CreateTin_3d("C:/FELIX_DATA/SCRARCH/tin","#","inside Shape.Z Hard_Line ;outside Shape.Z Hard_Line ;buffer Soft_Clip ","DELAUNAY")
arcpy.TinRaster_3d("TIN","C:/FELIX_DATA/SCRARCH/transition","FLOAT","LINEAR","CELLSIZE 1","1")
arcpy.MosaicToNewRaster_management("transition;flipped;dem","C:/FELIX_DATA/SCRARCH","smooth","#","32_BIT_FLOAT","#","1","FIRST","FIRST")
OUTPUT:
No comments:
Post a Comment