gdal_rasterize can miss points on the extent edges

[EliL]

Running gdal_rasterize -tr 100 100 -a Elevation -l town town.shp townraster.tif converts 12 points to 9 pixels (with values) in the output raster. It should convert it to 11 pixels (with values) in the output raster (1 point has a null value and should not get converted).

town.shp is a point layer with 12 points, 11 of those points have values in -a Elevation (1 is null), and in the output tif, 9 pixels have values. The far south and far east points do not have corresponding pixels with values.

Viewing the results in qgis, the points and pixels do not always line up correctly, perhaps it is the same cause for the extent edges to be missing pixels (i.e. those pixels would be in the next row or column of pixels).

Results should be reproducible with town.shp, attached.

[EliL]

small point shapefile

[mattbk]

I can confirm this behavior with the data provided and with my own. Discussion is occurring ​here. Pixel offset problem can be solved by shifting the raster by 0.5 pixels, but this is not a long-term solution.

[peifer]

Just to give another example: I would assume that the 5 points in test.csv [1] can be rasterised into a 3 x 3 raster with an extents of [xmin ymin xmax ymax] = [0 0 3 3]. This is obviously not true as one can see at [2]. Shifting the extents of the raster in SE direction by 0.5 pixel seems to be necessary to get the expected result, see [3].

[1]

$ cat test.csv
Longitude,Latitude,Name,Value
0.5,0.5,"First point",1
0.5,2.5,"Second point",2
2.5,2.5,"Third point",3
2.5,0.5,"Fourth point",4
1.5,1.5,"Fifth point",5

$ cat test.vrt
<OGRVRTDataSource>
    <OGRVRTLayer name="test">
        <SrcDataSource>test.csv</SrcDataSource>
        <GeometryType>wkbPoint</GeometryType>
        <LayerSRS>EPSG:4326</LayerSRS>
        <GeometryField encoding="PointFromColumns" x="Longitude" y="Latitude"/>
    </OGRVRTLayer>
</OGRVRTDataSource>


----

[2]

$ gdal_rasterize -l test test.vrt test.tif -a Value -tr 1 1 -te 0 0 3 3 &&
  gdal_translate -of aaigrid test.tif test.asc && cat test.asc
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 3, 3
0...10...20...30...40...50...60...70...80...90...100 - done.
ncols        3
nrows        3
xllcorner    0.000000000000
yllcorner    0.000000000000
cellsize     1.000000000000
 0 0 0
 0 2 0
 0 0 5

----

[3]

$ gdal_rasterize -l test test.vrt test.tif -a Value -tr 1 1 -te .5 -.5 3.5 2.5 &&
  gdal_translate -of aaigrid test.tif test.asc && cat test.asc
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 3, 3
0...10...20...30...40...50...60...70...80...90...100 - done.
ncols        3
nrows        3
xllcorner    0.500000000000
yllcorner    -0.500000000000
cellsize     1.000000000000
 2 0 3
 0 5 0
 1 0 4

[Even Rouault]

r22998 /trunk/ (3 files in 3 dirs): gdal_rasterize: fix half pixel shift when rasterizing points; make gdal_rasterize utility increase the computed raster extent by a half-pixel for point layers (#3774)

[warmerdam]

I have backported the llrasterize.cpp portion of this to 1.8 in time for 1.8.2.