diff --git a/raster/r.in.pdal/info.cpp b/raster/r.in.pdal/info.cpp index e6e62af587b..95cb0dc508b 100644 --- a/raster/r.in.pdal/info.cpp +++ b/raster/r.in.pdal/info.cpp @@ -8,11 +8,19 @@ * Read the COPYING file that comes with GRASS for details. * */ +#include +#include + +#include +#include #include "info.h" -#include -#ifdef PDAL_USE_NOSRS +extern "C" { +#include "projection.h" +} + +#ifdef R_IN_PDAL_USE_NOSRS void get_extent(struct StringList *infiles, double *min_x, double *max_x, double *min_y, double *max_y, double *min_z, double *max_z, bool nosrs) @@ -22,8 +30,9 @@ void get_extent(struct StringList *infiles, double *min_x, double *max_x, #endif { pdal::StageFactory factory; - bool first = 1; - + pdal::SpatialReference spatial_reference; + bool first = true; + double min_x_, max_x_, min_y_, max_y_, min_z_, max_z_; *min_x = *max_x = *min_y = *max_y = *min_z = *max_z = NAN; for (int i = 0; i < infiles->num_items; i++) { @@ -38,7 +47,7 @@ void get_extent(struct StringList *infiles, double *min_x, double *max_x, pdal::Options las_opts; pdal::Option las_opt("filename", infile); las_opts.add(las_opt); -#ifdef PDAL_USE_NOSRS +#ifdef R_IN_PDAL_USE_NOSRS if (nosrs) { pdal::Option nosrs_opt("nosrs", true); las_opts.add(nosrs_opt); @@ -52,35 +61,123 @@ void get_extent(struct StringList *infiles, double *min_x, double *max_x, catch (const std::exception &err) { G_fatal_error(_("PDAL error: %s"), err.what()); } + spatial_reference = las_reader.getSpatialReference(); + std::string dataset_wkt = spatial_reference.getWKT(); + bool proj_match = is_wkt_projection_same_as_loc(dataset_wkt.c_str()); + const pdal::LasHeader &las_header = las_reader.header(); + + min_x_ = las_header.minX(); + min_y_ = las_header.minY(); + min_z_ = las_header.minZ(); + max_x_ = las_header.maxX(); + max_y_ = las_header.maxY(); + max_z_ = las_header.maxZ(); +#ifdef R_IN_PDAL_USE_NOSRS + bool need_to_reproject = !nosrs && !proj_match && + spatial_reference.valid() && + !spatial_reference.empty(); +#else + bool need_to_reproject = !proj_match && spatial_reference.valid() && + !spatial_reference.empty(); +#endif + + if (need_to_reproject) { + get_reprojected_extent(spatial_reference, &min_x_, &max_x_, &min_y_, + &max_y_, &min_z_, &max_z_); + } if (first) { - *min_x = las_header.minX(); - *min_y = las_header.minY(); - *min_z = las_header.minZ(); - *max_x = las_header.maxX(); - *max_y = las_header.maxY(); - *max_z = las_header.maxZ(); - - first = 0; + *min_x = min_x_; + *min_y = min_y_; + *min_z = min_z_; + *max_x = max_x_; + *max_y = max_y_; + *max_z = max_z_; + + first = false; + } + else { + if (*min_x > min_x_) + *min_x = min_x_; + if (*min_y > min_y_) + *min_y = min_y_; + if (*min_z > min_z_) + *min_z = min_z_; + if (*max_x < max_x_) + *max_x = max_x_; + if (*max_y < max_y_) + *max_y = max_y_; + if (*max_z < max_z_) + *max_z = max_z_; + } + } +} +/* Reproject 4 points representing bounding box using PDAL pipeline */ +void get_reprojected_extent(pdal::SpatialReference &spatial_reference, + double *min_x, double *max_x, double *min_y, + double *max_y, double *min_z, double *max_z) +{ + pdal::PointTable table; + table.layout()->registerDim(pdal::Dimension::Id::X); + table.layout()->registerDim(pdal::Dimension::Id::Y); + table.layout()->registerDim(pdal::Dimension::Id::Z); + + pdal::PointViewPtr view(new pdal::PointView(table)); + + view->setField(pdal::Dimension::Id::X, 0, *min_x); + view->setField(pdal::Dimension::Id::Y, 0, *min_y); + view->setField(pdal::Dimension::Id::Z, 0, *min_z); + view->setField(pdal::Dimension::Id::X, 1, *min_x); + view->setField(pdal::Dimension::Id::Y, 1, *max_y); + view->setField(pdal::Dimension::Id::Z, 1, *min_z); + view->setField(pdal::Dimension::Id::X, 2, *max_x); + view->setField(pdal::Dimension::Id::Y, 2, *min_y); + view->setField(pdal::Dimension::Id::Z, 2, *max_z); + view->setField(pdal::Dimension::Id::X, 3, *max_x); + view->setField(pdal::Dimension::Id::Y, 3, *max_y); + view->setField(pdal::Dimension::Id::Z, 3, *max_z); + + pdal::BufferReader reader; + reader.addView(view); + + pdal::StageFactory factory; + pdal::Stage *reproject = factory.createStage("filters.reprojection"); + pdal::Options reproject_options; + reproject_options.add("in_srs", spatial_reference.getWKT()); + reproject_options.add("out_srs", location_projection_as_wkt(false)); + reproject->setOptions(reproject_options); + reproject->setInput(reader); + reproject->prepare(table); + reproject->execute(table); + + for (pdal::PointId i = 0; i < view->size(); ++i) { + + double x = view->getFieldAs(pdal::Dimension::Id::X, i); + double y = view->getFieldAs(pdal::Dimension::Id::Y, i); + double z = view->getFieldAs(pdal::Dimension::Id::Z, i); + if (i == 0) { + *min_x = *max_x = x; + *min_y = *max_y = y; + *min_z = *max_z = z; } else { - if (*min_x > las_header.minX()) - *min_x = las_header.minX(); - if (*min_y > las_header.minY()) - *min_y = las_header.minY(); - if (*min_z > las_header.minZ()) - *min_z = las_header.minZ(); - if (*max_x < las_header.maxX()) - *max_x = las_header.maxX(); - if (*max_y < las_header.maxY()) - *max_y = las_header.maxY(); - if (*max_z < las_header.maxZ()) - *max_z = las_header.maxZ(); + if (*min_x > x) + *min_x = x; + if (*min_y > y) + *min_y = y; + if (*min_z > z) + *min_z = z; + if (*max_x < x) + *max_x = x; + if (*max_y < y) + *max_y = y; + if (*max_z < z) + *max_z = z; } } } -#ifdef PDAL_USE_NOSRS +#ifdef R_IN_PDAL_USE_NOSRS void print_extent(struct StringList *infiles, bool nosrs) #else void print_extent(struct StringList *infiles) @@ -88,7 +185,7 @@ void print_extent(struct StringList *infiles) { double min_x, max_x, min_y, max_y, min_z, max_z; -#ifdef PDAL_USE_NOSRS +#ifdef R_IN_PDAL_USE_NOSRS get_extent(infiles, &min_x, &max_x, &min_y, &max_y, &min_z, &max_z, nosrs); #else get_extent(infiles, &min_x, &max_x, &min_y, &max_y, &min_z, &max_z); @@ -97,7 +194,7 @@ void print_extent(struct StringList *infiles) min_x, min_z, max_z); } -#ifdef PDAL_USE_NOSRS +#ifdef R_IN_PDAL_USE_NOSRS void print_lasinfo(struct StringList *infiles, bool nosrs) #else void print_lasinfo(struct StringList *infiles) @@ -123,7 +220,7 @@ void print_lasinfo(struct StringList *infiles) pdal::Options las_opts; pdal::Option las_opt("filename", infile); las_opts.add(las_opt); -#ifdef PDAL_USE_NOSRS +#ifdef R_IN_PDAL_USE_NOSRS if (nosrs) { pdal::Option nosrs_opt("nosrs", true); las_opts.add(nosrs_opt); diff --git a/raster/r.in.pdal/info.h b/raster/r.in.pdal/info.h index 09736bab543..4626a87af48 100644 --- a/raster/r.in.pdal/info.h +++ b/raster/r.in.pdal/info.h @@ -23,6 +23,7 @@ #include #include #include +#include #if defined(__clang__) #pragma clang diagnostic pop #endif @@ -31,7 +32,7 @@ #if (PDAL_VERSION_MAJOR >= 2 && PDAL_VERSION_MINOR > 4) || \ (PDAL_VERSION_MAJOR == 2 && PDAL_VERSION_MINOR == 4 && \ PDAL_VERSION_PATCH == 3) -#define PDAL_USE_NOSRS 1 +#define R_IN_R_IN_PDAL_USE_NOSRS 1 #endif extern "C" { @@ -40,7 +41,7 @@ extern "C" { #include "string_list.h" } -#ifdef PDAL_USE_NOSRS +#ifdef R_IN_PDAL_USE_NOSRS void get_extent(struct StringList *, double *, double *, double *, double *, double *, double *, bool); void print_extent(struct StringList *, bool); @@ -51,5 +52,7 @@ void get_extent(struct StringList *, double *, double *, double *, double *, void print_extent(struct StringList *); void print_lasinfo(struct StringList *); #endif +void get_reprojected_extent(pdal::SpatialReference &, double *, double *, + double *, double *, double *, double *); #endif // INFO_H diff --git a/raster/r.in.pdal/main.cpp b/raster/r.in.pdal/main.cpp index fe17e80cb08..0c22cbd7c32 100644 --- a/raster/r.in.pdal/main.cpp +++ b/raster/r.in.pdal/main.cpp @@ -250,11 +250,10 @@ int main(int argc, char *argv[]) reproject_flag->key = 'w'; reproject_flag->label = - _("Reproject to project's coordinate system if needed"); + _("Reproject to project's coordinate system if needed [deprecated]"); reproject_flag->description = - _("Reprojects input dataset to the coordinate system of" - " the GRASS project (by default only datasets with" - " matching coordinate system can be imported"); + _("This flag is deprecated and will be removed in a future release. " + "Input dataset will always be reprojected if needed."); reproject_flag->guisection = _("Projection"); // TODO: from the API it seems that also prj file path and proj string will @@ -446,7 +445,7 @@ int main(int argc, char *argv[]) /* If we print extent, there is no need to validate rest of the input */ if (print_extent_flag->answer) { -#ifdef PDAL_USE_NOSRS +#ifdef R_IN_PDAL_USE_NOSRS print_extent(&infiles, over_flag->answer); #else print_extent(&infiles); @@ -455,7 +454,7 @@ int main(int argc, char *argv[]) } if (print_info_flag->answer) { -#ifdef PDAL_USE_NOSRS +#ifdef R_IN_PDAL_USE_NOSRS print_lasinfo(&infiles, over_flag->answer); #else print_lasinfo(&infiles); @@ -463,6 +462,11 @@ int main(int argc, char *argv[]) exit(EXIT_SUCCESS); } + if (reproject_flag->answer) { + G_verbose_message( + _("Flag 'w' is deprecated and will be removed in a future release. " + "Input dataset will always be reprojected if needed.")); + } /* we could use rules but this gives more info and allows continuing */ if (set_region_flag->answer && !(extents_flag->answer || res_opt->answer || base_rast_res_flag->answer)) { @@ -515,7 +519,7 @@ int main(int argc, char *argv[]) if (extents_flag->answer) { double min_x, max_x, min_y, max_y, min_z, max_z; -#ifdef PDAL_USE_NOSRS +#ifdef R_IN_PDAL_USE_NOSRS get_extent(&infiles, &min_x, &max_x, &min_y, &max_y, &min_z, &max_z, over_flag->answer); #else @@ -718,6 +722,8 @@ int main(int argc, char *argv[]) std::vector readers; pdal::StageFactory factory; pdal::MergeFilter merge_filter; + pdal::SpatialReference spatial_reference; + bool need_to_reproject = false; /* loop of input files */ for (int i = 0; i < infiles.num_items; i++) { const char *infile = infiles.items[i]; @@ -730,7 +736,7 @@ int main(int argc, char *argv[]) pdal::Options las_opts; pdal::Option las_opt("filename", infile); las_opts.add(las_opt); -#ifdef PDAL_USE_NOSRS +#ifdef R_IN_PDAL_USE_NOSRS if (over_flag->answer) { pdal::Option nosrs_opt("nosrs", true); las_opts.add(nosrs_opt); @@ -744,8 +750,28 @@ int main(int argc, char *argv[]) infile); reader->setOptions(las_opts); readers.push_back(reader); + + // getting projection is possible only after prepare + if (!over_flag->answer) { + pdal::PointTable table; + reader->prepare(table); + spatial_reference = reader->getSpatialReference(); + if (spatial_reference.empty()) + G_fatal_error(_("The input dataset has undefined projection")); + std::string dataset_wkt = spatial_reference.getWKT(); + bool proj_match = + is_wkt_projection_same_as_loc(dataset_wkt.c_str()); + + if (!proj_match) + need_to_reproject = true; + } merge_filter.setInput(*reader); } + if (over_flag->answer) { + G_important_message(_("Overriding projection check and assuming" + " that the CRS of input matches" + " the project's CRS")); + } // we need to keep pointer to the last stage // merge filter puts the n readers into one stage, @@ -754,7 +780,7 @@ int main(int argc, char *argv[]) pdal::ReprojectionFilter reprojection_filter; // we reproject when requested regardless of the input projection - if (reproject_flag->answer) { + if (need_to_reproject) { G_message(_("Reprojecting the input to the project's CRS")); char *proj_wkt = location_projection_as_wkt(false); @@ -807,24 +833,6 @@ int main(int argc, char *argv[]) G_fatal_error(_("PDAL error: %s"), err.what()); } - // getting projection is possible only after prepare - if (over_flag->answer) { - G_important_message(_("Overriding projection check and assuming" - " that the CRS of input matches" - " the project's CRS")); - } - else if (!reproject_flag->answer) { - pdal::SpatialReference spatial_reference = - merge_filter.getSpatialReference(); - if (spatial_reference.empty()) - G_fatal_error(_("The input dataset has undefined projection")); - std::string dataset_wkt = spatial_reference.getWKT(); - bool proj_match = is_wkt_projection_same_as_loc(dataset_wkt.c_str()); - - if (!proj_match) - wkt_projection_mismatch_report(dataset_wkt.c_str()); - } - G_important_message(_("Running PDAL algorithms...")); // get the layout to see the dimensions diff --git a/raster/r.in.pdal/projection.c b/raster/r.in.pdal/projection.c index a2b17d8d50e..8c177b17f3a 100644 --- a/raster/r.in.pdal/projection.c +++ b/raster/r.in.pdal/projection.c @@ -88,9 +88,6 @@ void projection_mismatch_report(struct Cell_head cellhd, " in the coordinate reference system definitions," " use the -o flag to ignore them and use" " current project definition.\n")); - strcat(error_msg, - _("Consider generating a new project with 'project' parameter" - " from input data set.\n")); G_fatal_error("%s", error_msg); } diff --git a/raster/r.in.pdal/testsuite/data/points_3358.las b/raster/r.in.pdal/testsuite/data/points_3358.las new file mode 100644 index 00000000000..9e48b70f4af Binary files /dev/null and b/raster/r.in.pdal/testsuite/data/points_3358.las differ diff --git a/raster/r.in.pdal/testsuite/data/points_6346.las b/raster/r.in.pdal/testsuite/data/points_6346.las new file mode 100644 index 00000000000..b6390aa410c Binary files /dev/null and b/raster/r.in.pdal/testsuite/data/points_6346.las differ diff --git a/raster/r.in.pdal/testsuite/test_r_in_pdal_extent.py b/raster/r.in.pdal/testsuite/test_r_in_pdal_extent.py new file mode 100644 index 00000000000..ecaa0c31e6e --- /dev/null +++ b/raster/r.in.pdal/testsuite/test_r_in_pdal_extent.py @@ -0,0 +1,181 @@ +import os +import sys +import pathlib +import unittest +import shutil + +from grass.gunittest.case import TestCase +from grass.gunittest.main import test + +import grass.script as gs +from grass.tools import Tools + +# points.csv +# X,Y,Z +# 637510,222000,100 +# 637505,222012,101 +# 637512,222006,102 +# 637500,222004,103 +# 637503,222009,104 +# 637503,222009,106 +# 637504,222008,106 +# 637511,222001,107 +# 637510,222002,108 +# 637504,222006,109 + +# pipeline_3358.json +# { +# "pipeline": [ +# { +# "type": "readers.text", +# "filename": "points.csv" +# }, +# { +# "type": "writers.las", +# "filename": "data/points_3358.las", +# "a_srs": "EPSG:3358" +# } +# ] + +# pipeline_3358_6346.json +# [ +# { +# "type": "readers.las", +# "filename": "data/points_3358.las" +# }, +# { +# "type": "filters.reprojection", +# "out_srs": "EPSG:6346" +# }, +# { +# "type": "writers.las", +# "filename": "data/points_6346.las" +# } +# ] + +# pdal pipeline pipeline_3358.json +# pdal pipeline pipeline_3358_6346.json + + +class ExtentTest(TestCase): + """Testing extent with reprojection""" + + output = "pdal_output" + + @classmethod + def setUpClass(cls): + """Ensures expected computational region and generated data""" + cls.use_temp_region() + cls.runModule("g.region", n=1, s=0, e=1, w=0, res=1) + + cls.data_dir = os.path.join(pathlib.Path(__file__).parent.absolute(), "data") + cls.point_file_6346 = os.path.join(cls.data_dir, "points_6346.las") + cls.point_file_3358 = os.path.join(cls.data_dir, "points_3358.las") + + @classmethod + def tearDownClass(cls): + """Remove the temporary region and generated data""" + cls.del_temp_region() + + def tearDown(self): + """Remove the outputs created by the import + + This is executed after each test run. + """ + self.runModule( + "g.remove", + flags="f", + type="raster", + name=(self.output), + ) + + @unittest.skipIf(shutil.which("r.in.pdal") is None, "Cannot find r.in.pdal") + def test_reprojection_extent(self): + """Test extent is automatically reprojected from epsg 6346""" + tools = Tools() + extent = tools.r_in_pdal(input=self.point_file_6346, flags="g").text + extent = gs.parse_key_val(extent, vsep=" ", val_type=float) + expected = { + "n": 222012.882582, + "s": 222000.743357, + "e": 637511.895399, + "w": 637499.697426, + "b": 100, + "t": 109, + } + if sys.platform == "darwin": + delta = 0.3 + else: + delta = 1e-6 + for key, expected_value in expected.items(): + with self.subTest(key=key): + self.assertAlmostEqual(extent[key], expected_value, delta=delta) + + @unittest.skipIf(shutil.which("r.in.pdal") is None, "Cannot find r.in.pdal") + def test_no_reprojection_needed_extent(self): + """Test extent matches input 3358 data""" + tools = Tools() + extent = tools.r_in_pdal(input=self.point_file_3358, flags="g").text + extent = gs.parse_key_val(extent, vsep=" ", val_type=float) + self.assertAlmostEqual(extent["n"], 222012, places=6) + self.assertAlmostEqual(extent["s"], 222000, places=6) + self.assertAlmostEqual(extent["e"], 637512, places=6) + self.assertAlmostEqual(extent["w"], 637500, places=6) + self.assertAlmostEqual(extent["b"], 100, places=6) + self.assertAlmostEqual(extent["t"], 109, places=6) + + @unittest.skipIf(shutil.which("r.in.pdal") is None, "Cannot find r.in.pdal") + def test_override_extent(self): + """Test extent matches input 3358 data""" + tools = Tools() + extent = tools.r_in_pdal(input=self.point_file_3358, flags="go").text + extent = gs.parse_key_val(extent, vsep=" ", val_type=float) + self.assertAlmostEqual(extent["n"], 222012, places=6) + self.assertAlmostEqual(extent["s"], 222000, places=6) + self.assertAlmostEqual(extent["e"], 637512, places=6) + self.assertAlmostEqual(extent["w"], 637500, places=6) + self.assertAlmostEqual(extent["b"], 100, places=6) + self.assertAlmostEqual(extent["t"], 109, places=6) + + @unittest.skipIf(shutil.which("r.in.pdal") is None, "Cannot find r.in.pdal") + def test_no_reprojection_needed_flags_e(self): + """Test 3358 data is imported correctly""" + tools = Tools() + tools.g_region(n=1, s=0, e=1, w=0, res=1) + tools.r_in_pdal(input=self.point_file_3358, output=self.output, flags="e") + self.assertRasterExists(self.output) + rinfo = tools.r_info(map=self.output, flags="gr", format="json") + self.assertAlmostEqual(222012, rinfo["north"]) + self.assertAlmostEqual(222000, rinfo["south"]) + self.assertAlmostEqual(637512, rinfo["east"]) + self.assertAlmostEqual(637500, rinfo["west"]) + # The extent shows min=100, that point is filtered out + # because it's at the edge, this may not be ideal behavior + self.assertAlmostEqual(101, rinfo["min"]) + self.assertAlmostEqual(109, rinfo["max"]) + + @unittest.skipIf(shutil.which("r.in.pdal") is None, "Cannot find r.in.pdal") + def test_reprojection_needed_flags_en(self): + """Test 6346 data is imported correctly and test n flag""" + tools = Tools() + tools.g_region(n=1, s=0, e=1, w=0, res=1) + tools.r_in_pdal( + input=self.point_file_6346, output=self.output, resolution=2, flags="en" + ) + self.assertRasterExists(self.output) + rinfo = tools.r_info(map=self.output, flags="g", format="json") + region = tools.g_region(flags="p", format="json") + univar = tools.r_univar(map=self.output, format="json") + self.assertAlmostEqual(region["north"], rinfo["north"]) + self.assertAlmostEqual(region["south"], rinfo["south"]) + self.assertAlmostEqual(region["east"], rinfo["east"]) + self.assertAlmostEqual(region["west"], rinfo["west"]) + self.assertAlmostEqual(2, rinfo["nsres"]) + self.assertAlmostEqual(2, rinfo["ewres"]) + # No points here are filtered out + self.assertAlmostEqual(univar["min"], 100) + self.assertAlmostEqual(univar["max"], 109) + + +if __name__ == "__main__": + test()