[pdal] Crop SRS not working as expected

[Andrew Bell]

Hi Andreas,

If the SRS of the bounding region (you specified 4326) doesn't match that
of the input, it will be transformed to match that of the points before
cropping occurs, which probably means that what you specified as a box
becomes something else (rhombus?). If you want your crop region to be a
square in the projection of your points, you'll have to leave off the
"a_srs" and use the proper coordinates for the projection of the points.

You could also transform the points to 4326 *before* applying the crop
filter, in which case you can drop "a_srs" from the crop filter as well,
and then you'll get output in 4326.

Hope that helps,

On Tue, Aug 17, 2021 at 11:44 AM Andreas Yankopolus <andreas at yank.to> wrote:

> All—
>
> I’m cropping and processing USGS LPC data into WGS84 raster files, but I
> can’t get the cropping to work out correctly: My pdal crop stage specifies
> a square region in WGS8, but the resulting files display as rhombuses.
> Here’s a minimal example:
>
> {
>   "pipeline": [
>     "/home/ayank/Documents/geodata/lpc/*.laz",
>     {
>       "type": "filters.crop",
>       "bounds": "([-117.18,-117.13],[32.68,32.73])",
>       "a_srs": "EPSG:4326"
>     },
>     {
>       "type": "filters.range",
>       "limits": "Classification[2:2]"
>     },
>     {
>       "filename": "/home/ayank/Documents/geodata/export/sd_coord_test.laz",
>       "type": "writers.las",
>       "compression": "laszip",
>       "extra_dims": "all"
>     },
>     {
>
>   "filename": "/home/ayank/Documents/geodata/export/sd_coord_test.tiff",
>       "type": "writers.gdal",
>       "gdaldriver": "GTiff",
>       "dimension": "Z",
>       "output_type": "idw",
>       "resolution": "10"
>     }
>   ]
> }
>
> My expectation is that the output files, when read into QGIS with WGS84
> (EPSG:4326) as the project CRS, would render as squares. But they’re both a
> rhombus.
>
> I’ve also tried adding these two stages to the pipeline:
>
>     {
>       "type": "filters.reprojection",
>       "out_srs": "EPSG:4326"
>     },
>
>     {
>
>   "filename": "/home/ayank/Documents/geodata/export/sd_coord_test_4326.tiff",
>       "type": "writers.gdal",
>       "gdaldriver": "GTiff",
>       "dimension": "Z",
>       "output_type": "idw",
>       "resolution": "9.25925925925925925925925e-5"
>     }
>
> Running gdalinfo on the resulting TIFF shows it as not a square:
>
> Size is 552, 541
>
> Corner Coordinates:
> Upper Left  (-117.1805780,  32.7300307) (117d10'50.08"W, 32d43'48.11"N)
> Lower Left  (-117.1805780,  32.6799381) (117d10'50.08"W, 32d40'47.78"N)
> Upper Right (-117.1294669,  32.7300307) (117d 7'46.08"W, 32d43'48.11"N)
> Lower Right (-117.1294669,  32.6799381) (117d 7'46.08"W, 32d40'47.78"N)
> Center      (-117.1550224,  32.7049844) (117d 9'18.08"W, 32d42'17.94”N)
>
> It overlaps with the previous two output files perfectly.
>
> What am I doing wrong?
>
> My thinking has been to keep the LPC data in its original CRS to avoid
> dealing with scaling factors and offsets. It appears to combine NAD83 and
> survey feet.
>
> Thanks,
>
> Andreas
> _______________________________________________
> pdal mailing list
> pdal at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/pdal
>


-- 
Andrew Bell
andrew.bell.ia at gmail.com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/pdal/attachments/20210817/6110d43f/attachment.html>