[pdal] pdal tindex --merge problem with polygon clippping
Howard Butler
howard at hobu.co
Fri Aug 14 05:40:23 PDT 2015
Stefan,
Indeed there was some hackery in the tindex command before. I was using it for EPSG:4326 only, and I had made some changes to support just that.
I think tindex should be updated support what you want to do, as this is the common workflow pattern that gdaltindex uses which we are aiming to copy. One caveat with the merge scenario is that it will not be memory efficient for very large merges. You might need to do multi-stage operations if you are reading many billions of points to then merge. Currently, only the bounds of the WKT are used as a spatial filter (and passed down to OGR) to limit the amount of data to be merged. We should add another option to add typical OGR layer attribute filters as well.
Can you file a ticket with this stuff and we'll take it on?
Thanks,
Howard
> On Aug 14, 2015, at 6:01 AM, Stefan Ziegler <stefan.ziegler.de at gmail.com> wrote:
>
> Hi
>
> I created a small test case here:
>
> http://www.catais.org/tmp/pdal_tindex.zip
>
> It contains a small las file and the tileindex file. This is what I've done:
>
> 1) Create tileindex:
>
> find . -iname "*.las" | pdal tindex tileindex.shp --a_srs EPSG:21781 --t_srs EPSG:21781 --fast_boundary --stdin
>
> 2) Commented out line 161 to allow "--t_srs". Pdal now does *not* reproject data anymore to EPSG:4326. But when trying to write the las file it stops with the following error:
>
> PDAL: Unable to convert scaled value (2.04e+10) to int32 for dimension 'X' when writing LAS/LAZ file tmp.las.
>
> 3) Commented out writerOptions.add() methods for scales and offsets on line 445 - 450 and added setCommonOptions(writerOptions).
>
> Now I get the desired result by executing:
>
> pdal tindex --scale "1 1 1" --merge tileindex.shp --lyr_name tileindex --a_srs EPSG:21781 --t_srs EPSG:21781 tmp.las
>
> So, this works for me but definitely introduces some evil stuff :) See also attached patch to exactly see what I've done.
>
> regards
> Stefan
>
>
> On Fri, Aug 14, 2015 at 6:16 AM, Stefan Ziegler <stefan.ziegler.de at gmail.com> wrote:
> find . -iname "*.laz" | pdal tindex tileindex.shp --a_srs EPSG:21781 --t_srs EPSG:21781 --fast_boundary --stdin
>
> ogrinfo -al tileindex.shp:
>
> Layer name: tileindex
> Geometry: Polygon
> Feature Count: 6
> Extent: (592000.000000, 226000.000000) - (594000.000000, 230000.000000)
> Layer SRS WKT:
> PROJCS["CH1903_LV03",
> GEOGCS["GCS_CH1903",
> DATUM["CH1903",
> SPHEROID["Bessel_1841",6377397.155,299.1528128]],
> PRIMEM["Greenwich",0],
> UNIT["Degree",0.017453292519943295]],
> PROJECTION["Hotine_Oblique_Mercator_Azimuth_Center"],
> PARAMETER["latitude_of_center",46.95240555555556],
> PARAMETER["longitude_of_center",7.439583333333333],
> PARAMETER["azimuth",90],
> PARAMETER["scale_factor",1],
> PARAMETER["false_easting",600000],
> PARAMETER["false_northing",200000],
> PARAMETER["rectified_grid_angle",90],
> UNIT["Meter",1]]
> location: String (254.0)
> srs: String (254.0)
> modified: Date (10.0)
> created: Date (10.0)
> OGRFeature(tileindex):0
> location (String) = /home/stefan/tmp/lidar/srs/./LAS_593227.laz
> srs (String) = EPSG:21781
> modified (Date) = 2015/08/13
> created (Date) = 2015/08/13
> POLYGON ((593000 227000,593000 228000,594000 228000,594000 227000,593000 227000))
>
> pdal info LAS_593227.laz --metadata:
>
> {
> "filename": "LAS_593227.laz",
> "metadata":
> {
> "comp_spatialreference": "PROJCS[\"CH1903 / LV03\",GEOGCS[\"CH1903\",DATUM[\"CH1903\",SPHEROID[\"Bessel 1841\",6377397.155,299.1528128,AUTHORITY[\"EPSG\",\"7004\"]],TOWGS84[674.4,15.1,405.3,0,0,0,0],AUTHORITY[\"EPSG\",\"6149\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4149\"]],PROJECTION[\"Hotine_Oblique_Mercator_Azimuth_Center\"],PARAMETER[\"latitude_of_center\",46.95240555555556],PARAMETER[\"longitude_of_center\",7.439583333333333],PARAMETER[\"azimuth\",90],PARAMETER[\"rectified_grid_angle\",90],PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",600000],PARAMETER[\"false_northing\",200000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Y\",EAST],AXIS[\"X\",NORTH],AUTHORITY[\"EPSG\",\"21781\"]]",
> "compressed": true,
> "count": 14724002,
> "creation_doy": 224,
> "creation_year": 2015,
> "dataformat_id": 3,
> "dataoffset": 2181,
> "filesource_id": 0,
> "global_encoding": "AAA=",
> "header_size": 227,
> "major_version": 1,
> "maxx": 593999.98999999999,
> "maxy": 227999.98999999999,
> "maxz": 1667.6600000000001,
> "minor_version": 2,
> "minx": 593000,
> "miny": 227000,
> "minz": 622.80000000000007,
> "offset_x": 0,
> "offset_y": 0,
> "offset_z": 0,
> "project_id": "00000000-0000-0000-0000-000000000000",
> "scale_x": 0.01,
> "scale_y": 0.01,
> "scale_z": 0.01,
> "software_id": "PDAL 1.0.0.b1 (e412bd)",
> "spatialreference": "PROJCS[\"CH1903 / LV03\",GEOGCS[\"CH1903\",DATUM[\"CH1903\",SPHEROID[\"Bessel 1841\",6377397.155,299.1528128,AUTHORITY[\"EPSG\",\"7004\"]],TOWGS84[674.4,15.1,405.3,0,0,0,0],AUTHORITY[\"EPSG\",\"6149\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4149\"]],PROJECTION[\"Hotine_Oblique_Mercator_Azimuth_Center\"],PARAMETER[\"latitude_of_center\",46.95240555555556],PARAMETER[\"longitude_of_center\",7.439583333333333],PARAMETER[\"azimuth\",90],PARAMETER[\"rectified_grid_angle\",90],PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",600000],PARAMETER[\"false_northing\",200000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Y\",EAST],AXIS[\"X\",NORTH],AUTHORITY[\"EPSG\",\"21781\"]]",
> "system_id": "PDAL"
> },
> "pdal_version": "PDAL 1.0.0.b1 (e412bd)"
> }
>
> In the mergeFile() method on line 413 you assign the output srs for the reprojection stage. But this cannot be set by the user and will therefore always be EPSG:4326. As far as I understand.
>
> regards
> Stefan
>
>
>
>
> On Fri, Aug 14, 2015 at 12:10 AM, Andrew Bell <andrew.bell.ia at gmail.com> wrote:
> On Thu, Aug 13, 2015 at 2:46 PM, Stefan Ziegler <stefan.ziegler.de at gmail.com> wrote:
> Hi
>
> I'm trying to merge some las files with the following command:
>
> pdal tindex --merge /home/stefan/tmp/lidar/srs/tileindex.gpkg --lyr_name tileindex tmp.las
>
> The input las files have an srs assigned (EPSG:21781). But the merged one is EPSG:4326.
>
> Can you use ogrinfo to verify the SRS of the layer? Can you provide the command used to create the index file?
>
> Thanks,
>
> --
> Andrew Bell
> andrew.bell.ia at gmail.com
>
>
> <TIndexKernel.patch>_______________________________________________
> pdal mailing list
> pdal at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/pdal
More information about the pdal
mailing list