lasgrid half pixel shift with geotiffs

[Steven F]

Hi Martin,
I think lasgrid is causing a half pixel shift in geotiff files when the las header projection info contains the key: "key 1025 tiff_tag_location 0 count 1 value_offset 2 - GTRasterTypeGeoKey: RasterPixelIsPoint"
"RasterPixelIsArea" probably causes the same problem.

When the las file contains this geokey I think it writes an x and y pixel tie point of 0.5 into the geotiff header. The result is that programs like ArcGIS and ENVI load the raster in a position where many of the cells don't actually contain the las points being gridded. I noticed this when comparing the las points against the lasgrid geotiff in ArcMap. So hopefully this isn't just another ArcMap glitch, but I doubt it is. When I rewrote the las header vlr and reran lasgrid the tiff file was properly positioned.


Here is the vlr for the original file:

variable length header record 1 of 1:
  reserved             0
  user ID              'LASF_Projection'
  record ID            34735
  length after header  40
  description          ''
    GeoKeyDirectoryTag version 1.1.0 number of keys 4
      key 1024 tiff_tag_location 0 count 1 value_offset 1 - GTModelTypeGeoKey: ModelTypeProjected
      key 1025 tiff_tag_location 0 count 1 value_offset 2 - GTRasterTypeGeoKey: RasterPixelIsPoint
      key 3072 tiff_tag_location 0 count 1 value_offset 32613 - ProjectedCSTypeGeoKey: PCS_WGS84_UTM_zone_13N
      key 4099 tiff_tag_location 0 count 1 value_offset 9001 - VerticalUnitsGeoKey: Linear_Meter

Then I recode the vlr with: las2las -i original.las -remove_all_vlrs -wgs84 -utm 13N -o new.las

Which changes it to:

variable length header record 1 of 1:
  reserved             43707
  user ID              'LASF_Projection'
  record ID            34735
  length after header  40
  description          'by LAStools of Martin Isenburg'
    GeoKeyDirectoryTag version 1.1.0 number of keys 4
      key 1024 tiff_tag_location 0 count 1 value_offset 1 - GTModelTypeGeoKey: ModelTypeProjected
      key 3072 tiff_tag_location 0 count 1 value_offset 32613 - ProjectedCSTypeGeoKey: PCS_WGS84_UTM_zone_13N
      key 3076 tiff_tag_location 0 count 1 value_offset 9001 - ProjLinearUnitsGeoKey: Linear_Meter
      key 4099 tiff_tag_location 0 count 1 value_offset 9001 - VerticalUnitsGeoKey: Linear_Meter

Then I rerun lasgrid with the same settings: lasgrid -i new.las -step 0.5 -elevation -highest -o new.tif

The new tiff world file is exactly the same as original world file, but I think the offending geotag is no longer in the geotiff header because pixels have now moved to where I expect them. Maybe this behavior is appropriate, but it was unexpectedly causing a small georegistration problem for me. Others might be having the same problem, but are totally unaware of it since the difference is so small.  I just wanted to bring it to your attention .

Cheers,
Steve


Screen shots of las points over original and new tiffs:

[Martin Isenburg]

Hello Stephen,

excellent bug report and wonderfully documented! Thank you! 

And you are right on all counts. I undiscriminatorily added this GeoTIFF tag into the TIFF header whenever it happened to be in the LAS/LAZ header. Not sure why a file filled with points would contain a flag specifying

  key 1025 tiff_tag_location 0 count 1 value_offset 2 - GTRasterTypeGeoKey: RasterPixelIsPoint

Maybe to indicate the LAS/LAZ file was converted from a raster of elevation pixels and the flag informs whether the points represent the corners or the centers of the pixel grid cells? That's the only use case that I can imagine. Either way, in the next release any GTRasterTypeGeoKey in the LAS/LAZ file will be ignored and the TIF files will always specify the GTRasterTypeGeoKey with the correct value RasterPixelIsArea.

Cheers,

Martin

--
Download LAStools at
http://lastools.org
http://rapidlasso.com
Be social with LAStools at
http://facebook.com/LAStools
http://twitter.com/LAStools
http://linkedin.com/groups/LAStools-4408378
Manage your settings at
http://groups.google.com/group/lastools/subscribe

[Steven F]

Thanks Martin,
I'm not sure why my data provider included that tag in the las file, but it could be a default of their processing software which I think is Optech LMS.

Steve

[Evon Silvia]

I've seen the Blue Marble's GeoCalc SDK (which is Optech LMS's datum engine) include that field in their LAS file GeoTIFF encoding, so that's probably the source. If one considers the GTRasterTypeGeoKey field to be required by the GeoTIFF spec, one could argue that the RasterPixelIsPoint field would specify that the dataset is actually a vector (point) dataset.

I personally disagree with that reasoning and simple exclude the GTRasterTypeGeoKey, but I think that's where they're coming from.

Evon

[Beaty, Paul W (Paul Beaty)]

Both the RasterPixelIsArea and RasterPixelIsPoint options should be expected by the application software and handled correctly. Otherwise is limited because it application will always have a problem dealing with geospatial data stored as TIFF or using GeoTIFF tags. The application should expect both options and adjust the data accordingly.

 

From, remotesensing.org, the creater of GeoTIFF tags:

What is the effect of PixelIsArea vs. PixelIsPoint?

Setting the GTRasterTypeGeoKey value to RasterPixelIsPoint or RasterPixelIsArea alters how the raster coordinate space is to be interpreted. This is defined in section 2.5.2.2 of the GeoTIFF specification. In the case of PixelIsArea (default) a pixel is treated as an area and the raster coordinate (0,0) is the top left corner of the top left pixel. PixelIsPoint treats pixels as point samples with empty space between the "pixel" samples. In this case raster (0,0) is the location of the top left raster pixel.

The net effect is a half pixel offset when interpreting PixelIsPoint vs. PixelIsArea images. There has been some confusion on the intepretation of this value in the GeoTIFF community for some time and some held the opinion that this parameter should have no effect on the georeferencing of the image. While this matter has been settled, the result is that some software packages, particularly those built on GDAL 1.7.x or older incorrectly interprete the georeferenced location of PixelIsPoint images - placing them offset by half a pixel. One way of avoiding compatability problems is to produce GeoTIFF files with the default PixelIsArea interpretation.

 

Regards,

 

Paul Beaty

--
Download LAStools at