29

I'm using gdal_proximity to find the distance to the nearest major river across the USA (lower 48 states). I've projected the NHD+ network flowlines to Conus Albers (epsg:5070), selected rivers with stream order > 5, and rasterized, burning rivers as 255, no river as 0. This much is fine, but now I need to find distance to nearest river for sites within 50km. The input file is at 30m resolution on continental scale, so is very large, but the conversion should be a simple gdal_proximity command:

gdal_proximity.bat -values 255 -distunits GEO -maxdist 50000 -nodata -999 infile.tif outfile.tif -co COMPRESS=DEFLATE -co BIGTIFF=YES -co TILED=YES

This seems to -almost- work, but is producing a weird geometric pattern in the output (see image). The data that are present in the output have processed correctly. Can anyone suggest why so much of the output is missing?

Proximity to river

Edit: To test if this was caused by any of the optional parameters, I ran gdal_proximity again in this configuration:

gdal_proximity.bat H:\data\tmp\NHDplus_network_flowline_SO6plus.tif H:/data/tmp/NHDplus_network_flowline_SO6plus_proximity.tif -values 255 -maxdist 50000 -of GTiff

Which yielded essentially the same result:

Proximity to river, no optional parameters

My only thought is that it may be related to the size of the raster (~100 gb uncompressed. As far as I'm aware there isn't a limit to the size of a BigTiff, but maybe there's a limit to what gdal can analyse effectively?

Milad
  • 28
  • 11
R Rhodes
  • 391
  • 2
  • 5
  • 1
    what happens if you turn off tiled=YES? Also, does it work if you change from GEO to PIXEL? (The output might not be suitable, but it might narrow down the problem) – Steven Kay Jul 04 '17 at 17:09
  • Thanks for the suggestion - have added a reply to the original question. – R Rhodes Jul 07 '17 at 13:04
  • What resolution your infile.tif at? – shahryar Mar 01 '18 at 04:57
  • 2
    Can you try to read the data using GDAL in batches (lines) and see if the problem is the data itself or QGIS not being able to visualize it? A first step to find this problem is reducting the spatial extent to a sample AOI. – RutgerH Apr 24 '18 at 13:45

1 Answers1

3

I suspect that you are hitting a memory limit somewhere, possibly when RAM is exhausted and the OS dumps to a paging file. Monitor your system resources during the process. It is unclear to me why your results occur in curved swaths, but do ensure you have projected (saved) all the data into the same coordinate system.

Lets take a look at numeric datatypes in order to help this algorithm. The rasterized stream network need only contain binary values, so we can save on resources by using a Byte raster data type. Burn a value of 1 for streams and 0 for background:

gdal_rasterize -l streams -burn 1 -tr 50 50 -a_nodata 0 -te -2339101 311625 2227004 3134200 -ot Byte -of GTiff streams.shp streams.tif

Next, the proximity we are interested in is positive and less than or equal to 50,000m. An appropriate datatype is an unsigned 16-bit integer UInt16. Furthermore, if we set the 'no data' to the max 65535 we can retain a 0 value for the stream cells.

If necessary, you could also step down to an 8-bit unsigned integer UInt8 and still have ~200m proximity precision.

gdal_proximity.bat -srcband 1 -distunits GEO -values 1 -maxdist 50000 -nodata 65535 -ot UInt16 -of GTiff streams.tif proximity.tif

*Note that I've used a 50m cell size. The gdal_proximity consumed ~20GB RAM and took ~5min on my machine. If you are RAM limited, split the input raster into manageable sizes, as others noted.

gdal_proximity results

CrystallineEntity
  • 510
  • 6
  • 13