From 0d911b26c3fb7767ae3f4e9459f2b84d849e988d Mon Sep 17 00:00:00 2001 From: Howard Soh Date: Tue, 19 Nov 2024 19:21:27 +0000 Subject: [PATCH 01/12] #3012 Initial release --- .../libcode/vx_grid/search_3d_kdtree.cc | 328 ++++++++++++++++++ 1 file changed, 328 insertions(+) create mode 100644 internal/test_util/libcode/vx_grid/search_3d_kdtree.cc 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..bf2097cf1f --- /dev/null +++ b/internal/test_util/libcode/vx_grid/search_3d_kdtree.cc @@ -0,0 +1,328 @@ +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* +// ** 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 +// 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); +void llh_to_ecef(double lat, double lon, double alt, double *x_km, double *y_km, double *z_km); +void llh_to_ecef_radian(double lat_r, double lon_r, double alt, 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 << "\n" << argv[0] + << ": Must specify input NetCDF name to process!\n\n"; + if(argc < 3) cerr << "\n" << argv[0] + << ": Must specify latitude to search\n\n"; + if(argc < 4) cerr << "\n" << argv[0] + << ": Must specify longitude to search\n\n"; + if(argc < 5) cerr << "\n" << argv[0] + << ": Must specify altitude to search\n\n"; + + if(argc < 5) { + 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\n"; + 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) { + + 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 llh_to_ecef(double lat, double lon, double alt, double *x_km, double *y_km, double *z_km) { + llh_to_ecef_radian(lat*rad_per_deg, lon*rad_per_deg, alt, x_km, y_km, z_km); +} + +void llh_to_ecef_radian(double lat_r, double lon_r, double alt, double *x_km, double *y_km, double *z_km) { + 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) * cosLat * cos(lon_r) / 1000.; + *y_km = (radius * C + alt) * cosLat * sin(lon_r) / 1000.; + *z_km = (radius * S + alt) * sinLat / 1000.; +} + + +//////////////////////////////////////////////////////////////////////// + +/* +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://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_geodetic_to_ECEF_coordinates + +*/ From acaf115075b7656c27bcff841a3692ed545392f8 Mon Sep 17 00:00:00 2001 From: Howard Soh Date: Tue, 19 Nov 2024 19:21:50 +0000 Subject: [PATCH 02/12] #3012 Added search_3d_kdtree.cc --- .../test_util/libcode/vx_grid/Makefile.am | 21 ++++++- .../test_util/libcode/vx_grid/Makefile.in | 62 ++++++++++++++++--- 2 files changed, 75 insertions(+), 8 deletions(-) diff --git a/internal/test_util/libcode/vx_grid/Makefile.am b/internal/test_util/libcode/vx_grid/Makefile.am index a36bcc50e0..0f5a86f44e 100644 --- a/internal/test_util/libcode/vx_grid/Makefile.am +++ b/internal/test_util/libcode/vx_grid/Makefile.am @@ -10,7 +10,8 @@ include ${top_srcdir}/Make-include # Test programs -noinst_PROGRAMS = test_grid_area +noinst_PROGRAMS = test_grid_area \ + search_3d_kdtree test_grid_area_SOURCES = test_grid_area.cc test_grid_area_CPPFLAGS = ${MET_CPPFLAGS} @@ -27,3 +28,21 @@ 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_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..eb393ec7c2 100644 --- a/internal/test_util/libcode/vx_grid/Makefile.in +++ b/internal/test_util/libcode/vx_grid/Makefile.in @@ -88,7 +88,7 @@ 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) subdir = internal/test_util/libcode/vx_grid ACLOCAL_M4 = $(top_srcdir)/aclocal.m4 am__aclocal_m4_deps = $(top_srcdir)/configure.ac @@ -100,9 +100,16 @@ 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_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 +129,9 @@ 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)/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 +150,8 @@ 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) $(test_grid_area_SOURCES) +DIST_SOURCES = $(search_3d_kdtree_SOURCES) $(test_grid_area_SOURCES) am__can_run_installinfo = \ case $$AM_UPDATE_INFO_DIR in \ n|no|NO) false;; \ @@ -343,6 +352,24 @@ 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_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 +407,10 @@ $(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) + 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 +421,7 @@ 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)/test_grid_area-test_grid_area.Po@am__quote@ # am--include-marker $(am__depfiles_remade): @@ -412,6 +444,20 @@ 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` + 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 +596,8 @@ 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)/test_grid_area-test_grid_area.Po -rm -f Makefile distclean-am: clean-am distclean-compile distclean-generic \ distclean-tags @@ -596,7 +643,8 @@ 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)/test_grid_area-test_grid_area.Po -rm -f Makefile maintainer-clean-am: distclean-am maintainer-clean-generic From 246007e0fe3ff2926b1bc3723ea336d293432faf Mon Sep 17 00:00:00 2001 From: Howard Soh Date: Wed, 11 Dec 2024 16:26:17 +0000 Subject: [PATCH 03/12] #3012 Initial release --- .../libcode/vx_grid/search_3d_kdtree_api.cc | 160 ++++++++++++++++++ 1 file changed, 160 insertions(+) create mode 100644 internal/test_util/libcode/vx_grid/search_3d_kdtree_api.cc 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; + +} From 789551e23de5ed065498ad11823ed2c2bb8f0ba2 Mon Sep 17 00:00:00 2001 From: Howard Soh Date: Wed, 11 Dec 2024 16:28:16 +0000 Subject: [PATCH 04/12] #3012 Added check_llh_to_ecef and test_llh_to_ecef --- .../libcode/vx_grid/search_3d_kdtree.cc | 340 ++++++++++-------- 1 file changed, 185 insertions(+), 155 deletions(-) diff --git a/internal/test_util/libcode/vx_grid/search_3d_kdtree.cc b/internal/test_util/libcode/vx_grid/search_3d_kdtree.cc index bf2097cf1f..f047c58637 100644 --- a/internal/test_util/libcode/vx_grid/search_3d_kdtree.cc +++ b/internal/test_util/libcode/vx_grid/search_3d_kdtree.cc @@ -11,6 +11,8 @@ // // 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]) @@ -49,13 +51,13 @@ using IndexKDTree = atlas::util::IndexKDTree; static Geometry atlas_geometry; -IndexKDTree *build_tree(vector lat, vector lon, vector alt); -void llh_to_ecef(double lat, double lon, double alt, double *x_km, double *y_km, double *z_km); -void llh_to_ecef_radian(double lat_r, double lon_r, double alt, double *x_km, double *y_km, double *z_km); +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[]) { @@ -63,164 +65,160 @@ int main(int argc, char *argv[]) string lon_name = "glon"; string alt_name = "zalt"; - if(argc < 2) cerr << "\n" << argv[0] - << ": Must specify input NetCDF name to process!\n\n"; - if(argc < 3) cerr << "\n" << argv[0] - << ": Must specify latitude to search\n\n"; - if(argc < 4) cerr << "\n" << argv[0] - << ": Must specify longitude to search\n\n"; - if(argc < 5) cerr << "\n" << argv[0] - << ": Must specify altitude to search\n\n"; - + 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\n"; + << ", 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]); + NcFile * _ncFile = open_ncfile(argv[1]); + if (IS_INVALID_NC_P(_ncFile)) { + if (_ncFile) { + delete _ncFile; + _ncFile = (NcFile *)nullptr; + } + exit(1); + } -if(argc > 5) lat_name = argv[5]; -if(argc > 6) lon_name = argv[6]; -if(argc > 7) alt_name = argv[7]; + double lat = std::stod(argv[2]); + double lon = std::stod(argv[3]); + double alt = std::stod(argv[4]); -cout << "Requested lat=" << lat << ", lon=" << lon << ", alt=" << alt - << ", lat_name=" << lat_name << ", lon_name=" << lon_name - << ", alt_name=" << alt_name << "\n\n"; + if(argc > 5) lat_name = argv[5]; + if(argc > 6) lon_name = argv[6]; + if(argc > 7) alt_name = argv[7]; -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); + 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); + 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 (!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; -} + if (_ncFile) { + delete _ncFile; + _ncFile = (NcFile *)nullptr; + } // // done @@ -230,7 +228,9 @@ if (_ncFile) { } -IndexKDTree *build_tree(vector lat, vector lon, vector alt) { +//////////////////////////////////////////////////////////////////////// + +IndexKDTree *build_tree(vector &lat, vector &lon, vector &alt_m) { double x; double y; @@ -240,8 +240,7 @@ IndexKDTree *build_tree(vector lat, vector lon, vector a IndexKDTree *kdtree = new IndexKDTree(atlas_geometry); kdtree->reserve(count); for (int i=0; iinsert(point, n++); } @@ -251,11 +250,27 @@ IndexKDTree *build_tree(vector lat, vector lon, vector a return kdtree; } -void llh_to_ecef(double lat, double lon, double alt, double *x_km, double *y_km, double *z_km) { - llh_to_ecef_radian(lat*rad_per_deg, lon*rad_per_deg, alt, x_km, y_km, 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; + + 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_radian(double lat_r, double lon_r, double alt, double *x_km, double *y_km, double *z_km) { +//////////////////////////////////////////////////////////////////////// + +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); @@ -264,14 +279,28 @@ void llh_to_ecef_radian(double lat_r, double lon_r, double alt, double *x_km, do double C = 1./sqrt(cosLat*cosLat + FF * sinLat*sinLat); double S = C * FF; - *x_km = (radius * C + alt) * cosLat * cos(lon_r) / 1000.; - *y_km = (radius * C + alt) * cosLat * sin(lon_r) / 1000.; - *z_km = (radius * S + alt) * sinLat / 1000.; + *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 @@ -322,6 +351,7 @@ Flattening factor, f depends on the model you assume for your conversion. Typica 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 From 4ed997240d2d8b36211aa9e3014394dcd4db390e Mon Sep 17 00:00:00 2001 From: Howard Soh Date: Wed, 11 Dec 2024 16:30:03 +0000 Subject: [PATCH 05/12] #3012Added search_3d_kdtree_api --- .../test_util/libcode/vx_grid/Makefile.am | 20 ++++++- .../test_util/libcode/vx_grid/Makefile.in | 55 ++++++++++++++++++- 2 files changed, 71 insertions(+), 4 deletions(-) diff --git a/internal/test_util/libcode/vx_grid/Makefile.am b/internal/test_util/libcode/vx_grid/Makefile.am index 0f5a86f44e..5b67c45483 100644 --- a/internal/test_util/libcode/vx_grid/Makefile.am +++ b/internal/test_util/libcode/vx_grid/Makefile.am @@ -11,7 +11,8 @@ include ${top_srcdir}/Make-include # Test programs noinst_PROGRAMS = test_grid_area \ - search_3d_kdtree + search_3d_kdtree \ + search_3d_kdtree_api test_grid_area_SOURCES = test_grid_area.cc test_grid_area_CPPFLAGS = ${MET_CPPFLAGS} @@ -33,6 +34,23 @@ 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 \ diff --git a/internal/test_util/libcode/vx_grid/Makefile.in b/internal/test_util/libcode/vx_grid/Makefile.in index eb393ec7c2..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) search_3d_kdtree$(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 @@ -108,6 +109,13 @@ 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) test_grid_area_DEPENDENCIES = $(am__DEPENDENCIES_1) \ @@ -131,6 +139,7 @@ depcomp = $(SHELL) $(top_srcdir)/depcomp am__maybe_remake_depfiles = depfiles 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@) @@ -150,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 = $(search_3d_kdtree_SOURCES) $(test_grid_area_SOURCES) -DIST_SOURCES = $(search_3d_kdtree_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;; \ @@ -356,6 +367,23 @@ 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 \ @@ -411,6 +439,10 @@ search_3d_kdtree$(EXEEXT): $(search_3d_kdtree_OBJECTS) $(search_3d_kdtree_DEPEND @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) @@ -422,6 +454,7 @@ 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): @@ -458,6 +491,20 @@ search_3d_kdtree-search_3d_kdtree.obj: search_3d_kdtree.cc @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 @@ -597,6 +644,7 @@ clean-am: clean-generic clean-noinstPROGRAMS mostlyclean-am distclean: distclean-am -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 \ @@ -644,6 +692,7 @@ installcheck-am: maintainer-clean: maintainer-clean-am -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 From 44812ae73d35b7c7c666f7cfb281928811e6b5e0 Mon Sep 17 00:00:00 2001 From: Howard Soh Date: Wed, 11 Dec 2024 16:32:17 +0000 Subject: [PATCH 06/12] #3012 Reanmed point_lonlat to points_lonlat. Added points_XYZ & points_XYZ_km --- src/libcode/vx_grid/grid_base.cc | 51 ++++++++++++++------ src/libcode/vx_grid/unstructured_grid_defs.h | 19 ++++++-- 2 files changed, 49 insertions(+), 21 deletions(-) diff --git a/src/libcode/vx_grid/grid_base.cc b/src/libcode/vx_grid/grid_base.cc index d144aa3490..89b9e9a9c5 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: " << (0 < points_lonlat.size() ? "PointLonLat" : "PointXYZ") << "\n" + << " lat_checksum: " << lat_checksum << "\n" + << " lon_checksum: " << lon_checksum << "\n" + << " alt_checksum: " << alt_checksum << "\n" + << " max_distance_km: " << max_distance_km << "\n" ; } @@ -764,7 +765,12 @@ UnstructuredData *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; @@ -1690,13 +1696,26 @@ 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 + && us1->points_XYZ[0] == us2->points_lonlat[0] + && us1->points_XYZ[us1->n_face-1] == us2->points_XYZ[us2->n_face-1]; + && is_eq(us1->alt_checksum, us2->alt_checksum); + } + } } return status; diff --git a/src/libcode/vx_grid/unstructured_grid_defs.h b/src/libcode/vx_grid/unstructured_grid_defs.h index 68456d0d35..f5dab8d20a 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(); @@ -53,10 +53,19 @@ struct UnstructuredData { bool is_in_distance(double distance_km) const; void set_points(int count, double *_lon, double *_lat); void set_points(int count, const std::vector &); + void set_points(int count, double *_lon, double *_lat, 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(); + 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; }; From 6dc8fa7da29ef2f0b8501ef574cf61d62a6fd9e3 Mon Sep 17 00:00:00 2001 From: Howard Soh Date: Wed, 11 Dec 2024 16:42:18 +0000 Subject: [PATCH 07/12] #3012 Compare points_XYZ or points_latlon --- src/libcode/vx_grid/grid_base.cc | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/libcode/vx_grid/grid_base.cc b/src/libcode/vx_grid/grid_base.cc index 89b9e9a9c5..f16ee138a1 100644 --- a/src/libcode/vx_grid/grid_base.cc +++ b/src/libcode/vx_grid/grid_base.cc @@ -1711,10 +1711,15 @@ bool is_eq(const UnstructuredData * us1, const UnstructuredData * us2) } else { status = status - && us1->points_XYZ[0] == us2->points_lonlat[0] - && us1->points_XYZ[us1->n_face-1] == us2->points_XYZ[us2->n_face-1]; + && 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); } + } } } From 981dcc268eca8b23ea14de20d0b0a9573a4b88e2 Mon Sep 17 00:00:00 2001 From: Howard Soh Date: Wed, 11 Dec 2024 16:49:38 +0000 Subject: [PATCH 08/12] #3012 Commented out latlonalt_to_xy & xy_to_latlonalt --- src/libcode/vx_grid/unstructured_grid.h | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/libcode/vx_grid/unstructured_grid.h b/src/libcode/vx_grid/unstructured_grid.h index 3f46913472..786ac53625 100644 --- a/src/libcode/vx_grid/unstructured_grid.h +++ b/src/libcode/vx_grid/unstructured_grid.h @@ -56,6 +56,10 @@ class UnstructuredGrid : public GridRep { virtual void xy_to_latlon(double x, double y, double & lat, double & lon) const; + //void latlonalt_to_xy(double lat, double lon, double alt_m, double & x, double & y) const; + + //void xy_to_latlonalt(double x, double y, double & lat, double & lon, double &alt) const; + virtual double calc_area(int x, int y) const; virtual int nx() const; From f70e4ccfaae7e8a7f3d072f74d5ca0f9cdab03a6 Mon Sep 17 00:00:00 2001 From: Howard Soh Date: Wed, 11 Dec 2024 16:50:51 +0000 Subject: [PATCH 09/12] #3012 Reanmed point_lonlat to points_lonlat. Added points_XYZ & points_XYZ_km. Added llh_to_ecef APIs --- src/libcode/vx_grid/unstructured_grid.cc | 345 +++++++++++++++++++++-- 1 file changed, 318 insertions(+), 27 deletions(-) diff --git a/src/libcode/vx_grid/unstructured_grid.cc b/src/libcode/vx_grid/unstructured_grid.cc index 98cf8bccb9..9958a803e4 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, + 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"; } @@ -143,10 +161,42 @@ void UnstructuredGrid::latlon_to_xy(double lat, double lon, double &x, double &y //////////////////////////////////////////////////////////////////////// +//void UnstructuredGrid::latlonalt_to_xy(double lat, double lon, double alt_m, double &x, double &y) const { +// +// double x_km; +// double y_km; +// double z_km; +// llh_to_ecef(lat, lon, alt_m, &x_km, &y_km, &z_km); +// points_XYZ_km[i] = {x_km, y_km, z_km}; +// PointXYZ _pointXYZ(x_km, y_km, z_km); +// +// IndexKDTree::ValueList neighbor = Data.closest_points(lat, lon, 1, alt_m); +// size_t index(neighbor[0].payload()); +// double distance_km(neighbor[0].distance()/1000.); +// bool in_distance = Data.is_in_distance(distance_km); +// +// x = in_distance ? index : -1.0; +// y = 0; +// +// if(mlog.verbosity_level() >= UGRID_DEBUG_LEVEL) mlog +// << Debug(UGRID_DEBUG_LEVEL) << "UnstructuredGrid::latlonalt_to_xy() " +// << "input(lat,lon,alt)=(" << lat << ", " << lon <<", " << alt +// << ") ==> km(" << x_km << ", " << y_km << ", " << z_km +// << ") ==> (" << x << ", " << y << ") == (" +// << Data.points_XYZ[index].x() << ", " << Data.points_XYZ[index].y() +// << ", " << Data.points_XYZ[index].z() << ") distance= " << distance_km << "km, " +// << _pointXYZ.distance(Data.points_XYZ[index]) +// << " km" << (in_distance ? " " : ", rejected") << "\n"; +//} + + +//////////////////////////////////////////////////////////////////////// + + 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[x].y(); + lon = Data.points_lonlat[x].x(); if(mlog.verbosity_level() >= UGRID_DEBUG_LEVEL) mlog << Debug(UGRID_DEBUG_LEVEL) << "UnstructuredGrid::xy_to_latlon() " @@ -349,22 +399,73 @@ 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++); + 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++); + 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 +474,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 +487,12 @@ void UnstructuredData::copy_from(const UnstructuredData *us_data) { //////////////////////////////////////////////////////////////////////// +bool UnstructuredData::has_PointLatLon() const { + return (points_lonlat.size() > 0); +} + +//////////////////////////////////////////////////////////////////////// + bool UnstructuredData::is_in_distance(double distance_km) const { bool in_distance = is_eq(max_distance_km, bad_data_double) || (max_distance_km <= 0) @@ -397,14 +509,14 @@ void UnstructuredData::set_points(int count, double *_lon, 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 +526,203 @@ 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, double *_lon, double *_lat, double *_alt) { clear_data(); + double x_km; + double y_km; + double z_km; 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, 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) << "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) << 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 (size_t i=0; i < neighbor.size(); i++) { + int index = (int)neighbor[i].payload(); + distance_km = neighbor[i].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 (size_t i=0; i < neighbor.size(); i++) { + int index = (int)neighbor[i].payload(); + distance_km = neighbor[i].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"; + } + //if(mlog.verbosity_level() >= UGRID_DEBUG_LEVEL) { + // 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 (" + // << points_XYZ[0].x() << ", " << points_XYZ[0].y() << ", " << points_XYZ[0].z() + // << ") and (" << points_XYZ[last_i].x() << ", " << points_XYZ[last_i].y() + // << ", " << points_XYZ[last_i].z() << ")\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() { + 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, 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"; + +} From edeb3fd372250f44e8f341e169d4d9afc895acab Mon Sep 17 00:00:00 2001 From: MET Tools Test Account Date: Fri, 13 Dec 2024 18:28:41 +0000 Subject: [PATCH 10/12] Per #3012, update the .gitignore settings to ignore the newly created test executables in the test_util/libcode/vx_grid directory. --- internal/test_util/libcode/vx_grid/.gitignore | 2 ++ 1 file changed, 2 insertions(+) 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 From 24950993d6ae35e979ee97ea034a32674e565872 Mon Sep 17 00:00:00 2001 From: MET Tools Test Account Date: Fri, 13 Dec 2024 20:16:48 +0000 Subject: [PATCH 11/12] 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. --- src/libcode/vx_grid/grid_base.cc | 80 ++++++++++--------- src/libcode/vx_grid/unstructured_grid.cc | 82 ++++++-------------- src/libcode/vx_grid/unstructured_grid.h | 5 -- src/libcode/vx_grid/unstructured_grid_defs.h | 6 +- 4 files changed, 68 insertions(+), 105 deletions(-) diff --git a/src/libcode/vx_grid/grid_base.cc b/src/libcode/vx_grid/grid_base.cc index f16ee138a1..4592dc19b3 100644 --- a/src/libcode/vx_grid/grid_base.cc +++ b/src/libcode/vx_grid/grid_base.cc @@ -337,7 +337,7 @@ void UnstructuredData::dump() const { mlog << Debug(grid_debug_level) << "\nUnstructured Grid Data:\n" << " n_face: " << n_face << "\n" - << " points: " << (0 < points_lonlat.size() ? "PointLonLat" : "PointXYZ") << "\n" + << " points: " << (points_lonlat.empty() ? "PointXYZ" : "PointLonLat") << "\n" << " lat_checksum: " << lat_checksum << "\n" << " lon_checksum: " << lon_checksum << "\n" << " alt_checksum: " << alt_checksum << "\n" @@ -443,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; @@ -556,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; @@ -578,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; @@ -600,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; @@ -622,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; @@ -644,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; @@ -666,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; @@ -688,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; @@ -710,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; @@ -732,7 +732,7 @@ void GridInfo::set(const SemiLatLonData & data) clear(); -SemiLatLonData * D = (SemiLatLonData *) nullptr; +SemiLatLonData * D = nullptr; D = new SemiLatLonData; @@ -743,7 +743,7 @@ D = new SemiLatLonData; *D = data; -sl = D; D = (SemiLatLonData *) nullptr; +sl = D; D = nullptr; return; @@ -760,7 +760,9 @@ 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; @@ -772,7 +774,7 @@ else { D->set_points(data.n_face, data.points_XYZ); } us = D; -D = (UnstructuredData *)nullptr; +D = nullptr; } #endif @@ -1250,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(); @@ -1302,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); @@ -1346,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; // @@ -1732,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 9958a803e4..7a7ec26f98 100644 --- a/src/libcode/vx_grid/unstructured_grid.cc +++ b/src/libcode/vx_grid/unstructured_grid.cc @@ -49,7 +49,7 @@ 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, - string location); + const string &location); //////////////////////////////////////////////////////////////////////// @@ -161,42 +161,10 @@ void UnstructuredGrid::latlon_to_xy(double lat, double lon, double &x, double &y //////////////////////////////////////////////////////////////////////// -//void UnstructuredGrid::latlonalt_to_xy(double lat, double lon, double alt_m, double &x, double &y) const { -// -// double x_km; -// double y_km; -// double z_km; -// llh_to_ecef(lat, lon, alt_m, &x_km, &y_km, &z_km); -// points_XYZ_km[i] = {x_km, y_km, z_km}; -// PointXYZ _pointXYZ(x_km, y_km, z_km); -// -// IndexKDTree::ValueList neighbor = Data.closest_points(lat, lon, 1, alt_m); -// size_t index(neighbor[0].payload()); -// double distance_km(neighbor[0].distance()/1000.); -// bool in_distance = Data.is_in_distance(distance_km); -// -// x = in_distance ? index : -1.0; -// y = 0; -// -// if(mlog.verbosity_level() >= UGRID_DEBUG_LEVEL) mlog -// << Debug(UGRID_DEBUG_LEVEL) << "UnstructuredGrid::latlonalt_to_xy() " -// << "input(lat,lon,alt)=(" << lat << ", " << lon <<", " << alt -// << ") ==> km(" << x_km << ", " << y_km << ", " << z_km -// << ") ==> (" << x << ", " << y << ") == (" -// << Data.points_XYZ[index].x() << ", " << Data.points_XYZ[index].y() -// << ", " << Data.points_XYZ[index].z() << ") distance= " << distance_km << "km, " -// << _pointXYZ.distance(Data.points_XYZ[index]) -// << " km" << (in_distance ? " " : ", rejected") << "\n"; -//} - - -//////////////////////////////////////////////////////////////////////// - - void UnstructuredGrid::xy_to_latlon(double x, double y, double &lat, double &lon) const { - lat = Data.points_lonlat[x].y(); - lon = Data.points_lonlat[x].x(); + lat = (double) Data.points_lonlat[x].y(); + lon = (double) Data.points_lonlat[x].x(); if(mlog.verbosity_level() >= UGRID_DEBUG_LEVEL) mlog << Debug(UGRID_DEBUG_LEVEL) << "UnstructuredGrid::xy_to_latlon() " @@ -207,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.; @@ -215,6 +184,7 @@ double UnstructuredGrid::calc_area(int x, int y) const { } + //////////////////////////////////////////////////////////////////////// @@ -403,7 +373,8 @@ void UnstructuredData::build_tree() { for (int i=0; iinsert(pointLL, n++); + kdtree->insert(pointLL, n); + n++; lat_checksum += (i+1) * points_lonlat[i].y(); lon_checksum += (i+1) * points_lonlat[i].x(); } @@ -417,7 +388,8 @@ void UnstructuredData::build_tree() { llh_to_ecef(points_XYZ[i].y(), points_XYZ[i].x(), points_XYZ[i].z(), &x_km, &y_km, &z_km); points_XYZ_km.push_back({x_km, y_km, z_km}); - kdtree->insert(points_XYZ_km[i], n++); + kdtree->insert(points_XYZ_km[i], n); + n++; lat_checksum += (i+1) * y_km; lon_checksum += (i+1) * x_km; alt_checksum += (i+1) * z_km; @@ -488,7 +460,7 @@ void UnstructuredData::copy_from(const UnstructuredData *us_data) { //////////////////////////////////////////////////////////////////////// bool UnstructuredData::has_PointLatLon() const { - return (points_lonlat.size() > 0); + return (!points_lonlat.empty()); } //////////////////////////////////////////////////////////////////////// @@ -504,7 +476,7 @@ 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(); @@ -526,13 +498,10 @@ void UnstructuredData::set_points(int count, double *_lon, double *_lat) { //////////////////////////////////////////////////////////////////////// -void UnstructuredData::set_points(int count, double *_lon, double *_lat, double *_alt) { +void UnstructuredData::set_points(int count, const double *_lon, const double *_lat, const double *_alt) { clear_data(); - double x_km; - double y_km; - double z_km; n_face = count; points_XYZ.reserve(count); for (int i=0; i= UGRID_DEBUG_LEVEL) { - // 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 (" - // << points_XYZ[0].x() << ", " << points_XYZ[0].y() << ", " << points_XYZ[0].z() - // << ") and (" << points_XYZ[last_i].x() << ", " << points_XYZ[last_i].y() - // << ", " << points_XYZ[last_i].z() << ")\n"; - //} } } @@ -690,7 +650,7 @@ void UnstructuredData::test_kdtree() { // - 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() { +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"); @@ -715,7 +675,7 @@ void llh_to_ecef(double lat, double lon, double alt_m, double *x_km, double *y_k //////////////////////////////////////////////////////////////////////// -void check_llh_to_ecef(double lat, double lon, double alt_m, double true_x_km, double true_y_km, double true_z_km, string location) { +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; @@ -726,3 +686,5 @@ void check_llh_to_ecef(double lat, double lon, double alt_m, double true_x_km, d << ") 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 786ac53625..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); // @@ -56,10 +55,6 @@ class UnstructuredGrid : public GridRep { virtual void xy_to_latlon(double x, double y, double & lat, double & lon) const; - //void latlonalt_to_xy(double lat, double lon, double alt_m, double & x, double & y) const; - - //void xy_to_latlonalt(double x, double y, double & lat, double & lon, double &alt) const; - virtual double calc_area(int x, int y) const; virtual int nx() const; diff --git a/src/libcode/vx_grid/unstructured_grid_defs.h b/src/libcode/vx_grid/unstructured_grid_defs.h index f5dab8d20a..2ec67504a6 100644 --- a/src/libcode/vx_grid/unstructured_grid_defs.h +++ b/src/libcode/vx_grid/unstructured_grid_defs.h @@ -51,9 +51,9 @@ 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, double *_lon, double *_lat, double *_alt); + 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 &); @@ -61,7 +61,7 @@ struct UnstructuredData { void clear_data(); bool has_PointLatLon() const; void test_kdtree(); - void test_llh_to_ecef(); + 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; From c561df534f8a3f7d3c2165385f2bc84fec443709 Mon Sep 17 00:00:00 2001 From: MET Tools Test Account Date: Fri, 13 Dec 2024 20:55:21 +0000 Subject: [PATCH 12/12] Per #3012, couple more minor SonarQube tweaks --- src/libcode/vx_grid/unstructured_grid.cc | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/libcode/vx_grid/unstructured_grid.cc b/src/libcode/vx_grid/unstructured_grid.cc index 7a7ec26f98..1e318f5075 100644 --- a/src/libcode/vx_grid/unstructured_grid.cc +++ b/src/libcode/vx_grid/unstructured_grid.cc @@ -163,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 = (double) Data.points_lonlat[x].y(); - lon = (double) Data.points_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() " @@ -604,8 +604,8 @@ void UnstructuredData::test_kdtree() { lon = points_lonlat[idx].x(); cout << " - search index=" << idx << " (" << lat << ", " << lon << ")\n"; IndexKDTree::ValueList neighbor = closest_points(lon, lat, closest_n); - for (auto &x : neighbor) { - auto index = (int) x.payload(); + 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"; @@ -629,8 +629,8 @@ void UnstructuredData::test_kdtree() { 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 (auto &x : neighbor) { - auto index = (int) x.payload(); + 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 ("