gdal_grid question

[Barrett Higman]

Hi all,

 

I have a shapefile with 5.3 million points relating to heights that I want to use gdal_grid to turn into a DEM. The points are derived from a 20m contour interval shapefile so I would like to end up with a 20m x 20m pixel resolution on the DEM. The gdal_grid command I am using is:

 

gdal_grid -ot UInt16 -of GTiff -zfield "ALTITUDE" -l contour_nodes C:\contour_nodes.shp C:\dem_20m --config GDAL_NUM_THREADS ALL_CPUS

 

However, without any input regarding grid output size it automatically creates a DEM with a pixel size of 400m x 400m. I can't tell it to use a particular cell size, only the overall output grid size (using -outsize x x).

 

Can anyone tell me if there is an easy way to determine the appropriate grid size (number of cells and lines) needed to produce the desired cell sizes rather than just undertaking some rough measuring and maths? ie. I would like this to be precise and produce 20m x 20m cells not something that is around about that.

 

Hope this makes sense :)

 

Oh and if anyone had any ideas about how to speed up the process (the 400m x 400m DEM took several hours, the 20m x 20m will take forever!) that would be great too.

 

Cheers,

Barrett

 

 

 

Barrett Higman

GIS Officer

Alpine Shire Council

 

[Shaun Kolomeitz]

Barrett,

Have you had a look at gdal_rasterize ?

Will that do the trick for you ?

Cheers,

Shaun

-- 
You received this message because you are subscribed to the Google Groups "Australian QGIS User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to australian-qgis-use...@googlegroups.com.
To post to this group, send an email to australian-qg...@googlegroups.com.
Visit this group at http://groups.google.com/group/australian-qgis-user-group.
For more options, visit https://groups.google.com/groups/opt_out.
 
 

[Luke]

Perhaps create the grid first in QGIS using research tools and then intersect with your points

Then convert vector grid to DEM

--

Luke Bassett
Melbourne Australia

[Barrett Higman]

Thanks for the reply Shaun, I must admit I have no experience using gdal_rasterize but I don't think it undertakes any type of interpolation? My points are not regular having coming from a contour dataset so I think that gdal command would simply leave gaps everywhere there are no points. Might be wrong though?

 

Cheers,

Barrett

 

Barrett Higman

GIS Officer

Alpine Shire Council

 

[Barrett Higman]

That is a very interesting idea Luke, have you tried it before?

 

The gdal_grid documentation does state that it creates a regular grid from scattered data so I hadn't thought of making the data regular first... I'll give it a go and see what happens.

 

Cheers,

Barrett

 

Barrett Higman

GIS Officer

Alpine Shire Council

 

From: Luke [mailto:coolha...@gmail.com]
Sent: Thursday, 1 August 2013 9:11 AM
To: australian-qg...@googlegroups.com
Subject: Re: gdal_grid question

 

Perhaps create the grid first in QGIS using research tools and then intersect with your points

[Michael Sumner]

gdal_grid has arguments to specify the grid:

[-txe xmin xmax] [-tye ymin ymax] [-outsize xsize ysize]

It would probably help to add the creation options for GeoTIFF -co
TILED=YES -co COMPRESS=LZW, and I would experiment with subsets first
and with the -a option for the interpolation algorithm.

http://www.gdal.org/gdal_grid.html#gdal_grid_algorithms

Also, are you sure that integer output is what you want? I would
experiment there too. In terms of the speed I would test out different
subsamplings first to get a feel for the result, and also see the new
GDAL_NUM_THREADS option for GDAL 1.10.

HTH, I've not really used this utility much, it's a hard thing to get right.

Cheers, Mike.

> --
> You received this message because you are subscribed to the Google Groups
> "Australian QGIS User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to australian-qgis-use...@googlegroups.com.
> To post to this group, send an email to
> australian-qg...@googlegroups.com.
> Visit this group at
> http://groups.google.com/group/australian-qgis-user-group.
> For more options, visit https://groups.google.com/groups/opt_out.
>
>



--

Michael Sumner
Hobart, Australia
e-mail: mdsu...@gmail.com

[Barrett Higman]

First off, thank you for the replies, especially to Shaun for pointing me in the direction of a 64bit version of GDAL: http://www.gisinternals.com/sdk/ and also to Mike for some improvements in output compression, etc.

I have managed to produce a very acceptable output with plenty of mucking around with the available algorithms, however, I was wondering if anyone has any specific knowledge regarding the best use of the available algorithms and their various settings? I have done a lot of trial and error but can't try every combination so it would be good if there was some actual science behind the settings I chose!

Any thoughts would be appreciated.

Cheers,
Barrett



Barrett Higman
GIS Officer
Alpine Shire Council

Great Alpine Road
PO Box 139
Bright Victoria 3741
Telephone: (03) 5755 0552
Mobile: 0439 368 438
Facsimile: (03) 5755 1811
Web: www.alpineshire.vic.gov.au

-----Original Message-----
From: Michael Sumner [mailto:mdsu...@gmail.com]
Sent: Thursday, 1 August 2013 9:32 AM
To: australian-qg...@googlegroups.com
Subject: Re: gdal_grid question

[Michael Sumner]

The best I've seen for contours is the DEST algorithm, but the only
software I have for that is Manifold GIS:

http://onlinelibrary.wiley.com/doi/10.1029/2004JF000150/abstract

It's probably worth encouraging the smart ones to get something like
that working in GDAL, or perhaps it can be done with existing free
tools some how. :)

Actually, there's GRASS - I haven't use this much but it's apparently
easy to install/use from OSGeo4W:
http://grasswiki.osgeo.org/wiki/Contour_lines_to_DEM#r.surf.contour

[Barrett Higman]

I have finally got my 20m DEM but I'll have to admit that I cheated somewhat... after working with gdal_grid (both in QGIS and in commandline) as well as GRASS on a subset of the data I determined that the algorithm settings I wanted were the following:

gdal_grid -of GTiff -co "TILED=YES" -co "COMPRESS=LZW" -zfield "ALTITUDE" -a invdist:power:2.0smoothing:2.0:radius1=300:radius2=300 -l contour_nodes -outsize 4373 5868 C:\contour_nodes.shp C:\dem_20m --config GDAL_NUM_THREADS ALL_CPUS

That is a lot of lines and cells!!

I tried running this on the full dataset but after 6 hours of processing on a machine with quadcore and 8Gb RAM it was still only a quarter of the way through. GRASS simply wouldn't import all the points, it got to about 54% and threw an error.

So (here's the confession), I tried the process out on FME using the RasterDEMGenerator and I had a 20m DEM half an hour later, the majority of that time was spent simply reading in the 5+ million points.

I had considered breaking the point file up into separate, more manageable, subsets, processing those using gdal_grid then stitching it all back together using gdal_merge.py, however, I realised that this would probably result in issues whereby points on the edge of each subset could not be interpolated using the adjoining points on the next subset and therefore introduce errors in each and every subset. If anyone would like to correct my thinking here feel free!

Any comments are welcome, but I just thought I'd provide some closure to this thread :)

Cheers,
Barrett

Barrett Higman
GIS Officer
Alpine Shire Council