If you have to to process data delivered in plain ASCII files you have probably run into concave hull problem. Contrary to convex hull problem GIS tools usually don’t have some builtin function that can help you with this issue. There is a good demonstration of PostGIS ST_ConcaveHull function
This is just to illustrate multiple possible solution to concave hull problem.
I this post presented scenario is not generic and it works best for regularly spaced points or almost regularly spaced points. With some adjustments can be applied for irregularly scattered points. For now I’ll describe solution for regularly spaced points. Although this procedure is not general I hope it will eventually help somebody.
The idea is to create raster/grid file from point data and to identify areas (pixels) that doesn’t contain any point. The easiest way is to produce grid in the way that value of each pixel represents number of points. For regularly spaced data if everything is setup well we should get only values 0 and 1.
From the created grid just generate contour with value 1. It represents concave hull for the regularly spaced data. If you have almost regularly spaced data you can apply the same principle and just generate contour with value 0.1 (this also works for regularly spacing) . To illustrate difference here is the figure with three contours with values 0.1, 0.5 and 1.
I’ll will start with simple ASCII file with data organized into three columns: x, y, z. That is very common data format that many professional use to share point data. Here is just few lines from used file dtm_predicted_tides.csv (in this scenario file does not contain headers).
446034.00, 6707606.00, 113.92 445036.00, 6707606.00, 113.95 445038.00, 6707606.00, 113.98 445040.00, 6707606.00, 114.00 ...First we have to create VRT (virtual format) file because GDAL cannot directly read and “understand” ASCII file. VRT file details you can find here.
For this application (no headers in ASCII file) VRT file will be (dem.vrt):
<OGRVRTDataSource>
<OGRVRTLayer name="dtm_predicted_tides">
<SrcDataSource>dtm_predicted_tides.csv</SrcDataSource>
<GeometryType>wkbPoint</GeometryType>
<GeometryField encoding="PointFromColumns" x="field_1" y="field_2" z="field_3"/>
</OGRVRTLayer>
</OGRVRTDataSource>
To test it use command ogrinfo:
M:\__gdal_test>ogrinfo -al -so dem.vrt
INFO: Open of `dem.vrt'
using driver `VRT' successful.
Layer name: dtm_predicted_tides
Geometry: Point
Feature Count: 433752
Extent: (444448.000000, 6707606.000000) - (447314.000000, 6709738.000000)
Layer SRS WKT:
(unknown)
field_1: String (0.0)
field_2: String (0.0)
field_3: String (0.0)
From this info data you will need extent to calculate some parameters that you will need in next step. To assure that generated contour line is closed it is necessary to extend area for at least one pixel value. The size of pixel is determined by resolution of original data. For regularly spaced points this value is equal to distance of two neighboring points along the easting and northing axes. For almost regularly spaced data you have to estimate this value bearing in mind the proposed scenario for concave hull estimation.
Taking into account generated extent and resolution of 2m it can be easily calculated that the generated grid should be sized 1433 x 1066. To extend area and to assure that point is located in the middle of generated pixel it is required to extend area for one and half pixel and in this case that is 3m.
Now we have all values necessary to run command. Notice that the size of grid is also changed because original extent was changed.
gdal_grid -a count:radius1=1:radius2=1
-txe 444445 447317 -tye 6707603 6709741
-outsize 1436 1069
-of GTiff -ot Float32
-l dtm_predicted_tides dem.vrt concavehull_ext2.tiff
Usage: gdal_grid [--help-general] [--formats]
[-ot {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/
CInt16/CInt32/CFloat32/CFloat64}]
[-of format] [-co "NAME=VALUE"]
[-zfield field_name]
[-a_srs srs_def] [-spat xmin ymin xmax ymax]
[-clipsrc |WKT|datasource|spat_extent]
[-clipsrcsql sql_statement] [-clipsrclayer layer]
[-clipsrcwhere expression]
[-l layername]* [-where expression] [-sql select_statement]
[-txe xmin xmax] [-tye ymin ymax] [-outsize xsize ysize]
[-a algorithm[:parameter1=value1]*] [-q]
<src_datasource> <dst_filename>
Available algorithms and parameters with their's defaults:
Inverse distance to a power (default)
invdist:power=2.0:smoothing=0.0:radius1=0.0:radius2=0.0:angle=0.0:max_points=0:min_points=0:nodata=0.0
Moving average
average:radius1=0.0:radius2=0.0:angle=0.0:min_points=0:nodata=0.0
Nearest neighbor
nearest:radius1=0.0:radius2=0.0:angle=0.0:nodata=0.0
Various data metrics
:radius1=0.0:radius2=0.0:angle=0.0:min_points=0:nodata=0.0
possible metrics are:
minimum
maximum
range
count
average_distance
average_distance_pts
Notice that we have used count metrics with the boath radius of search ellipse set to 1. The intention was to get one point per pixel. With data that are not regularly spaced this is not possible to get so in that case you should just play a little bit firstly with pixel size and then by search ellipse radius. At the end just interpolate... lets say 0.1 contour and you can get satisfactory results. .... Use gdalinfo to check generated raster
M:\__gdal_test>gdalinfo concavehull_ext2.tiff
Driver: GTiff/GeoTIFF
Files: concavehull_ext2.tiff
concavehull_ext2.tiff.aux.xml
Size is 1436, 1069
Coordinate System is:
PROJCS["ED_1950_UTM_Zone_31N",
GEOGCS["GCS_European_1950",
DATUM["European_Datum_1950",
SPHEROID["International_1924",6378388.0,297.0]],
PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]],
PROJECTION["Transverse_Mercator"],
PARAMETER["False_Easting",500000.0],
PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",3.0],
PARAMETER["Scale_Factor",0.9996],
PARAMETER["Latitude_Of_Origin",0.0],
UNIT["Meter",1.0],
AUTHORITY["EPSG","23031"]]
Origin = (444445.000000000000000,6707603.000000000000000)
Pixel Size = (2.000000000000000,2.000000000000000)
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( 444445.000, 6707603.000) ( 1d59'19.46"E, 60d29'57.45"N)
Lower Left ( 444445.000, 6709741.000) ( 1d59'17.31"E, 60d31' 6.55"N)
Upper Right ( 447317.000, 6707603.000) ( 2d 2'27.63"E, 60d29'58.84"N)
Lower Right ( 447317.000, 6709741.000) ( 2d 2'25.59"E, 60d31' 7.94"N)
Center ( 445881.000, 6708672.000) ( 2d 0'52.50"E, 60d30'32.70"N)
Band 1 Block=1436x1 Type=Float32, ColorInterp=Gray
Min=0.000 Max=1.000
Minimum=0.000, Maximum=1.000, Mean=0.283, StdDev=0.450
Metadata:
STATISTICS_MAXIMUM=1
STATISTICS_MEAN=0,28255913031469
STATISTICS_MINIMUM=0
STATISTICS_STDDEV=0,45024393416031
Pixel size and corner coordinates should be as it was previously calculated. You won’t get info about statistics and coordinate system but for now that is not important.
In this moment we are ready to create concave hull or to be more precise we are able to generate contours. The following command will generate shape file with 0.1 contour and attribute elev that holds contour value.
gdal_contour -a elev concavehull_ext2.tiff concavehull_ext2_z.shp -fl 0.1
Usage: gdal_contour [-bHere is generated concave hull for my case. Gray pixels contain no points and for the other pixels count value is equal to 1 (that was intention from the beginning).] [-a ] [-3d] [-inodata] [-snodata n] [-f ] [-i ] [-off ] [-fl ...] [-nln ] [-q]
gdal_contour -a elev concavehull_ext2.tiff contours.shp -fl 0.1 0.5 0.7If you need to visualise points you can use ogr2ogr to create shape file from vrt.
ogr2ogr dem.shp dem.vrtHere is segment with metric raster, points and concave hull (contour=1).
<LayerSRS>EPSG:23031</LayerSRS>The previous is sufficient to propagate coordinate system definition to all generated files.
The whole process is very simple and can be summarized into few lines:
1. create vrt file and check it (ogrinfo); optionaly crate shape file ogr2ogr
2. create metric raster gdal_grid (check it with gdalinfo)
3. generate concave hull (gdal_contour)
... just play with your data and enjoy! :-)
No comments:
Post a Comment