From c7e4345c384caad0c68c37b80843903897a87697 Mon Sep 17 00:00:00 2001 From: Howard Soh Date: Fri, 13 Dec 2024 14:37:42 -0700 Subject: [PATCH] Feature #3012 ugrid for 3D points (PointXYZ) (#3035) * #3012 Initial release * #3012 Added search_3d_kdtree.cc * #3012 Initial release * #3012 Added check_llh_to_ecef and test_llh_to_ecef * #3012Added search_3d_kdtree_api * #3012 Reanmed point_lonlat to points_lonlat. Added points_XYZ & points_XYZ_km * #3012 Compare points_XYZ or points_latlon * #3012 Commented out latlonalt_to_xy & xy_to_latlonalt * #3012 Reanmed point_lonlat to points_lonlat. Added points_XYZ & points_XYZ_km. Added llh_to_ecef APIs * Per #3012, update the .gitignore settings to ignore the newly created test executables in the test_util/libcode/vx_grid directory. * Per #3012, changes to the vx_grid library to address SonarQube code smells in New Code and hopefully reduce the Overall number of code smells. * Per #3012, couple more minor SonarQube tweaks --------- Co-authored-by: Howard Soh Co-authored-by: MET Tools Test Account --- internal/test_util/libcode/vx_grid/.gitignore | 2 + .../test_util/libcode/vx_grid/Makefile.am | 39 +- .../test_util/libcode/vx_grid/Makefile.in | 111 +++++- .../libcode/vx_grid/search_3d_kdtree.cc | 358 ++++++++++++++++++ .../libcode/vx_grid/search_3d_kdtree_api.cc | 160 ++++++++ src/libcode/vx_grid/grid_base.cc | 134 ++++--- src/libcode/vx_grid/unstructured_grid.cc | 309 +++++++++++++-- src/libcode/vx_grid/unstructured_grid.h | 1 - src/libcode/vx_grid/unstructured_grid_defs.h | 21 +- 9 files changed, 1040 insertions(+), 95 deletions(-) create mode 100644 internal/test_util/libcode/vx_grid/search_3d_kdtree.cc create mode 100644 internal/test_util/libcode/vx_grid/search_3d_kdtree_api.cc diff --git a/internal/test_util/libcode/vx_grid/.gitignore b/internal/test_util/libcode/vx_grid/.gitignore index e09bcd4b7a..b2cd6a7bc9 100644 --- a/internal/test_util/libcode/vx_grid/.gitignore +++ b/internal/test_util/libcode/vx_grid/.gitignore @@ -1,4 +1,6 @@ test_grid_area +search_3d_kdtree +search_3d_kdtree_api *.o *.a .deps diff --git a/internal/test_util/libcode/vx_grid/Makefile.am b/internal/test_util/libcode/vx_grid/Makefile.am index a36bcc50e0..5b67c45483 100644 --- a/internal/test_util/libcode/vx_grid/Makefile.am +++ b/internal/test_util/libcode/vx_grid/Makefile.am @@ -10,7 +10,9 @@ include ${top_srcdir}/Make-include # Test programs -noinst_PROGRAMS = test_grid_area +noinst_PROGRAMS = test_grid_area \ + search_3d_kdtree \ + search_3d_kdtree_api test_grid_area_SOURCES = test_grid_area.cc test_grid_area_CPPFLAGS = ${MET_CPPFLAGS} @@ -27,3 +29,38 @@ test_grid_area_LDADD = \ -lvx_log \ $(UGRID_DEP_LIBS) \ -lz -lm -lproj + +search_3d_kdtree_SOURCES = search_3d_kdtree.cc +search_3d_kdtree_CPPFLAGS = ${MET_CPPFLAGS} +search_3d_kdtree_LDFLAGS = -L. ${MET_LDFLAGS} +search_3d_kdtree_LDADD = \ + $(UGRID_MET_LIBS) \ + -lvx_nc_util \ + -lvx_geodesy \ + -lvx_nav \ + -lvx_util_math \ + -lvx_util \ + -lvx_math \ + -lvx_cal \ + -lvx_log \ + $(UGRID_DEP_LIBS) \ + -lz -lm -lproj \ + -lnetcdf_c++4 -lnetcdf + +search_3d_kdtree_api_SOURCES = search_3d_kdtree_api.cc +search_3d_kdtree_api_CPPFLAGS = ${MET_CPPFLAGS} +search_3d_kdtree_api_LDFLAGS = -L. ${MET_LDFLAGS} +search_3d_kdtree_api_LDADD = \ + $(UGRID_MET_LIBS) \ + -lvx_nc_util \ + -lvx_grid \ + -lvx_geodesy \ + -lvx_nav \ + -lvx_util_math \ + -lvx_util \ + -lvx_math \ + -lvx_cal \ + -lvx_log \ + $(UGRID_DEP_LIBS) \ + -lz -lm -lproj \ + -lnetcdf_c++4 -lnetcdf diff --git a/internal/test_util/libcode/vx_grid/Makefile.in b/internal/test_util/libcode/vx_grid/Makefile.in index 9f7c65ed86..2eba7d2cfa 100644 --- a/internal/test_util/libcode/vx_grid/Makefile.in +++ b/internal/test_util/libcode/vx_grid/Makefile.in @@ -88,7 +88,8 @@ PRE_UNINSTALL = : POST_UNINSTALL = : build_triplet = @build@ host_triplet = @host@ -noinst_PROGRAMS = test_grid_area$(EXEEXT) +noinst_PROGRAMS = test_grid_area$(EXEEXT) search_3d_kdtree$(EXEEXT) \ + search_3d_kdtree_api$(EXEEXT) subdir = internal/test_util/libcode/vx_grid ACLOCAL_M4 = $(top_srcdir)/aclocal.m4 am__aclocal_m4_deps = $(top_srcdir)/configure.ac @@ -100,9 +101,23 @@ CONFIG_HEADER = $(top_builddir)/config.h CONFIG_CLEAN_FILES = CONFIG_CLEAN_VPATH_FILES = PROGRAMS = $(noinst_PROGRAMS) +am_search_3d_kdtree_OBJECTS = \ + search_3d_kdtree-search_3d_kdtree.$(OBJEXT) +search_3d_kdtree_OBJECTS = $(am_search_3d_kdtree_OBJECTS) +am__DEPENDENCIES_1 = +search_3d_kdtree_DEPENDENCIES = $(am__DEPENDENCIES_1) \ + $(am__DEPENDENCIES_1) +search_3d_kdtree_LINK = $(CXXLD) $(AM_CXXFLAGS) $(CXXFLAGS) \ + $(search_3d_kdtree_LDFLAGS) $(LDFLAGS) -o $@ +am_search_3d_kdtree_api_OBJECTS = \ + search_3d_kdtree_api-search_3d_kdtree_api.$(OBJEXT) +search_3d_kdtree_api_OBJECTS = $(am_search_3d_kdtree_api_OBJECTS) +search_3d_kdtree_api_DEPENDENCIES = $(am__DEPENDENCIES_1) \ + $(am__DEPENDENCIES_1) +search_3d_kdtree_api_LINK = $(CXXLD) $(AM_CXXFLAGS) $(CXXFLAGS) \ + $(search_3d_kdtree_api_LDFLAGS) $(LDFLAGS) -o $@ am_test_grid_area_OBJECTS = test_grid_area-test_grid_area.$(OBJEXT) test_grid_area_OBJECTS = $(am_test_grid_area_OBJECTS) -am__DEPENDENCIES_1 = test_grid_area_DEPENDENCIES = $(am__DEPENDENCIES_1) \ $(am__DEPENDENCIES_1) test_grid_area_LINK = $(CXXLD) $(AM_CXXFLAGS) $(CXXFLAGS) \ @@ -122,7 +137,10 @@ am__v_at_1 = DEFAULT_INCLUDES = -I.@am__isrc@ -I$(top_builddir) depcomp = $(SHELL) $(top_srcdir)/depcomp am__maybe_remake_depfiles = depfiles -am__depfiles_remade = ./$(DEPDIR)/test_grid_area-test_grid_area.Po +am__depfiles_remade = \ + ./$(DEPDIR)/search_3d_kdtree-search_3d_kdtree.Po \ + ./$(DEPDIR)/search_3d_kdtree_api-search_3d_kdtree_api.Po \ + ./$(DEPDIR)/test_grid_area-test_grid_area.Po am__mv = mv -f AM_V_lt = $(am__v_lt_@AM_V@) am__v_lt_ = $(am__v_lt_@AM_DEFAULT_V@) @@ -141,8 +159,10 @@ AM_V_CXXLD = $(am__v_CXXLD_@AM_V@) am__v_CXXLD_ = $(am__v_CXXLD_@AM_DEFAULT_V@) am__v_CXXLD_0 = @echo " CXXLD " $@; am__v_CXXLD_1 = -SOURCES = $(test_grid_area_SOURCES) -DIST_SOURCES = $(test_grid_area_SOURCES) +SOURCES = $(search_3d_kdtree_SOURCES) $(search_3d_kdtree_api_SOURCES) \ + $(test_grid_area_SOURCES) +DIST_SOURCES = $(search_3d_kdtree_SOURCES) \ + $(search_3d_kdtree_api_SOURCES) $(test_grid_area_SOURCES) am__can_run_installinfo = \ case $$AM_UPDATE_INFO_DIR in \ n|no|NO) false;; \ @@ -343,6 +363,41 @@ test_grid_area_LDADD = \ $(UGRID_DEP_LIBS) \ -lz -lm -lproj +search_3d_kdtree_SOURCES = search_3d_kdtree.cc +search_3d_kdtree_CPPFLAGS = ${MET_CPPFLAGS} +search_3d_kdtree_LDFLAGS = -L. ${MET_LDFLAGS} +search_3d_kdtree_LDADD = \ + $(UGRID_MET_LIBS) \ + -lvx_nc_util \ + -lvx_geodesy \ + -lvx_nav \ + -lvx_util_math \ + -lvx_util \ + -lvx_math \ + -lvx_cal \ + -lvx_log \ + $(UGRID_DEP_LIBS) \ + -lz -lm -lproj \ + -lnetcdf_c++4 -lnetcdf + +search_3d_kdtree_api_SOURCES = search_3d_kdtree_api.cc +search_3d_kdtree_api_CPPFLAGS = ${MET_CPPFLAGS} +search_3d_kdtree_api_LDFLAGS = -L. ${MET_LDFLAGS} +search_3d_kdtree_api_LDADD = \ + $(UGRID_MET_LIBS) \ + -lvx_nc_util \ + -lvx_grid \ + -lvx_geodesy \ + -lvx_nav \ + -lvx_util_math \ + -lvx_util \ + -lvx_math \ + -lvx_cal \ + -lvx_log \ + $(UGRID_DEP_LIBS) \ + -lz -lm -lproj \ + -lnetcdf_c++4 -lnetcdf + all: all-am .SUFFIXES: @@ -380,6 +435,14 @@ $(am__aclocal_m4_deps): clean-noinstPROGRAMS: -test -z "$(noinst_PROGRAMS)" || rm -f $(noinst_PROGRAMS) +search_3d_kdtree$(EXEEXT): $(search_3d_kdtree_OBJECTS) $(search_3d_kdtree_DEPENDENCIES) $(EXTRA_search_3d_kdtree_DEPENDENCIES) + @rm -f search_3d_kdtree$(EXEEXT) + $(AM_V_CXXLD)$(search_3d_kdtree_LINK) $(search_3d_kdtree_OBJECTS) $(search_3d_kdtree_LDADD) $(LIBS) + +search_3d_kdtree_api$(EXEEXT): $(search_3d_kdtree_api_OBJECTS) $(search_3d_kdtree_api_DEPENDENCIES) $(EXTRA_search_3d_kdtree_api_DEPENDENCIES) + @rm -f search_3d_kdtree_api$(EXEEXT) + $(AM_V_CXXLD)$(search_3d_kdtree_api_LINK) $(search_3d_kdtree_api_OBJECTS) $(search_3d_kdtree_api_LDADD) $(LIBS) + test_grid_area$(EXEEXT): $(test_grid_area_OBJECTS) $(test_grid_area_DEPENDENCIES) $(EXTRA_test_grid_area_DEPENDENCIES) @rm -f test_grid_area$(EXEEXT) $(AM_V_CXXLD)$(test_grid_area_LINK) $(test_grid_area_OBJECTS) $(test_grid_area_LDADD) $(LIBS) @@ -390,6 +453,8 @@ mostlyclean-compile: distclean-compile: -rm -f *.tab.c +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/search_3d_kdtree-search_3d_kdtree.Po@am__quote@ # am--include-marker +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/search_3d_kdtree_api-search_3d_kdtree_api.Po@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/test_grid_area-test_grid_area.Po@am__quote@ # am--include-marker $(am__depfiles_remade): @@ -412,6 +477,34 @@ am--depfiles: $(am__depfiles_remade) @AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ @am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(CXXCOMPILE) -c -o $@ `$(CYGPATH_W) '$<'` +search_3d_kdtree-search_3d_kdtree.o: search_3d_kdtree.cc +@am__fastdepCXX_TRUE@ $(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(search_3d_kdtree_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT search_3d_kdtree-search_3d_kdtree.o -MD -MP -MF $(DEPDIR)/search_3d_kdtree-search_3d_kdtree.Tpo -c -o search_3d_kdtree-search_3d_kdtree.o `test -f 'search_3d_kdtree.cc' || echo '$(srcdir)/'`search_3d_kdtree.cc +@am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) $(DEPDIR)/search_3d_kdtree-search_3d_kdtree.Tpo $(DEPDIR)/search_3d_kdtree-search_3d_kdtree.Po +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ $(AM_V_CXX)source='search_3d_kdtree.cc' object='search_3d_kdtree-search_3d_kdtree.o' libtool=no @AMDEPBACKSLASH@ +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ +@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(search_3d_kdtree_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o search_3d_kdtree-search_3d_kdtree.o `test -f 'search_3d_kdtree.cc' || echo '$(srcdir)/'`search_3d_kdtree.cc + +search_3d_kdtree-search_3d_kdtree.obj: search_3d_kdtree.cc +@am__fastdepCXX_TRUE@ $(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(search_3d_kdtree_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT search_3d_kdtree-search_3d_kdtree.obj -MD -MP -MF $(DEPDIR)/search_3d_kdtree-search_3d_kdtree.Tpo -c -o search_3d_kdtree-search_3d_kdtree.obj `if test -f 'search_3d_kdtree.cc'; then $(CYGPATH_W) 'search_3d_kdtree.cc'; else $(CYGPATH_W) '$(srcdir)/search_3d_kdtree.cc'; fi` +@am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) $(DEPDIR)/search_3d_kdtree-search_3d_kdtree.Tpo $(DEPDIR)/search_3d_kdtree-search_3d_kdtree.Po +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ $(AM_V_CXX)source='search_3d_kdtree.cc' object='search_3d_kdtree-search_3d_kdtree.obj' libtool=no @AMDEPBACKSLASH@ +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ +@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(search_3d_kdtree_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o search_3d_kdtree-search_3d_kdtree.obj `if test -f 'search_3d_kdtree.cc'; then $(CYGPATH_W) 'search_3d_kdtree.cc'; else $(CYGPATH_W) '$(srcdir)/search_3d_kdtree.cc'; fi` + +search_3d_kdtree_api-search_3d_kdtree_api.o: search_3d_kdtree_api.cc +@am__fastdepCXX_TRUE@ $(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(search_3d_kdtree_api_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT search_3d_kdtree_api-search_3d_kdtree_api.o -MD -MP -MF $(DEPDIR)/search_3d_kdtree_api-search_3d_kdtree_api.Tpo -c -o search_3d_kdtree_api-search_3d_kdtree_api.o `test -f 'search_3d_kdtree_api.cc' || echo '$(srcdir)/'`search_3d_kdtree_api.cc +@am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) $(DEPDIR)/search_3d_kdtree_api-search_3d_kdtree_api.Tpo $(DEPDIR)/search_3d_kdtree_api-search_3d_kdtree_api.Po +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ $(AM_V_CXX)source='search_3d_kdtree_api.cc' object='search_3d_kdtree_api-search_3d_kdtree_api.o' libtool=no @AMDEPBACKSLASH@ +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ +@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(search_3d_kdtree_api_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o search_3d_kdtree_api-search_3d_kdtree_api.o `test -f 'search_3d_kdtree_api.cc' || echo '$(srcdir)/'`search_3d_kdtree_api.cc + +search_3d_kdtree_api-search_3d_kdtree_api.obj: search_3d_kdtree_api.cc +@am__fastdepCXX_TRUE@ $(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(search_3d_kdtree_api_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT search_3d_kdtree_api-search_3d_kdtree_api.obj -MD -MP -MF $(DEPDIR)/search_3d_kdtree_api-search_3d_kdtree_api.Tpo -c -o search_3d_kdtree_api-search_3d_kdtree_api.obj `if test -f 'search_3d_kdtree_api.cc'; then $(CYGPATH_W) 'search_3d_kdtree_api.cc'; else $(CYGPATH_W) '$(srcdir)/search_3d_kdtree_api.cc'; fi` +@am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) $(DEPDIR)/search_3d_kdtree_api-search_3d_kdtree_api.Tpo $(DEPDIR)/search_3d_kdtree_api-search_3d_kdtree_api.Po +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ $(AM_V_CXX)source='search_3d_kdtree_api.cc' object='search_3d_kdtree_api-search_3d_kdtree_api.obj' libtool=no @AMDEPBACKSLASH@ +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ +@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(search_3d_kdtree_api_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o search_3d_kdtree_api-search_3d_kdtree_api.obj `if test -f 'search_3d_kdtree_api.cc'; then $(CYGPATH_W) 'search_3d_kdtree_api.cc'; else $(CYGPATH_W) '$(srcdir)/search_3d_kdtree_api.cc'; fi` + test_grid_area-test_grid_area.o: test_grid_area.cc @am__fastdepCXX_TRUE@ $(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(test_grid_area_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT test_grid_area-test_grid_area.o -MD -MP -MF $(DEPDIR)/test_grid_area-test_grid_area.Tpo -c -o test_grid_area-test_grid_area.o `test -f 'test_grid_area.cc' || echo '$(srcdir)/'`test_grid_area.cc @am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) $(DEPDIR)/test_grid_area-test_grid_area.Tpo $(DEPDIR)/test_grid_area-test_grid_area.Po @@ -550,7 +643,9 @@ clean: clean-am clean-am: clean-generic clean-noinstPROGRAMS mostlyclean-am distclean: distclean-am - -rm -f ./$(DEPDIR)/test_grid_area-test_grid_area.Po + -rm -f ./$(DEPDIR)/search_3d_kdtree-search_3d_kdtree.Po + -rm -f ./$(DEPDIR)/search_3d_kdtree_api-search_3d_kdtree_api.Po + -rm -f ./$(DEPDIR)/test_grid_area-test_grid_area.Po -rm -f Makefile distclean-am: clean-am distclean-compile distclean-generic \ distclean-tags @@ -596,7 +691,9 @@ install-ps-am: installcheck-am: maintainer-clean: maintainer-clean-am - -rm -f ./$(DEPDIR)/test_grid_area-test_grid_area.Po + -rm -f ./$(DEPDIR)/search_3d_kdtree-search_3d_kdtree.Po + -rm -f ./$(DEPDIR)/search_3d_kdtree_api-search_3d_kdtree_api.Po + -rm -f ./$(DEPDIR)/test_grid_area-test_grid_area.Po -rm -f Makefile maintainer-clean-am: distclean-am maintainer-clean-generic diff --git a/internal/test_util/libcode/vx_grid/search_3d_kdtree.cc b/internal/test_util/libcode/vx_grid/search_3d_kdtree.cc new file mode 100644 index 0000000000..f047c58637 --- /dev/null +++ b/internal/test_util/libcode/vx_grid/search_3d_kdtree.cc @@ -0,0 +1,358 @@ +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* +// ** Copyright UCAR (c) 1992 - 2024 +// ** University Corporation for Atmospheric Research (UCAR) +// ** National Center for Atmospheric Research (NCAR) +// ** Research Applications Lab (RAL) +// ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* + +// ./search_3d_kdtree /d1/personal/dadriaan/projects/NRL/PyIRI/pyiri_f4_2020.nc -23.114 334.77 100. +// ./search_3d_kdtree /d1/personal/dadriaan/projects/NRL/PyIRI/pyiri_f4_2020.nc 34.0522 -118.40806 0. glat glon zalt +// +// lat=-23.114, lon=334.77, alt=100 (meters) +// ==> X: 5309.352 km, Y: -2501.788 km, Z: -2488.375 km +// lat=34.0522, lon=-118.40806, alt=0. (meters) +// ==> X: -2516.715 km, Y: -4653.00 km, Z: 3551.245 km +// index=0 -9.59874, 287.326, 100] to [812813, 6.22979e+06, 1.09673e+06]) +// ==> X: 1873.072 km, Y: -6004.143 km, Z: -1056.531 km +// index=1152 -27.9711, 107.326, 100] to [-5.30871e+06, -2.98241e+06, -1.89224e+06]) +// ==> X: -1678.837 km, Y: 5381.522 km, Z: -2973.724 km +// index=2303 75.6264, 287.326, 100] to [-803942, -6.16181e+06, 1.43314e+06]) +// ==> X: 473.024 km, Y: -1516.283 km, Z: 6156.59 km + +//////////////////////////////////////////////////////////////////////// + +#include +#include +#include +#include + +#include + +#include "atlas/grid/Grid.h" // PointXYZ +#include "atlas/util/Geometry.h" +#include "atlas/util/KDTree.h" + +#include "vx_util.h" +#include "vx_grid.h" +#include "nav.h" +#include "nc_utils_core.h" +#include "nc_utils.h" + +using namespace std; +using namespace netCDF; + +using PointXYZ = atlas::PointXYZ; +using PointLonLat = atlas::PointLonLat; +using Geometry = atlas::Geometry; +using IndexKDTree = atlas::util::IndexKDTree; + +//////////////////////////////////////////////////////////////////////// + +static Geometry atlas_geometry; + +IndexKDTree *build_tree(vector &lat, vector &lon, vector &alt_m); +void check_llh_to_ecef(); +void llh_to_ecef(double lat, double lon, double alt_m, double *x_km, double *y_km, double *z_km); +void test_llh_to_ecef(double lat, double lon, double alt_m, double true_x_km, double true_y_km, double true_z_km, string location); + +//////////////////////////////////////////////////////////////////////// + +int main(int argc, char *argv[]) +{ + + string lat_name = "glat"; + string lon_name = "glon"; + string alt_name = "zalt"; + + if(argc < 2) cerr << argv[0] << ": Must specify input NetCDF name to process!\n"; + if(argc < 3) cerr << argv[0] << ": Must specify latitude to search\n"; + if(argc < 4) cerr << argv[0] << ": Must specify longitude to search\n"; + if(argc < 5) { + cerr << argv[0] << ": Must specify altitude to search\n"; + cerr << "\n Usage: " << argv[0] + << " netcdf_name lat_deg lon_deg alt_meters \n"; + cerr << " def_lat_var_name=" << lat_name + << ", def_lon_var_name=" << lon_name + << ", def_alt_var_name=" << alt_name << "\n"; + + if(argc < 2) check_llh_to_ecef(); + + exit(1); + } + + NcFile * _ncFile = open_ncfile(argv[1]); + if (IS_INVALID_NC_P(_ncFile)) { + if (_ncFile) { + delete _ncFile; + _ncFile = (NcFile *)nullptr; + } + exit(1); + } + + double lat = std::stod(argv[2]); + double lon = std::stod(argv[3]); + double alt = std::stod(argv[4]); + + if(argc > 5) lat_name = argv[5]; + if(argc > 6) lon_name = argv[6]; + if(argc > 7) alt_name = argv[7]; + + cout << "Requested lat=" << lat << ", lon=" << lon << ", alt=" << alt + << ", lat_name=" << lat_name << ", lon_name=" << lon_name + << ", alt_name=" << alt_name << "\n\n"; + + NcVar lat_var = get_nc_var(_ncFile, lat_name); + NcVar lon_var = get_nc_var(_ncFile, lon_name); + NcVar alt_var = get_nc_var(_ncFile, alt_name); + + const int nlat = get_data_size(&lat_var); + const int nlon = get_data_size(&lon_var); + const int nalt = get_data_size(&alt_var); + vector lat_values(nlat); + vector lon_values(nlon); + vector alt_values(nalt); + + if (!get_nc_data(&lat_var, lat_values.data())) { + cerr << "\n" + << " Fail to read " << lat_name << "\n\n"; + } + if (!get_nc_data(&lon_var, lon_values.data())) { + cerr << "\n" + << " Fail to read " << lon_name << "\n\n"; + } + if (!get_nc_data(&alt_var, alt_values.data())) { + cerr << "\n" + << " Fail to read " << alt_name << "\n\n"; + } + + double x,y,z; + double x_f,y_f,z_f; + double x_m,y_m,z_m; + double x_l,y_l,z_l; + + llh_to_ecef(lat, lon, alt, &x, &y, &z); + PointXYZ target_point(x, y, z); + + int idx = 0; + llh_to_ecef(lat_values[idx], lon_values[idx], alt_values[idx], &x_f, &y_f, &z_f); + PointXYZ target_point_first(x_f, y_f, z_f); + + idx = lat_values.size()/2; + llh_to_ecef(lat_values[idx], lon_values[idx], alt_values[idx], &x_m, &y_m, &z_m); + PointXYZ target_point_mid(x_m, y_m, z_m); + + idx = lat_values.size() - 1; + llh_to_ecef(lat_values[idx], lon_values[idx], alt_values[idx], &x_l, &y_l, &z_l); + PointXYZ target_point_last(x_l, y_l, z_l); + + //bool in_distance; + int closest_n = 5; + double distance_km; + IndexKDTree *kdtree = build_tree(lat_values, lon_values, alt_values); + + idx = 0; + cout << "The first point from the input file: " + << " index=" << idx << " " << lat_values[idx] << ", " << lon_values[idx] << ", " << alt_values[idx] << "] to [" << x_f << ", " << y_f << ", " << z_f << "])\n"; + IndexKDTree::ValueList neighbor = kdtree->closestPoints(target_point_first, closest_n); + for (size_t i=0; i < neighbor.size(); i++) { + size_t index(neighbor[i].payload()); + int index2 = (int)neighbor[i].payload(); + distance_km = neighbor[i].distance() / 1000.; + //in_distance = is_in_distance(distance_km); + llh_to_ecef(lat_values[index], lon_values[index], alt_values[index], &x, &y, &z); + cout << " - " << i << ": index=" << index << " distance=" << distance_km << " from [" + << lat_values[index] << ", " << lon_values[index] << ", " << alt_values[index] << "] (" + << x << ", " << y << ", " << z << ")\n"; + } + cout << "\n"; + + idx = lat_values.size()/2; + cout << "The mid point from the input file: " + << " index=" << idx << " " << lat_values[idx] << ", " << lon_values[idx] << ", " << alt_values[idx] << "] to [" << x_m << ", " << y_m << ", " << z_m << "])\n"; + neighbor = kdtree->closestPoints(target_point_mid, closest_n); + for (size_t i=0; i < neighbor.size(); i++) { + size_t index(neighbor[i].payload()); + int index2 = (int)neighbor[i].payload(); + distance_km = neighbor[i].distance() / 1000.; + //in_distance = is_in_distance(distance_km); + llh_to_ecef(lat_values[index], lon_values[index], alt_values[index], &x, &y, &z); + cout << " - " << i << ": index=" << index << " distance=" << distance_km << " from [" + << lat_values[index] << ", " << lon_values[index] << ", " << alt_values[index] << "] (" + << x << ", " << y << ", " << z << ")\n"; + } + cout << "\n"; + + idx = lat_values.size() - 1; + cout << "The last point from the input file: " + << " index=" << idx << " " << lat_values[idx] << ", " << lon_values[idx] << ", " << alt_values[idx] << "] to [" << x_m << ", " << y_m << ", " << z_m << "])\n"; + neighbor = kdtree->closestPoints(target_point_last, closest_n); + for (size_t i=0; i < neighbor.size(); i++) { + size_t index(neighbor[i].payload()); + int index2 = (int)neighbor[i].payload(); + distance_km = neighbor[i].distance() / 1000.; + //in_distance = is_in_distance(distance_km); + llh_to_ecef(lat_values[index], lon_values[index], alt_values[index], &x, &y, &z); + cout << " - " << i << ": index=" << index << " distance=" << distance_km << " from [" + << lat_values[index] << ", " << lon_values[index] << ", " << alt_values[index] << "] (" + << x << ", " << y << ", " << z << ")\n"; + } + cout << "\n"; + + cout << "The requested point: " + << lat << ", " << lon << ", " << alt + << "] to [" << x << ", " << y << ", " << z << "])\n"; + neighbor = kdtree->closestPoints(target_point, closest_n); + for (size_t i=0; i < neighbor.size(); i++) { + size_t index(neighbor[i].payload()); + int index2 = (int)neighbor[i].payload(); + distance_km = neighbor[i].distance() / 1000.; + //in_distance = is_in_distance(distance_km); + llh_to_ecef(lat_values[index], lon_values[index], alt_values[index], &x, &y, &z); + cout << " - " << i << ": index=" << index << " distance=" << distance_km << " from [" + << lat_values[index] << ", " << lon_values[index] << ", " << alt_values[index] << "] (" + << x << ", " << y << ", " << z << ")\n"; + } + cout << "\n"; + + if (_ncFile) { + delete _ncFile; + _ncFile = (NcFile *)nullptr; + } + + // + // done + // + + return 0; + +} + +//////////////////////////////////////////////////////////////////////// + +IndexKDTree *build_tree(vector &lat, vector &lon, vector &alt_m) { + + double x; + double y; + double z; + int count = lat.size(); + atlas::idx_t n = 0; + IndexKDTree *kdtree = new IndexKDTree(atlas_geometry); + kdtree->reserve(count); + for (int i=0; iinsert(point, n++); + } + + kdtree->build(); + + return kdtree; +} + +//////////////////////////////////////////////////////////////////////// + +void check_llh_to_ecef() { + double lat; + double lon; + double alt_m; + double true_x_km; + double true_y_km; + double true_z_km; + + test_llh_to_ecef(34.0522, -118.40806, 0.0, -2516.715, -4653.003, 3551.245, "LA"); + test_llh_to_ecef(-9.59874, 287.326, 100, 1873.072, -6004.143, -1056.531, "First Point"); + test_llh_to_ecef(-27.9711, 107.326, 100, -1678.837, 5381.522, -2973.724, "Middle Point"); + test_llh_to_ecef(75.6264, 287.326, 100, 473.024, -1516.283, 6156.59, "Last Point"); +} + +//////////////////////////////////////////////////////////////////////// + +void llh_to_ecef(double lat, double lon, double alt_m, double *x_km, double *y_km, double *z_km) { + const double lat_r = lat*rad_per_deg; + const double lon_r = lon*rad_per_deg; + const double radius = 6378137.0; // Radius of the Earth (in meters) + const double flat_factor = 1.0/298.257223563; // Flattening factor WGS84 Model + double cosLat = cos(lat_r); + double sinLat = sin(lat_r); + double FF = (1.0-flat_factor)*(1.0-flat_factor); + double C = 1./sqrt(cosLat*cosLat + FF * sinLat*sinLat); + double S = C * FF; + + *x_km = (radius * C + alt_m) * cosLat * cos(lon_r) / 1000.; + *y_km = (radius * C + alt_m) * cosLat * sin(lon_r) / 1000.; + *z_km = (radius * S + alt_m) * sinLat / 1000.; +} + +//////////////////////////////////////////////////////////////////////// + +void test_llh_to_ecef(double lat, double lon, double alt_m, double true_x_km, double true_y_km, double true_z_km, string location) { + double x_km; + double y_km; + double z_km; + + llh_to_ecef(lat, lon, alt_m, &x_km, &y_km, &z_km); + cout << location << ": (" << lat << ", " << lon << ", " << alt_m + << ") => (" << x_km << ", " << y_km << ", " << z_km + << ") Diff: (" << (true_x_km - x_km) << ", " << (true_y_km - y_km) << ", " << (true_z_km - z_km) << ")\n"; + +} + +//////////////////////////////////////////////////////////////////////// + + +/* +https://stackoverflow.com/questions/10473852/convert-latitude-and-longitude-to-point-in-3d-space + +def llarToWorld(lat, lon, alt, rad): + # see: http://www.mathworks.de/help/toolbox/aeroblks/llatoecefposition.html + f = 0 # flattening + ls = atan((1 - f)**2 * tan(lat)) # lambda + + x = rad * cos(ls) * cos(lon) + alt * cos(lat) * cos(lon) + y = rad * cos(ls) * sin(lon) + alt * cos(lat) * sin(lon) + z = rad * sin(ls) + alt * sin(lat) + + return c4d.Vector(x, y, z) + +def LLHtoECEF(lat, lon, alt): + # see http://www.mathworks.de/help/toolbox/aeroblks/llatoecefposition.html + + rad = np.float64(6378137.0) # Radius of the Earth (in meters) + f = np.float64(1.0/298.257223563) # Flattening factor WGS84 Model + cosLat = np.cos(lat) + sinLat = np.sin(lat) + FF = (1.0-f)**2 + C = 1/np.sqrt(cosLat**2 + FF * sinLat**2) + S = C * FF + + x = (rad * C + alt)*cosLat * np.cos(lon) + y = (rad * C + alt)*cosLat * np.sin(lon) + z = (rad * S + alt)*sinLat + + return (x, y, z) + + +Comparison output: finding ECEF for Los Angeles, CA (34.0522, -118.40806, 0 elevation) +My code: +X = -2516715.36114 meters or -2516.715 km +Y = -4653003.08089 meters or -4653.003 km +Z = 3551245.35929 meters or 3551.245 km + +Your Code: +X = -2514072.72181 meters or -2514.072 km +Y = -4648117.26458 meters or -4648.117 km +Z = 3571424.90261 meters or 3571.424 km + +Although in your earth rotation environment your function will produce right geographic region for display, it will NOT give the right ECEF equivalent coordinates. As you can see some of the parameters vary by as much as 20 KM which is rather a large error. + +Flattening factor, f depends on the model you assume for your conversion. Typical, model is WGS 84; however, there are other models. + +Personally, I like to use this link to Naval Postgraduate School for sanity checks on my conversions. + +https://www.oc.nps.edu/oc2902w/coord/llhxyz.html +https://www.oc.nps.edu/oc2902w/coord/llhxyz.htm + +https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_geodetic_to_ECEF_coordinates + +*/ diff --git a/internal/test_util/libcode/vx_grid/search_3d_kdtree_api.cc b/internal/test_util/libcode/vx_grid/search_3d_kdtree_api.cc new file mode 100644 index 0000000000..9bff64ba21 --- /dev/null +++ b/internal/test_util/libcode/vx_grid/search_3d_kdtree_api.cc @@ -0,0 +1,160 @@ +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* +// ** Copyright UCAR (c) 1992 - 2024 +// ** University Corporation for Atmospheric Research (UCAR) +// ** National Center for Atmospheric Research (NCAR) +// ** Research Applications Lab (RAL) +// ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* + +// ./search_3d_kdtree_api /d1/personal/dadriaan/projects/NRL/PyIRI/pyiri_f4_2020.nc -23.114 334.77 100. +// ./search_3d_kdtree_api /d1/personal/dadriaan/projects/NRL/PyIRI/pyiri_f4_2020.nc 34.0522 -118.40806 0. glat glon zalt +// +// lat=-23.114, lon=334.77, alt=100 (meters) +// ==> X: 5309.352 km, Y: -2501.788 km, Z: -2488.375 km +// lat=34.0522, lon=-118.40806, alt=0. (meters) +// ==> X: -2516.715 km, Y: -4653.00 km, Z: 3551.245 km +// index=0 -9.59874, 287.326, 100] to [812813, 6.22979e+06, 1.09673e+06]) +// ==> X: 1873.072 km, Y: -6004.143 km, Z: -1056.531 km +// index=1152 -27.9711, 107.326, 100] to [-5.30871e+06, -2.98241e+06, -1.89224e+06]) +// ==> X: -1678.837 km, Y: 5381.522 km, Z: -2973.724 km +// index=2303 75.6264, 287.326, 100] to [-803942, -6.16181e+06, 1.43314e+06]) +// ==> X: 473.024 km, Y: -1516.283 km, Z: 6156.59 km + +//////////////////////////////////////////////////////////////////////// + +#include +#include +#include +#include + +#include + +#include "atlas/util/KDTree.h" + +#include "vx_util.h" +#include "vx_grid.h" +#include "nav.h" +#include "nc_utils_core.h" +#include "nc_utils.h" + +using namespace std; +using namespace netCDF; + +using IndexKDTree = atlas::util::IndexKDTree; + +//////////////////////////////////////////////////////////////////////// + + +//////////////////////////////////////////////////////////////////////// + +extern void llh_to_ecef(double lat, double lon, double alt_m, double *x_km, double *y_km, double *z_km); + +//////////////////////////////////////////////////////////////////////// + +int main(int argc, char *argv[]) +{ + + string lat_name = "glat"; + string lon_name = "glon"; + string alt_name = "zalt"; + + if(argc < 2) cerr << argv[0] << ": Must specify input NetCDF name to process!\n"; + if(argc < 3) cerr << argv[0] << ": Must specify latitude to search\n"; + if(argc < 4) cerr << argv[0] << ": Must specify longitude to search\n"; + if(argc < 5) { + cerr << argv[0] << ": Must specify altitude to search\n"; + + cerr << "\n Usage: " << argv[0] + << " netcdf_name lat_deg lon_deg alt_meters \n"; + cerr << " def_lat_var_name=" << lat_name + << ", def_lon_var_name=" << lon_name + << ", def_alt_var_name=" << alt_name << "\n"; + + if(argc < 2) { + UnstructuredData ugrid_data1; + ugrid_data1.test_llh_to_ecef(); + } + exit(1); + + } + + NcFile * _ncFile = open_ncfile(argv[1]); + if (IS_INVALID_NC_P(_ncFile)) { + if (_ncFile) { + delete _ncFile; + _ncFile = (NcFile *)nullptr; + } + exit(1); + } + + double lat = std::stod(argv[2]); + double lon = std::stod(argv[3]); + double alt = std::stod(argv[4]); + + if(argc > 5) lat_name = argv[5]; + if(argc > 6) lon_name = argv[6]; + if(argc > 7) alt_name = argv[7]; + + cout << "Requested lat=" << lat << ", lon=" << lon << ", alt=" << alt + << ", lat_name=" << lat_name << ", lon_name=" << lon_name + << ", alt_name=" << alt_name << "\n\n"; + + NcVar lat_var = get_nc_var(_ncFile, lat_name); + NcVar lon_var = get_nc_var(_ncFile, lon_name); + NcVar alt_var = get_nc_var(_ncFile, alt_name); + + const int nlat = get_data_size(&lat_var); + const int nlon = get_data_size(&lon_var); + const int nalt = get_data_size(&alt_var); + vector lat_values(nlat); + vector lon_values(nlon); + vector alt_values(nalt); + + if (!get_nc_data(&lat_var, lat_values.data())) { + cerr << "\n" + << " Fail to read " << lat_name << "\n\n"; + } + if (!get_nc_data(&lon_var, lon_values.data())) { + cerr << "\n" + << " Fail to read " << lon_name << "\n\n"; + } + if (!get_nc_data(&alt_var, alt_values.data())) { + cerr << "\n" + << " Fail to read " << alt_name << "\n\n"; + } + + UnstructuredData ugrid_data; + ugrid_data.set_points(nlat, lon_values.data(), lat_values.data(), alt_values.data()); + ugrid_data.test_kdtree(); + + double x,y,z; + double distance_km; + int closest_n = 5; + + cout << "The requested point: " + << lat << ", " << lon << ", " << alt + << "] to [" << x << ", " << y << ", " << z << "])\n"; + IndexKDTree::ValueList neighbor = ugrid_data.closest_points(lat, lon, closest_n, alt); + for (size_t i=0; i < neighbor.size(); i++) { + size_t index(neighbor[i].payload()); + int index2 = (int)neighbor[i].payload(); + distance_km = neighbor[i].distance() / 1000.; + llh_to_ecef(lat_values[index], lon_values[index], alt_values[index], &x, &y, &z); + cout << " - " << i << ": index=" << index << " distance=" << distance_km << " from [" + << lat_values[index] << ", " << lon_values[index] << ", " << alt_values[index] << "] (" + << x << ", " << y << ", " << z << ")\n"; + } + cout << "\n"; + + if (_ncFile) { + delete _ncFile; + _ncFile = (NcFile *)nullptr; + } + + // + // done + // + + return 0; + +} diff --git a/src/libcode/vx_grid/grid_base.cc b/src/libcode/vx_grid/grid_base.cc index d144aa3490..4592dc19b3 100644 --- a/src/libcode/vx_grid/grid_base.cc +++ b/src/libcode/vx_grid/grid_base.cc @@ -316,31 +316,32 @@ void UnstructuredData::clear() { } - //////////////////////////////////////////////////////////////////////// void UnstructuredData::clear_data() { n_face = n_node = n_edge = 0; - point_lonlat.clear(); - lat_checksum = lon_checksum = 0.; + points_lonlat.clear(); + points_XYZ.clear(); + points_XYZ_km.clear(); + lat_checksum = lon_checksum = alt_checksum = 0.; if (kdtree) { delete kdtree; kdtree = nullptr; } } - //////////////////////////////////////////////////////////////////////// - void UnstructuredData::dump() const { mlog << Debug(grid_debug_level) << "\nUnstructured Grid Data:\n" << " n_face: " << n_face << "\n" - << " lat_checksum: " << lat_checksum << "\n" - << " lon_checksum: " << lon_checksum << "\n" - << " max_distance_km: " << max_distance_km << "\n" + << " points: " << (points_lonlat.empty() ? "PointXYZ" : "PointLonLat") << "\n" + << " lat_checksum: " << lat_checksum << "\n" + << " lon_checksum: " << lon_checksum << "\n" + << " alt_checksum: " << alt_checksum << "\n" + << " max_distance_km: " << max_distance_km << "\n" ; } @@ -442,18 +443,18 @@ void GridInfo::clear() { -if ( lc ) { delete lc; lc = (const LambertData *) nullptr; }; -if ( st ) { delete st; st = (const StereographicData *) nullptr; }; -if ( ll ) { delete ll; ll = (const LatLonData *) nullptr; }; -if ( rll ) { delete rll; rll = (const RotatedLatLonData *) nullptr; }; -if ( m ) { delete m; m = (const MercatorData *) nullptr; }; -if ( g ) { delete g; g = (const GaussianData *) nullptr; }; -if ( gi ) { delete gi; gi = (const GoesImagerData *) nullptr; }; -if ( la ) { delete la; la = (const LaeaData *) nullptr; }; -if ( tc ) { delete tc; tc = (const TcrmwData *) nullptr; }; -if ( sl ) { delete sl; sl = (const SemiLatLonData *) nullptr; }; +if ( lc ) { delete lc; lc = (const LambertData *) nullptr; } +if ( st ) { delete st; st = (const StereographicData *) nullptr; } +if ( ll ) { delete ll; ll = (const LatLonData *) nullptr; } +if ( rll ) { delete rll; rll = (const RotatedLatLonData *) nullptr; } +if ( m ) { delete m; m = (const MercatorData *) nullptr; } +if ( g ) { delete g; g = (const GaussianData *) nullptr; } +if ( gi ) { delete gi; gi = (const GoesImagerData *) nullptr; } +if ( la ) { delete la; la = (const LaeaData *) nullptr; } +if ( tc ) { delete tc; tc = (const TcrmwData *) nullptr; } +if ( sl ) { delete sl; sl = (const SemiLatLonData *) nullptr; } #ifdef WITH_UGRID -if ( us ) { delete us; us = (const UnstructuredData *) nullptr; }; +if ( us ) { delete us; us = (const UnstructuredData *) nullptr; } #endif return; @@ -555,13 +556,13 @@ void GridInfo::set(const LambertData & data) clear(); -LambertData * D = (LambertData *) nullptr; +LambertData * D = nullptr; D = new LambertData; memcpy(D, &data, sizeof(data)); -lc = D; D = (LambertData *) nullptr; +lc = D; D = nullptr; return; @@ -577,13 +578,13 @@ void GridInfo::set(const StereographicData & data) clear(); -StereographicData * D = (StereographicData *) nullptr; +StereographicData * D = nullptr; D = new StereographicData; memcpy(D, &data, sizeof(data)); -st = D; D = (StereographicData *) nullptr; +st = D; D = nullptr; return; @@ -599,13 +600,13 @@ void GridInfo::set(const LatLonData & data) clear(); -LatLonData * D = (LatLonData *) nullptr; +LatLonData * D = nullptr; D = new LatLonData; memcpy(D, &data, sizeof(data)); -ll = D; D = (LatLonData *) nullptr; +ll = D; D = nullptr; return; @@ -621,13 +622,13 @@ void GridInfo::set(const RotatedLatLonData & data) clear(); -RotatedLatLonData * D = (RotatedLatLonData *) nullptr; +RotatedLatLonData * D = nullptr; D = new RotatedLatLonData; memcpy(D, &data, sizeof(data)); -rll = D; D = (RotatedLatLonData *) nullptr; +rll = D; D = nullptr; return; @@ -643,13 +644,13 @@ void GridInfo::set(const MercatorData & data) clear(); -MercatorData * D = (MercatorData *) nullptr; +MercatorData * D = nullptr; D = new MercatorData; memcpy(D, &data, sizeof(data)); -m = D; D = (MercatorData *) nullptr; +m = D; D = nullptr; return; @@ -665,13 +666,13 @@ void GridInfo::set(const GaussianData & data) clear(); -GaussianData * D = (GaussianData *) nullptr; +GaussianData * D = nullptr; D = new GaussianData; memcpy(D, &data, sizeof(data)); -g = D; D = (GaussianData *) nullptr; +g = D; D = nullptr; return; @@ -687,13 +688,13 @@ void GridInfo::set(const GoesImagerData & data) clear(); -GoesImagerData * D = (GoesImagerData *) nullptr; +GoesImagerData * D = nullptr; D = new GoesImagerData; memcpy(D, &data, sizeof(data)); -gi = D; D = (GoesImagerData *) nullptr; +gi = D; D = nullptr; return; @@ -709,13 +710,13 @@ void GridInfo::set(const LaeaData & data) clear(); -LaeaData * D = (LaeaData *) nullptr; +LaeaData * D = nullptr; D = new LaeaData; memcpy(D, &data, sizeof(data)); -la = D; D = (LaeaData *) nullptr; +la = D; D = nullptr; return; @@ -731,7 +732,7 @@ void GridInfo::set(const SemiLatLonData & data) clear(); -SemiLatLonData * D = (SemiLatLonData *) nullptr; +SemiLatLonData * D = nullptr; D = new SemiLatLonData; @@ -742,7 +743,7 @@ D = new SemiLatLonData; *D = data; -sl = D; D = (SemiLatLonData *) nullptr; +sl = D; D = nullptr; return; @@ -759,14 +760,21 @@ void GridInfo::set(const UnstructuredData & data) clear(); -UnstructuredData *D = new UnstructuredData; +UnstructuredData * D = nullptr; + +D = new UnstructuredData; D->n_edge = data.n_edge; D->n_node = data.n_node; D->max_distance_km = data.max_distance_km; -D->set_points(data.n_face, data.point_lonlat); +if (data.has_PointLatLon()) { + D->set_points(data.n_face, data.points_lonlat); +} +else { + D->set_points(data.n_face, data.points_XYZ); +} us = D; -D = (UnstructuredData *)nullptr; +D = nullptr; } #endif @@ -1244,7 +1252,8 @@ if ( (nx_new < 2) || ( ny_new < 2) ) { } Grid g_new; -double lat_ll, lon_ll; +double lat_ll; +double lon_ll; GridInfo info_new = info(); @@ -1296,7 +1305,8 @@ if ( info_new.lc ) { } else if ( info_new.m ) { MercatorData m_new = *(info_new.m); - double lat_ur, lon_ur; + double lat_ur; + double lon_ur; xy_to_latlon(x_ll + nx_new - 1, y_ll + ny_new - 1, lat_ur, lon_ur); @@ -1340,8 +1350,10 @@ Grid Grid::subset_center(double lat_center, double lon_center, int nx_new, int n // subset_ll does sanity checking on nx_new and ny_new, so we don't have to do it here // -int ix_ll, iy_ll; -double dx_center, dy_center; +int ix_ll; +int iy_ll; +double dx_center; +double dy_center; // @@ -1690,13 +1702,31 @@ bool is_eq(const UnstructuredData * us1, const UnstructuredData * us2) bool status = false; if (us1 && us2) { if (us1 == us2) status = true; - else status = us1->n_face == us2->n_face - && us1->n_node == us2->n_node - && us1->n_edge == us2->n_edge - && us1->point_lonlat[0] == us2->point_lonlat[0] - && (us1->n_face > 0 && us1->point_lonlat[us1->n_face-1] == us2->point_lonlat[us2->n_face-1]) - && is_eq(us1->lat_checksum, us2->lat_checksum) - && is_eq(us1->lon_checksum, us2->lon_checksum); + else { + status = us1->n_face == us2->n_face + && us1->n_node == us2->n_node + && us1->n_edge == us2->n_edge + && us1->has_PointLatLon() == us2->has_PointLatLon() + && is_eq(us1->lat_checksum, us2->lat_checksum) + && is_eq(us1->lon_checksum, us2->lon_checksum); + if (status && (us1->n_face > 0)) { + if (us1->has_PointLatLon()) { + status = status + && us1->points_lonlat[0] == us2->points_lonlat[0] + && us1->points_lonlat[us1->n_face-1] == us2->points_lonlat[us2->n_face-1]; + } + else { + status = status + && is_eq(us1->points_XYZ[0].x(), us2->points_XYZ[0].x()) + && is_eq(us1->points_XYZ[0].y(), us2->points_XYZ[0].y()) + && is_eq(us1->points_XYZ[0].z(), us2->points_XYZ[0].z()) + && is_eq(us1->points_XYZ[us1->n_face-1].x(), us2->points_XYZ[us2->n_face-1].x()) + && is_eq(us1->points_XYZ[us1->n_face-1].y(), us2->points_XYZ[us2->n_face-1].y()) + && is_eq(us1->points_XYZ[us1->n_face-1].z(), us2->points_XYZ[us2->n_face-1].z()) + && is_eq(us1->alt_checksum, us2->alt_checksum); + } + } + } } return status; @@ -1708,7 +1738,7 @@ return status; //////////////////////////////////////////////////////////////////////// -int ll_func(double x_center, int N) +static int ll_func(double x_center, int N) { diff --git a/src/libcode/vx_grid/unstructured_grid.cc b/src/libcode/vx_grid/unstructured_grid.cc index 98cf8bccb9..1e318f5075 100644 --- a/src/libcode/vx_grid/unstructured_grid.cc +++ b/src/libcode/vx_grid/unstructured_grid.cc @@ -28,16 +28,29 @@ using namespace std; +using PointXYZ = atlas::PointXYZ; using PointLonLat = atlas::PointLonLat; using Geometry = atlas::Geometry; using IndexKDTree = atlas::util::IndexKDTree; //////////////////////////////////////////////////////////////////////// +constexpr double EARTH_RADIUS = 6378137.0; // Radius of the Earth (in meters) +constexpr double FLAT_FACTOR = 1.0/298.257223563; // Flattening factor WGS84 Model + + static const int UGRID_DEBUG_LEVEL = 9; static atlas::Geometry atlas_geometry; +//////////////////////////////////////////////////////////////////////// + +void llh_to_ecef(double lat, double lon, double alt_m, + double *x_km, double *y_km, double *z_km); +void check_llh_to_ecef(double lat, double lon, double alt_m, + double true_x_km, double true_y_km, double true_z_km, + const string &location); + //////////////////////////////////////////////////////////////////////// @@ -103,7 +116,12 @@ void UnstructuredGrid::set_from_data(const UnstructuredData &data) { Data.n_node = data.n_node; Data.max_distance_km = data.max_distance_km; - Data.set_points(Nx, data.point_lonlat); + if (data.has_PointLatLon()) { + Data.set_points(Nx, data.points_lonlat); + } + else { + Data.set_points(Nx, data.points_XYZ); + } } @@ -120,9 +138,9 @@ void UnstructuredGrid::set_max_distance_km(double max_distance) { void UnstructuredGrid::latlon_to_xy(double lat, double lon, double &x, double &y) const { - PointLonLat _point_lonlat(lon, lat); + PointLonLat point_lonlat(lon, lat); - IndexKDTree::ValueList neighbor = Data.kdtree->closestPoints(_point_lonlat, 1); + IndexKDTree::ValueList neighbor = Data.closest_points(lat, lon, 1); size_t index(neighbor[0].payload()); double distance_km(neighbor[0].distance()/1000.); bool in_distance = Data.is_in_distance(distance_km); @@ -133,9 +151,9 @@ void UnstructuredGrid::latlon_to_xy(double lat, double lon, double &x, double &y if(mlog.verbosity_level() >= UGRID_DEBUG_LEVEL) mlog << Debug(UGRID_DEBUG_LEVEL) << "UnstructuredGrid::latlon_to_xy() " << "input=(" << lon << ", " << lat << ") ==> (" << x << ", " << y << ") == (" - << Data.point_lonlat[index].x() << ", " << Data.point_lonlat[index].y() + << Data.points_lonlat[index].x() << ", " << Data.points_lonlat[index].y() << ") distance= " << distance_km << "km, " - << _point_lonlat.distance(Data.point_lonlat[index]) + << point_lonlat.distance(Data.points_lonlat[index]) << " degree" << (in_distance ? " " : ", rejected") << "\n"; } @@ -145,8 +163,8 @@ void UnstructuredGrid::latlon_to_xy(double lat, double lon, double &x, double &y void UnstructuredGrid::xy_to_latlon(double x, double y, double &lat, double &lon) const { - lat = Data.point_lonlat[x].y(); - lon = Data.point_lonlat[x].x(); + lat = Data.points_lonlat[nint(x)].y(); + lon = Data.points_lonlat[nint(x)].x(); if(mlog.verbosity_level() >= UGRID_DEBUG_LEVEL) mlog << Debug(UGRID_DEBUG_LEVEL) << "UnstructuredGrid::xy_to_latlon() " @@ -157,6 +175,7 @@ void UnstructuredGrid::xy_to_latlon(double x, double y, double &lat, double &lon //////////////////////////////////////////////////////////////////////// + double UnstructuredGrid::calc_area(int x, int y) const { double area = 0.; @@ -165,6 +184,7 @@ double UnstructuredGrid::calc_area(int x, int y) const { } + //////////////////////////////////////////////////////////////////////// @@ -349,22 +369,75 @@ void UnstructuredData::build_tree() { atlas::idx_t n = 0; kdtree = new IndexKDTree(atlas_geometry); kdtree->reserve(n_face); - for (int i=0; iinsert(pointLL, n++); - lat_checksum += (i+1) * point_lonlat[i].y(); - lon_checksum += (i+1) * point_lonlat[i].x(); + if (has_PointLatLon()) { + for (int i=0; iinsert(pointLL, n); + n++; + lat_checksum += (i+1) * points_lonlat[i].y(); + lon_checksum += (i+1) * points_lonlat[i].x(); + } + } + else { + double x_km; + double y_km; + double z_km; + points_XYZ_km.reserve(n_face); + for (int i=0; iinsert(points_XYZ_km[i], n); + n++; + lat_checksum += (i+1) * y_km; + lon_checksum += (i+1) * x_km; + alt_checksum += (i+1) * z_km; + } } kdtree->build(); + ConcatString cs; + if (get_env("MET_TEST_UGRID_KDTREE", cs)) test_kdtree(); } //////////////////////////////////////////////////////////////////////// +IndexKDTree::ValueList UnstructuredData::closest_points(const double &lat, const double &lon, + const size_t &k, const double &alt_m) const { + const string method_name = "UnstructuredData::closest_points() --> "; + if (has_PointLatLon()) { + PointLonLat point_lonlat(lon, lat); + if (!is_eq(alt_m, bad_data_double)) { + mlog << Debug(4) << method_name << "ignored the altitude (" << alt_m << ")\n"; + } + return kdtree->closestPoints(point_lonlat, k); + } + else { + double x_km; + double y_km; + double z_km; + double _alt_m = alt_m; + if (is_eq(_alt_m, bad_data_double)) { + _alt_m = EARTH_RADIUS; + mlog << Debug(4) << method_name << "Set the altitude to earth radius (" << EARTH_RADIUS << ")\n"; + } + llh_to_ecef(lat, lon, _alt_m, &x_km, &y_km, &z_km); + PointXYZ point_XYZ(x_km, y_km, z_km); + return kdtree->closestPoints(point_XYZ, k); + } +}; + +//////////////////////////////////////////////////////////////////////// + void UnstructuredData::copy_from(const UnstructuredData &us_data) { - set_points(us_data.n_face, us_data.point_lonlat); + if (us_data.has_PointLatLon()) { + set_points(us_data.n_face, us_data.points_lonlat); + } + else { + set_points(us_data.n_face, us_data.points_XYZ); + } n_edge = us_data.n_edge; n_node = us_data.n_node; max_distance_km = us_data.max_distance_km; @@ -373,7 +446,12 @@ void UnstructuredData::copy_from(const UnstructuredData &us_data) { //////////////////////////////////////////////////////////////////////// void UnstructuredData::copy_from(const UnstructuredData *us_data) { - set_points(us_data->n_face, us_data->point_lonlat); + if (us_data->has_PointLatLon()) { + set_points(us_data->n_face, us_data->points_lonlat); + } + else { + set_points(us_data->n_face, us_data->points_XYZ); + } n_edge = us_data->n_edge; n_node = us_data->n_node; max_distance_km = us_data->max_distance_km; @@ -381,6 +459,12 @@ void UnstructuredData::copy_from(const UnstructuredData *us_data) { //////////////////////////////////////////////////////////////////////// +bool UnstructuredData::has_PointLatLon() const { + return (!points_lonlat.empty()); +} + +//////////////////////////////////////////////////////////////////////// + bool UnstructuredData::is_in_distance(double distance_km) const { bool in_distance = is_eq(max_distance_km, bad_data_double) || (max_distance_km <= 0) @@ -392,19 +476,19 @@ bool UnstructuredData::is_in_distance(double distance_km) const { //////////////////////////////////////////////////////////////////////// -void UnstructuredData::set_points(int count, double *_lon, double *_lat) { +void UnstructuredData::set_points(int count, const double *_lon, const double *_lat) { clear_data(); n_face = count; - point_lonlat.reserve(count); + points_lonlat.reserve(count); for (int i=0; i= UGRID_DEBUG_LEVEL) mlog << Debug(UGRID_DEBUG_LEVEL) << "UnstructuredData::set_points(int, double *, double *) first (" - << point_lonlat[0].x() << ", " << point_lonlat[0].y() << ") and last (" - << point_lonlat[count-1].x() << ", " << point_lonlat[count-1].y() << ") from (" + << points_lonlat[0].x() << ", " << points_lonlat[0].y() << ") and last (" + << points_lonlat[count-1].x() << ", " << points_lonlat[count-1].y() << ") from (" << _lon[0] << ", " << _lat[0] << ") and (" << _lon[count-1] << ", " << _lat[count-1] << ")\n"; @@ -414,24 +498,193 @@ void UnstructuredData::set_points(int count, double *_lon, double *_lat) { //////////////////////////////////////////////////////////////////////// -void UnstructuredData::set_points(int count, const std::vector &ptLonLat) { +void UnstructuredData::set_points(int count, const double *_lon, const double *_lat, const double *_alt) { clear_data(); n_face = count; - point_lonlat.reserve(count); + points_XYZ.reserve(count); + for (int i=0; i= UGRID_DEBUG_LEVEL) mlog + << Debug(UGRID_DEBUG_LEVEL) << "UnstructuredData::set_points(int, double *lon, double *lat, double *alt) first (" + << points_XYZ[0].x() << ", " << points_XYZ[0].y() << ", " << points_XYZ[0].z() << ") and last (" + << points_XYZ[count-1].x() << ", " << points_XYZ[count-1].y() << ") from (" + << _lon[0] << ", " << _lat[0] << ", " << _alt[0] << ") and (" + << _lon[count-1] << ", " << _lat[count-1] << ")\n"; + + build_tree(); + +} + +//////////////////////////////////////////////////////////////////////// + +void UnstructuredData::set_points(int count, const std::vector &pointsLL) { + static const string method_name + = "UnstructuredData::set_points(int, std::vector &) -> "; + + clear_data(); + + if (count != pointsLL.size()) { + mlog << Warning << "\n" << method_name + << " The count argument (" << count << ") does not match with vector (" + << pointsLL.size() << ")\n\n"; + } + + n_face = count; + points_lonlat.reserve(count); for (int i=0; i= UGRID_DEBUG_LEVEL) mlog - << Debug(UGRID_DEBUG_LEVEL) << "UnstructuredData::set_points(int, std::vector &) first: (" - << point_lonlat[0].x() << ", " << point_lonlat[0].y() << ") and last (" - << point_lonlat[count-1].x() << ", " << point_lonlat[count-1].y() << ") from (" - << ptLonLat[0].x() << ", " << ptLonLat[0].y() << ") and (" - << ptLonLat[count-1].x() << ", " << ptLonLat[count-1].y() << ")\n"; + << Debug(UGRID_DEBUG_LEVEL) << method_name << " first: (" + << points_lonlat[0].x() << ", " << points_lonlat[0].y() << ") and last (" + << points_lonlat[count-1].x() << ", " << points_lonlat[count-1].y() << ") from (" + << pointsLL[0].x() << ", " << pointsLL[0].y() << ") and (" + << pointsLL[count-1].x() << ", " << pointsLL[count-1].y() << ")\n"; + + build_tree(); + +} + +//////////////////////////////////////////////////////////////////////// + +void UnstructuredData::set_points(int count, const std::vector &pointsXYZ) { + static const string method_name + = "UnstructuredData::set_points(int, std::vector &) -> "; + + clear_data(); + + if (count != pointsXYZ.size()) { + mlog << Warning << "\n" << method_name + << " The count argument (" << count << ") does not match with vector (" + << pointsXYZ.size() << ")\n\n"; + } + + n_face = count; + points_XYZ.reserve(count); + for (int i=0; i= UGRID_DEBUG_LEVEL) { + int last_i = count - 1; + mlog << Debug(UGRID_DEBUG_LEVEL) << method_name + << "first: (" << points_XYZ[0].x() << ", " << points_XYZ[0].y() << ", " + << points_XYZ[0].z() << ") and last (" << points_XYZ[last_i].x() << ", " + << points_XYZ[last_i].y() << ", " << points_XYZ[last_i].z() << ") from (" + << pointsXYZ[0].x() << ", " << pointsXYZ[0].y() << ", " << pointsXYZ[0].z() + << ") and (" << pointsXYZ[last_i].x() << ", " << pointsXYZ[last_i].y() + << ", " << pointsXYZ[last_i].z() << ")\n"; + } build_tree(); } //////////////////////////////////////////////////////////////////////// + +void UnstructuredData::test_kdtree() { + static const string method_name + = "UnstructuredData::test_kdtree() -> "; + + int closest_n = 5; + double lat; + double lon; + double distance_km; + vector indices = {0, n_face/2, (n_face - 1)}; + if (has_PointLatLon()) { + if (n_face != points_lonlat.size()) { + mlog << Warning << "\n" << method_name + << " The count argument (" << n_face << ") does not match with the vector (latlon=" + << points_lonlat.size() << ", XYZ=" << points_XYZ.size() << ")\n\n"; + } + for (int idx : indices) { + lat = points_lonlat[idx].y(); + lon = points_lonlat[idx].x(); + cout << " - search index=" << idx << " (" << lat << ", " << lon << ")\n"; + IndexKDTree::ValueList neighbor = closest_points(lon, lat, closest_n); + for (const auto &x : neighbor) { + int index = x.payload(); + distance_km = x.distance() / 1000.; + cout << " + closest index=" << index << " distance=" << distance_km << " from (" + << points_lonlat[index].y() << ", " << points_lonlat[index].x() << ")\n"; + } + cout << "\n"; + } + } + else { + double x_km; + double y_km; + double z_km; + double alt_m; + if (n_face != points_XYZ.size()) { + mlog << Warning << "\n" << method_name + << " The count argument (" << n_face << ") does not match with the vector (XYZ=" + << points_XYZ.size() << ", latlon=" << points_lonlat.size() << ")\n\n"; + } + for (int idx : indices) { + lat = points_XYZ[idx].y(); + lon = points_XYZ[idx].x(); + alt_m = points_XYZ[idx].z(); + cout << " - search index=" << idx << " (" << lat << ", " << lon << ", " << alt_m << ")\n"; + IndexKDTree::ValueList neighbor = closest_points(lat, lon, closest_n, alt_m); + for (const auto &x : neighbor) { + int index = x.payload(); + distance_km = x.distance() / 1000.; + llh_to_ecef(lat, lon, alt_m, &x_km, &y_km, &z_km); + cout << " + closest index=" << index << " distance=" << distance_km << " from (" + << points_XYZ[index].y() << ", " << points_XYZ[index].x() << ", " << points_XYZ[index].z() + << ") km: [" << x_km << ", " << y_km << ", " << z_km << "]\n"; + } + cout << "\n"; + } + } + +} + +//////////////////////////////////////////////////////////////////////// +// Called by internal/test_util/libcode/vx_grid/search_3d_kdtree_api.cc +// Input: /d1/personal/dadriaan/projects/NRL/PyIRI/pyiri_f4_2020.nc +// - search index= 0 (-9.59874, 287.326, 100) ==> ( 1873.072, -6004.143, -1056.531) km +// - search index=1152 (-27.9711, 107.326, 100) ==> (-1678.837, 5381.522, -2973.724) km +// - search index=2303 ( 75.6264, 287.326, 100) ==> ( 473.024, -1516.283, 6156.59 ) km + +void UnstructuredData::test_llh_to_ecef() const { + check_llh_to_ecef( 34.0522, -118.40806, 0., -2516.715, -4653.003, 3551.245, " LA"); + check_llh_to_ecef(-9.59874, 287.326, 100., 1873.072, -6004.143, -1056.531, " First Point"); + check_llh_to_ecef(-27.9711, 107.326, 100., -1678.837, 5381.522, -2973.724, " Middle Point"); + check_llh_to_ecef( 75.6264, 287.326, 100., 473.024, -1516.283, 6156.59, " Last Point"); +} + +//////////////////////////////////////////////////////////////////////// + +void llh_to_ecef(double lat, double lon, double alt_m, double *x_km, double *y_km, double *z_km) { + const double lat_r = lat*rad_per_deg; + const double lon_r = lon*rad_per_deg; + double cosLat = cos(lat_r); + double sinLat = sin(lat_r); + double FF = (1.0-FLAT_FACTOR)*(1.0-FLAT_FACTOR); + double C = 1./sqrt(cosLat*cosLat + FF * sinLat*sinLat); + double S = C * FF; + + *x_km = (EARTH_RADIUS * C + alt_m) * cosLat * cos(lon_r) / 1000.; + *y_km = (EARTH_RADIUS * C + alt_m) * cosLat * sin(lon_r) / 1000.; + *z_km = (EARTH_RADIUS * S + alt_m) * sinLat / 1000.; +} + +//////////////////////////////////////////////////////////////////////// + +void check_llh_to_ecef(double lat, double lon, double alt_m, double true_x_km, double true_y_km, double true_z_km, const string &location) { + double x_km; + double y_km; + double z_km; + + llh_to_ecef(lat, lon, alt_m, &x_km, &y_km, &z_km); + cout << location << ": (" << lat << ", " << lon << ", " << alt_m + << ") => (" << x_km << ", " << y_km << ", " << z_km + << ") Diff: (" << (true_x_km - x_km) << ", " << (true_y_km - y_km) << ", " << (true_z_km - z_km) << ")\n"; + +} + +//////////////////////////////////////////////////////////////////////// diff --git a/src/libcode/vx_grid/unstructured_grid.h b/src/libcode/vx_grid/unstructured_grid.h index 3f46913472..b21e7d32c8 100644 --- a/src/libcode/vx_grid/unstructured_grid.h +++ b/src/libcode/vx_grid/unstructured_grid.h @@ -45,7 +45,6 @@ class UnstructuredGrid : public GridRep { void clear(); void set_from_data(const UnstructuredData &); - //void set_max_distance_deg(double max_distance); void set_max_distance_km(double max_distance); // diff --git a/src/libcode/vx_grid/unstructured_grid_defs.h b/src/libcode/vx_grid/unstructured_grid_defs.h index 68456d0d35..2ec67504a6 100644 --- a/src/libcode/vx_grid/unstructured_grid_defs.h +++ b/src/libcode/vx_grid/unstructured_grid_defs.h @@ -29,12 +29,9 @@ //////////////////////////////////////////////////////////////////////// - - - struct UnstructuredData { - const char * name; // not allocated + const char *name; // not allocated int n_face; int n_edge; @@ -42,8 +39,11 @@ struct UnstructuredData { double max_distance_km; // This should be set after calling set_points() double lat_checksum; double lon_checksum; + double alt_checksum; - std::vector point_lonlat; + std::vector points_lonlat; + std::vector points_XYZ; // lat_deg, lon_der, alt_meters + std::vector points_XYZ_km; // x_km, y_km, z_km atlas::util::IndexKDTree *kdtree; UnstructuredData(); @@ -51,12 +51,21 @@ struct UnstructuredData { void build_tree(); bool is_in_distance(double distance_km) const; - void set_points(int count, double *_lon, double *_lat); + void set_points(int count, const double *_lon, const double *_lat); void set_points(int count, const std::vector &); + void set_points(int count, const double *_lon, const double *_lat, const double *_alt); + void set_points(int count, const std::vector &); void copy_from(const UnstructuredData *); void copy_from(const UnstructuredData &); void clear(); void clear_data(); + bool has_PointLatLon() const; + void test_kdtree(); + void test_llh_to_ecef() const; + atlas::util::IndexKDTree::ValueList closest_points( + const double &lat, const double &lon, const size_t &k, + const double &alt_m=bad_data_double) const; + void dump() const; };