User Tools

Site Tools


wiki:exercise_createdem

Create a dem from point data

Exercise:
Script: open by kate ~/ost4sem/exercise/basic_adv_gdalogr/dem/dem.sh &
Data: ~/ost4sem/exercise/basic_adv_gdalogr/dem/M0227-[1-4]N.H29
Directory: ~/ost4sem/exercise/basicadvgdalogr/dem

Before running the script migrate to ~/ost4sem/exercise/basicadvgdalogr/dem and return to the following questions:

Do the M0227-[1-4]N.H29 files have a header?
How many lines do the files have?
Is it the same for all the files?
Plot one file in 3D using gnuplot.
Sample the data using awk.
Plot all the files in gnuplot.
Can we create a dem from this file?
Are all the files geographically contiguous.?
Is there an overlap between the files?
Which is the gdal command that allows me to create a surface from points data?


Merge all the files into one file:

  • Get rid of the header lines.
  • Merge the first and second columns.
rm point.asc 
for file in M* ; do
    awk '{if (NF<=3) print $1 $2,$3 }'  $file  >> point.asc  
done

At this point the file point.asc includes all the previous files and can be sorted based on the first column.

sort -k 1,1   point.asc  > point_s.asc 

The point that will have the same coordinates will be closed and can be deleted. Create a header for the csv file.

rm point.csv
echo "\"X\",\"Y\",\"Z\"" > point.csv 

Why there is the “\” simbol?

Create an AWK script able to delete the numbers that are equal. Do not copy/paste the awk script from the wiki, copy it from the source script. THIS CAN BE ALSO DONE WITH A BASH COMMAND search.

awk '{  
    if(NR==1) { 
	    old=$1 
            z=$2 
            print substr($1,1,6) "," substr($1,7,7) "," $2  
	} 
    else 
	{
	    if($1!=old) {
		    print substr(old,1,6) "," substr(old,7,7) "," z ; 
		    old=$1 ; 
		    z=$2}
	}
} 
END { 
    print substr(old,1,6) "," substr(old,7,7) "," z ; 
    old=$1 ; 
    z=$2  
}'  point_s.asc  >> point.csv

The file point.csv can be used as an input file to create a point ogr format.
Read carefully “READING COMMA SEPARATED VALUES” at http://www.gdal.org/gdal_grid.html
and explain how it is working

echo "<OGRVRTDataSource>" > point.vrt
echo "    <OGRVRTLayer name=\"point\">" >>  point.vrt
echo "        <SrcDataSource>point.csv</SrcDataSource>"  >>  point.vrt 
echo "	      <GeometryType>wkbPoint</GeometryType>"  >>  point.vrt  
echo "	      <GeometryField encoding=\"PointFromColumns\" x=\"X\" y=\"Y\" z=\"Z\"/> "  >>  point.vrt 
echo "    </OGRVRTLayer>"  >>   point.vrt  
echo "</OGRVRTDataSource>"  >>   point.vrt 

You have created a point file in ogr format.
You can use gdal_grid to create a regular grid (raster).
Two ways are performed using two interpolation methods.

gdal_grid -a average:radius1=25:radius2=25:angle=0.0:min_points=1:nodata=-9 -outsize 100 100 -of EHdr -ot Float32 -l point point.vrt 25avg.tif 

Open the two 25avg.tif and by openev.
Read the image characteristic by gdalinfo.
Make a new image with pixel size = 100

wiki/exercise_createdem.txt · Last modified: 2017/12/05 22:53 (external edit)