[gdal-dev] ASCII Gridded XYZ and "Hectare Raster" data (STATPOP and STATENT)

[Even Rouault]

On jeudi 25 avril 2019 00:26:35 CEST Stefan Keller wrote:
> Hi,
> 
> Am 22. Apr. 2019 23:49 wrote regarding gdal_grid:
> > Is there a solution for just setting output pixel size in map units
> 
> By guessing the "-outsize 65 75" (no. of X and Y lines) I finally got
> a result which contained a pixel size of 100, -100 meters plus
> sensible values:
> 
> # gdal_grid -ot UInt16 -of GTiff -zfield "C3" -a
> nearest:radius1=71:radius2=71 -outsize 65 75 in.gpkg out.tif
> 
> => But, I'm still looking for a hint to set pixel sizes without
> guessing or calculating the "outsize" from the extents of the input
> file.

This is left as an exercice to the contributor to implement a -tr option 
similar to gdalwarp / gdal_translate.

> 
> Then, I observed another problem: One obviously has no choice in
> gdal_grid than to choose an interpolation method. As one can see from
> the command above, I've chosen "nearest neighbour" as interpolation
> method (option "-a nearest:radius1=71:radius2=71").

Why 71 ? If you're grid size is 100 meters, then less than 50 meters would be 
a good radius
You would also need to grow the raster extent by a half-pixel size, so that 
the middle of each raster point matches exactly an exact x,y coordinate of 
your CSV (GDAL convention for georeferencing is the coordinate of the top-left 
corner of a pixel area, not its center)
By default, gdal_grid takes the extent of the vector file to set it as the 
extent of the raster file, and doesn't do this half-pixel extension, which is 
the likely cause for the effect you observed.

Example:

given grid.csv with

x,y,z
100,500,1
300,500,3
100,600,4
200,600,5
300,600,6

# note the missing value at 200,500

$ ogr2ogr grid.gpkg grid.csv -oo X_POSSIBLE_NAMES=x -oo Y_POSSIBLE_NAMES=y

# Note the extra small radius since the center of a raster grid matches a
# input vector point, or none
$ gdal_grid -a nearest:radius1=1:radius2=1:nodata=-99  grid.gpkg  grid.tif \
		-txe 50 350 -tye 650 450 -outsize 3 2 -zfield z	

$ gdal_translate grid.tif /vsistdout/ -of aaigrid
ncols        3
nrows        2
xllcorner    50.000000000000
yllcorner    450.000000000000
cellsize     100.000000000000
NODATA_value  -99
 4.0 5 6
 1 -99 3


Even

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com