From b5562ad105795f92ac06bd7c06ba5b33bfb5f04f Mon Sep 17 00:00:00 2001 From: John Halley Gotway Date: Mon, 24 Jul 2023 14:21:22 -0600 Subject: [PATCH] Update main_v11.1-ref after dtcenter/MET#2612 (#2621) Co-authored-by: John Halley Gotway Co-authored-by: Seth Linden Co-authored-by: jprestop Co-authored-by: Daniel Adriaansen Co-authored-by: John and Cindy Co-authored-by: Howard Soh Co-authored-by: George McCabe <23407799+georgemccabe@users.noreply.github.com> Co-authored-by: hsoh-u Co-authored-by: MET Tools Test Account Co-authored-by: Seth Linden Co-authored-by: lisagoodrich <33230218+lisagoodrich@users.noreply.github.com> Co-authored-by: davidalbo Co-authored-by: Lisa Goodrich Co-authored-by: metplus-bot <97135045+metplus-bot@users.noreply.github.com> Co-authored-by: j-opatz <59586397+j-opatz@users.noreply.github.com> Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> Co-authored-by: Jonathan Vigh Co-authored-by: Tracy Hertneky <39317287+hertneky@users.noreply.github.com> Fix Python environment issue (#2407) fix definitions of G172 and G220 based on comments in NOAA-EMC/NCEPLIBS-w3emc#157. (#2406) fix #2380 develop override (#2382) fix #2408 develop empty config (#2410) fix #2390 develop compile zlib (#2404) fix #2412 develop climo (#2422) fix #2437 develop convert (#2439) fix for develop, for #2437, forgot one reference to the search_parent for a dictionary lookup. fix #2452 develop airnow (#2454) fix #2449 develop pdf (#2464) fix #2402 develop sonarqube (#2468) fix #2426 develop buoy (#2475) fix 2518 dtypes appf docs (#2519) fix 2531 compilation errors (#2533) fix #2531 compilation_errors_configure (#2535) fix 2596 main v11.1 rpath compilation (#2614) --- .github/jobs/build_docker_image.sh | 6 + .github/workflows/testing.yml | 4 +- docs/Users_Guide/appendixB.rst | 4 + .../scripts/installation/compile_MET_all.sh | 55 +++- .../installation/config/install_met_env.narya | 27 ++ .../test_unit/xml/unit_plot_data_plane.xml | 4 - src/basic/vx_config/celltype_to_string.cc | 56 ++-- src/libcode/vx_data2d_grib2/data2d_grib2.cc | 149 +++++++-- src/libcode/vx_data2d_nccf/nccf_file.cc | 309 +++++++++++++++++- src/libcode/vx_grid/find_grid_by_name.cc | 19 +- src/libcode/vx_grid/find_grid_by_name.h | 20 +- src/libcode/vx_grid/grid_base.cc | 105 ++---- src/libcode/vx_grid/grid_base.h | 17 +- src/libcode/vx_grid/laea_grid.cc | 208 +++++------- src/libcode/vx_grid/laea_grid.h | 32 +- src/libcode/vx_grid/laea_grid_defs.h | 97 +----- src/libcode/vx_nc_util/grid_output.cc | 128 +------- src/tools/other/pb2nc/pb2nc.cc | 6 +- 18 files changed, 685 insertions(+), 561 deletions(-) create mode 100644 internal/scripts/installation/config/install_met_env.narya diff --git a/.github/jobs/build_docker_image.sh b/.github/jobs/build_docker_image.sh index f1f023e303..81ce302b44 100755 --- a/.github/jobs/build_docker_image.sh +++ b/.github/jobs/build_docker_image.sh @@ -17,3 +17,9 @@ if [ $? != 0 ]; then cat ${GITHUB_WORKSPACE}/docker_build.log exit 1 fi + +# Copy the log directory from the image +id=$(docker create ${DOCKERHUB_TAG}) +time_command docker cp $id:/met/logs met_logs +mv met_logs/*.log ${GITHUB_WORKSPACE}/. +docker rm -v $id diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml index 338b1d2c43..7a9a1f28dc 100644 --- a/.github/workflows/testing.yml +++ b/.github/workflows/testing.yml @@ -88,9 +88,9 @@ jobs: MET_BASE_REPO: ${{ needs.job_control.outputs.met_base_repo }} MET_BASE_TAG: ${{ needs.job_control.outputs.met_base_tag }} - - name: Copy Docker build log into logs directory + - name: Copy all build log files into logs directory if: always() - run: cp ${GITHUB_WORKSPACE}/docker_build.log ${RUNNER_WORKSPACE}/logs/ + run: cp ${GITHUB_WORKSPACE}/*.log ${RUNNER_WORKSPACE}/logs/ - name: Push Docker Image run: .github/jobs/push_docker_image.sh diff --git a/docs/Users_Guide/appendixB.rst b/docs/Users_Guide/appendixB.rst index ea3120c10e..a4d155a04f 100644 --- a/docs/Users_Guide/appendixB.rst +++ b/docs/Users_Guide/appendixB.rst @@ -11,6 +11,8 @@ The following map projections are currently supported in MET: * Lambert Conformal Projection +* Lambert Azimuthal Equal Area Projection + * Polar Stereographic Projection (Northern) * Polar Stereographic Projection (Southern) @@ -46,6 +48,8 @@ As an example of specifying a Lambert grid, suppose you have a northern hemisphe To grid = "lambert 614 428 12.190 -133.459 -95.0 12.19058 6367.47 25.0 N"; +For a Lambert Azimuthal Equal Area grid, grid specification strings are not supported. + For a Polar Stereographic grid, the syntax is .. code-block:: none diff --git a/internal/scripts/installation/compile_MET_all.sh b/internal/scripts/installation/compile_MET_all.sh index 47333ca25e..2adf80a45f 100644 --- a/internal/scripts/installation/compile_MET_all.sh +++ b/internal/scripts/installation/compile_MET_all.sh @@ -15,10 +15,10 @@ # either PYTHON_MODULE or by setting PYTHON_NAME and PYTHON_VERSION: # - PYTHON_MODULE (only used if USE_MODULES=TRUE) - format is the name # of the Python module to load followed by an underscore and then the -# version number (e.g. python_3.8.6, The script will then run "module -# load python/3.8.6") +# version number (e.g. python_3.10.4, The script will then run "module +# load python/3.10.4") # - PYTHON_NAME = python (or e.g. python3, etc.) -# - PYTHON_VERSION = 3.8.6 +# - PYTHON_VERSION = 3.10.4 # # For a description of these and other variables, visit the MET # downloads page under "Sample Script For Compiling External @@ -64,6 +64,14 @@ # MET_HDF, MET_HDFEOS, MET_FREETYPEINC, MET_FREETYPELIB, # MET_CAIROINC, MET_CAIROLIB. # +# Users can speed up the compilation of MET and its dependent libraries +# by adding the following to their environment configuration file: +# export MAKE_ARGS=-j # +# replacing the # with the number of cores to use (integer) or simply +# specifying: +# export MAKE_ARGS=-j +# with no integer argument to start as many processes in parallel as +# possible. #================================================ # print command, run it, then error and exit if non-zero value is returned @@ -276,6 +284,7 @@ echo "COMPILER = $COMPILER" echo "COMPILER_FAMILY = $COMPILER_FAMILY" echo "COMPILER_VERSION = $COMPILER_VERSION" COMPILER_MAJOR_VERSION=`echo $COMPILER_VERSION | cut -d'.' -f1` +COMPILER_MINOR_VERSION=`echo $COMPILER_VERSION | cut -d'.' -f2` echo echo "USE_MODULES = ${USE_MODULES}" @@ -337,6 +346,22 @@ echo "export F77=${F77}" echo "export F90=${F90}" echo +# Figure out what kind of OS is being used +unameOut="$(uname -s)" +case "${unameOut}" in + Linux*) machine=Linux;; + Darwin*) machine=Mac;; + CYGWIN*) machine=Cygwin;; + MINGW*) machine=MinGw;; + *) machine="UNKNOWN:${unameOut}" +esac + +if [[ $machine == "Mac" ]]; then + sed_inline="sed -i ''" +else + sed_inline="sed -i''" +fi + # Load Python module if [ ${USE_MODULES} = "TRUE" ]; then @@ -469,17 +494,19 @@ if [ $COMPILE_G2CLIB -eq 1 ]; then rm -rf ${LIB_DIR}/g2clib/g2clib* tar -xf ${TAR_DIR}/g2clib*.tar -C ${LIB_DIR}/g2clib cd ${LIB_DIR}/g2clib/g2clib* + # Sed commands use double-quotes to support variable expansion. - sed -i "s|INC=.*|INC=-I${LIB_DIR}/include -I${LIB_DIR}/include/jasper|g" makefile + $sed_inline "s|INC=.*|INC=-I${LIB_DIR}/include -I${LIB_DIR}/include/jasper|g" makefile # Allow other compilers besides gcc - sed -i "s|CC=gcc|CC=${CC}|g" makefile + $sed_inline "s|CC=gcc|CC=${CC}|g" makefile # remove -D__64BIT__ flag because compiling with it has # shown issues with GRIB/GRIB2 files that are over 2GB in size # This flag was removed in g2clib 1.6.4 # so this can be removed if the version is updated - sed -i 's/-D__64BIT__//g' makefile + $sed_inline 's/-D__64BIT__//g' makefile + echo "cd `pwd`" # g2clib appears to compile but causes failure compiling MET if -j argument is used # so exclude it from this call @@ -541,7 +568,7 @@ fi # Compile NetCDF if [ $COMPILE_NETCDF -eq 1 ]; then - + echo echo "Compiling HDF5 at `date`" mkdir -p ${LIB_DIR}/hdf5 @@ -563,13 +590,18 @@ if [ $COMPILE_NETCDF -eq 1 ]; then echo "cd `pwd`" run_cmd "./configure --prefix=${LIB_DIR} CFLAGS=-fPIC CXXFLAGS=-fPIC LDFLAGS=-L${LIB_DIR}/lib CPPFLAGS=-I${LIB_DIR}/include > netcdf-c.configure.log 2>&1" run_cmd "make ${MAKE_ARGS} install > netcdf-c.make_install.log 2>&1" - + echo echo "Compiling NetCDF-CXX at `date`" tar -xzf ${TAR_DIR}/netcdf-cxx*.tar.gz -C ${LIB_DIR}/netcdf cd ${LIB_DIR}/netcdf/netcdf-cxx* echo "cd `pwd`" - run_cmd "./configure --prefix=${LIB_DIR} LDFLAGS=-L${LIB_DIR}/lib CPPFLAGS=-I${LIB_DIR}/include > netcdf-cxx.configure.log 2>&1" + configure_lib_args="" + if [[ $machine == "Mac" ]]; then + configure_lib_args="-lhdf5_hl -lhdf5 -lz" + fi + run_cmd "./configure --prefix=${LIB_DIR} LDFLAGS=-L${LIB_DIR}/lib CPPFLAGS=-I${LIB_DIR}/include LIBS=\"${LIBS} ${configure_lib_args}\" > netcdf-cxx.configure.log 2>&1" + run_cmd "make ${MAKE_ARGS} install > netcdf-cxx.make_install.log 2>&1" fi @@ -665,6 +697,11 @@ export MET_PYTHON_BIN_EXE=${MET_PYTHON_BIN_EXE:=${MET_PYTHON}/bin/python3} export MET_PYTHON_LD export MET_PYTHON_CC export LDFLAGS="-Wl,--disable-new-dtags" + +if [[ $machine == "Mac" ]]; then + export LDFLAGS="" +fi + # https://www.gnu.org/software/bash/manual/html_node/Shell-Parameter-Expansion.html # ${parameter:+word} # If parameter is null or unset, nothing is substituted, otherwise the expansion of word is substituted. diff --git a/internal/scripts/installation/config/install_met_env.narya b/internal/scripts/installation/config/install_met_env.narya new file mode 100644 index 0000000000..7a4fd0834d --- /dev/null +++ b/internal/scripts/installation/config/install_met_env.narya @@ -0,0 +1,27 @@ +export TEST_BASE=/Volumes/d1/jpresto/projects/MET/MET_installation/11.1.0 +export COMPILER=gnu_12.3.0 +export MET_SUBDIR=${TEST_BASE} +export MET_TARBALL=v11.1.0.tar.gz +export USE_MODULES=FALSE +export MET_PYTHON=/opt/local +export MET_PYTHON_CC=-I${MET_PYTHON}/Library/Frameworks/Python.framework/Versions/3.10/include/python3.10 +export MET_PYTHON_LD=`python3-config --ldflags --embed` +export SET_D64BIT=FALSE +export MAKE_ARGS=-j +export CXXFLAGS="-std=c++11" +#export FFLAGS="-std=legacy" + +# If you've already compiled these and don't need to compile them again, set the following; +# Otherwise, if you do need to install these libraries, comment out the variables below +export EXTERNAL_LIBS=${TEST_BASE}/external_libs +#export MET_GRIB2CLIB=${EXTERNAL_LIBS}/lib +#export MET_GRIB2CINC=${EXTERNAL_LIBS}/include +#export GRIB2CLIB_NAME=-lgrib2c +#export MET_BUFRLIB=${EXTERNAL_LIBS}/lib +#export BUFRLIB_NAME=-lbufr +#export MET_NETCDF=${EXTERNAL_LIBS} +#export MET_HDF5=${EXTERNAL_LIBS} +#export MET_GSL=${EXTERNAL_LIBS} +#export LIB_JASPER=${EXTERNAL_LIBS}/lib +#export LIB_LIBPNG=${EXTERNAL_LIBS}/lib +#export LIB_Z=${EXTERNAL_LIBS}/lib \ No newline at end of file diff --git a/internal/test_unit/xml/unit_plot_data_plane.xml b/internal/test_unit/xml/unit_plot_data_plane.xml index 7c2a68b0b2..e9e5ae965f 100644 --- a/internal/test_unit/xml/unit_plot_data_plane.xml +++ b/internal/test_unit/xml/unit_plot_data_plane.xml @@ -522,9 +522,6 @@ - &MET_BIN;/plot_data_plane diff --git a/src/basic/vx_config/celltype_to_string.cc b/src/basic/vx_config/celltype_to_string.cc index 6270c79e8d..a1646688b3 100644 --- a/src/basic/vx_config/celltype_to_string.cc +++ b/src/basic/vx_config/celltype_to_string.cc @@ -24,12 +24,12 @@ //////////////////////////////////////////////////////////////////////// -#include +using namespace std; -#include "celltype_to_string.h" +#include -using namespace std; +#include "celltype_to_string.h" //////////////////////////////////////////////////////////////////////// @@ -39,7 +39,7 @@ ConcatString celltype_to_string(const CellType t) { -const char * s = (const char *) nullptr; +const char * s = (const char *) 0; switch ( t ) { @@ -73,7 +73,7 @@ switch ( t ) { } // switch -return ConcatString (s); +return ( ConcatString (s) ); } @@ -85,33 +85,33 @@ bool string_to_celltype(const char * text, CellType & t) { - if ( strcmp(text, "integer" ) == 0 ) { t = integer; return true; } -else if ( strcmp(text, "floating_point" ) == 0 ) { t = floating_point; return true; } -else if ( strcmp(text, "boolean" ) == 0 ) { t = boolean; return true; } -else if ( strcmp(text, "cell_mark" ) == 0 ) { t = cell_mark; return true; } -else if ( strcmp(text, "op_add" ) == 0 ) { t = op_add; return true; } - -else if ( strcmp(text, "op_multiply" ) == 0 ) { t = op_multiply; return true; } -else if ( strcmp(text, "op_divide" ) == 0 ) { t = op_divide; return true; } -else if ( strcmp(text, "op_subtract" ) == 0 ) { t = op_subtract; return true; } -else if ( strcmp(text, "op_power" ) == 0 ) { t = op_power; return true; } -else if ( strcmp(text, "op_square" ) == 0 ) { t = op_square; return true; } - -else if ( strcmp(text, "op_negate" ) == 0 ) { t = op_negate; return true; } -else if ( strcmp(text, "op_store" ) == 0 ) { t = op_store; return true; } -else if ( strcmp(text, "op_recall" ) == 0 ) { t = op_recall; return true; } -else if ( strcmp(text, "identifier" ) == 0 ) { t = identifier; return true; } -else if ( strcmp(text, "user_func" ) == 0 ) { t = user_func; return true; } - -else if ( strcmp(text, "builtin_func" ) == 0 ) { t = builtin_func; return true; } -else if ( strcmp(text, "local_var" ) == 0 ) { t = local_var; return true; } -else if ( strcmp(text, "character_string") == 0 ) { t = character_string; return true; } -else if ( strcmp(text, "no_cell_type" ) == 0 ) { t = no_cell_type; return true; } + if ( strcmp(text, "integer" ) == 0 ) { t = integer; return ( true ); } +else if ( strcmp(text, "floating_point" ) == 0 ) { t = floating_point; return ( true ); } +else if ( strcmp(text, "boolean" ) == 0 ) { t = boolean; return ( true ); } +else if ( strcmp(text, "cell_mark" ) == 0 ) { t = cell_mark; return ( true ); } +else if ( strcmp(text, "op_add" ) == 0 ) { t = op_add; return ( true ); } + +else if ( strcmp(text, "op_multiply" ) == 0 ) { t = op_multiply; return ( true ); } +else if ( strcmp(text, "op_divide" ) == 0 ) { t = op_divide; return ( true ); } +else if ( strcmp(text, "op_subtract" ) == 0 ) { t = op_subtract; return ( true ); } +else if ( strcmp(text, "op_power" ) == 0 ) { t = op_power; return ( true ); } +else if ( strcmp(text, "op_square" ) == 0 ) { t = op_square; return ( true ); } + +else if ( strcmp(text, "op_negate" ) == 0 ) { t = op_negate; return ( true ); } +else if ( strcmp(text, "op_store" ) == 0 ) { t = op_store; return ( true ); } +else if ( strcmp(text, "op_recall" ) == 0 ) { t = op_recall; return ( true ); } +else if ( strcmp(text, "identifier" ) == 0 ) { t = identifier; return ( true ); } +else if ( strcmp(text, "user_func" ) == 0 ) { t = user_func; return ( true ); } + +else if ( strcmp(text, "builtin_func" ) == 0 ) { t = builtin_func; return ( true ); } +else if ( strcmp(text, "local_var" ) == 0 ) { t = local_var; return ( true ); } +else if ( strcmp(text, "character_string") == 0 ) { t = character_string; return ( true ); } +else if ( strcmp(text, "no_cell_type" ) == 0 ) { t = no_cell_type; return ( true ); } // // nope // -return false; +return ( false ); } diff --git a/src/libcode/vx_data2d_grib2/data2d_grib2.cc b/src/libcode/vx_data2d_grib2/data2d_grib2.cc index 04fd761d5c..6a3f40574d 100644 --- a/src/libcode/vx_data2d_grib2/data2d_grib2.cc +++ b/src/libcode/vx_data2d_grib2/data2d_grib2.cc @@ -41,6 +41,7 @@ using namespace std; //////////////////////////////////////////////////////////////////////// double scaled2dbl(int scale_factor, int scale_value); +int parse_int4(g2int); //////////////////////////////////////////////////////////////////////// // @@ -1294,45 +1295,101 @@ void MetGrib2DataFile::read_grib2_grid( gribfield *gfld) { // Lambert Azimuthal Equal Area else if ( gfld->igdtnum == 140 ) { - const g2int * p = gfld->igdtmpl; - - ScanMode = p[16]; + ScanMode = gfld->igdtmpl[16]; - // build an LaeaGrib2Data struct with the projection information - LaeaGrib2Data laea; + // build an LaeaData struct with the projection information + LaeaData laea; laea.name = laea_proj_type; laea.spheroid_name = "Grib template"; - int earth_shape_int = p[0]; - if(earth_shape_int == 4) { - laea.radius_km = 0; - laea.equatorial_radius_km = 0.5*6378.1370; - laea.polar_radius_km = 0.5*6356.752314; - laea.is_sphere = false; - } - else { - mlog << Error << "\nMetGrib2DataFile::read_grib2_grid() -> " - << "unsupported earth shape value of " << earth_shape_int << "!\n\n"; - exit(1); + + // earth shape + // Reference: https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table3-2.shtml + int earth_shape_int = gfld->igdtmpl[0]; + switch(earth_shape_int){ + case 0: + laea.radius_km = 6367.47; + laea.equatorial_radius_km = 0; + laea.polar_radius_km = 0; + laea.is_sphere = true; + break; + + case 1: + laea.radius_km = scaled2dbl(gfld->igdtmpl[1], gfld->igdtmpl[2]) / 1000.0; + laea.equatorial_radius_km = 0; + laea.polar_radius_km = 0; + laea.is_sphere = true; + break; + + case 2: + laea.radius_km = 0; + laea.equatorial_radius_km = 6378.1600; + laea.polar_radius_km = 6356.7750; + laea.is_sphere = false; + break; + + case 3: + laea.radius_km = 0; + laea.equatorial_radius_km = scaled2dbl(gfld->igdtmpl[3], gfld->igdtmpl[4]); + laea.polar_radius_km = scaled2dbl(gfld->igdtmpl[5], gfld->igdtmpl[6]); + laea.is_sphere = false; + break; + + case 4: + laea.radius_km = 0; + laea.equatorial_radius_km = 6378.1370; + laea.polar_radius_km = 6356.752314; + laea.is_sphere = false; + break; + + case 6: + laea.radius_km = 6371.2290; + laea.equatorial_radius_km = 0; + laea.polar_radius_km = 0; + laea.is_sphere = true; + break; + + case 7: + laea.radius_km = 0; + laea.equatorial_radius_km = scaled2dbl(gfld->igdtmpl[3], gfld->igdtmpl[4]) / 1000.0; + laea.polar_radius_km = scaled2dbl(gfld->igdtmpl[5], gfld->igdtmpl[6]) / 1000.0; + laea.is_sphere = false; + break; + + default: + mlog << Error << "\nMetGrib2DataFile::read_grib2_grid() -> " + << "unsupported earth shape value of " << earth_shape_int << "!\n\n"; + exit(1); } - laea.nx = p[7]; - laea.ny = p[8]; - laea.lat_first = (double)p[9] / 1000000.0; - // TODO: Suspect a bug in g2clib - // laea.lon_first = -1.0*rescale_lon( (double)p[10] / 1000000.0 ); + + // + // MET#2565: + // Fix a bug when parsing laea.lon_first and laea.central_lon + // in gfld->igdtmpl[10] and gfld->igdtmpl[12] for GRIB2 UKV data. + // These longitudes are negative which the GRIB2C library does not + // parse properly. The parse_int4() function is a workaround to + // reprocesses the data as unsigned chars and check for negative + // values, mimicing logic from the int4() function in wgrib2. // - // grib_dump: longitudeOfFirstGridPointInDegrees = -17.1171; - // wgrib2: Lon1 2164.600777 - // 3-43=129,5,47,201 - // g2clib-1.6.0 grid_templates.h - // 3.140: Lambert Azimuthal Equal Area Projection - // BAD? {140, 17, 0, {1,1,4,1,4,1,4,4,4,-4,4,4,4,1,4,4,1} }, - // FIX? {140, 17, 0, {1,1,4,1,4,1,4,4,4,-4,4,-4,4,1,4,4,1} }, - cout << "TODO LAEA p[10] should be 17.1171 = "<< (double)p[10] / 1000000.0 << "\n"; - laea.lon_first = 17.1171; - laea.standard_lat = (double)p[11] / 1000000.0; - laea.central_lon = -1.0*rescale_lon( (double)p[12] / 1000000.0 ); - laea.dx_km = (double)p[14] / 1000000.0; - laea.dy_km = (double)p[15] / 1000000.0; + // The grib_dump utility already handles negative longitudes well. + // The wgrib2 utility does NOT when populating gfld->igdtmpl. + // However, the wgrib2 -get_byte and -get_int options do work, + // as shown below: + // wgrib2 ukv_agl_temperature_1.5_12.grib2 \ + // -get_byte 3 43 4 -get_int 3 43 1 \ + // -get_byte 3 51 4 -get_int 3 51 1 + // + // The longitudes are correctly parsed as -17.117129 and -2.5, respectively. + // 1:0:3-43=129,5,47,201:3-43=-17117129:3-51=128,38,37,160:3-51=-2500000 + // + + laea.nx = gfld->igdtmpl[7]; + laea.ny = gfld->igdtmpl[8]; + laea.lat_first = (double)gfld->igdtmpl[9] / 1000000.0; + laea.lon_first = -1.0*rescale_lon( (double)parse_int4(gfld->igdtmpl[10]) / 1000000.0); + laea.standard_lat = (double)gfld->igdtmpl[11] / 1000000.0; + laea.central_lon = -1.0*rescale_lon( (double)parse_int4(gfld->igdtmpl[12]) / 1000000.0); + laea.dx_km = (double)gfld->igdtmpl[14] / 1000000.0; + laea.dy_km = (double)gfld->igdtmpl[15] / 1000000.0; // store the grid information @@ -1567,3 +1624,27 @@ double scaled2dbl(int scale_factor, int scale_value) { } //////////////////////////////////////////////////////////////////////// + +int parse_int4(g2int i) { + unsigned char c[4]; + unsigned long n = i; + + c[0] = (n >> 24) & 0xFF; + c[1] = (n >> 16) & 0xFF; + c[2] = (n >> 8) & 0xFF; + c[3] = n & 0xFF; + + // convert unsigned char to signed integer + int i_val; + if(c[0] & 0x80) { + i_val = -(((c[0] & 0x7f) << 24) + (c[1] << 16) + (c[2] << 8) + c[3]); + } + else { + i_val = (c[0] << 24) + (c[1] << 16) + (c[2] << 8) + c[3]; + } + + return(i_val); +} + +//////////////////////////////////////////////////////////////////////// + diff --git a/src/libcode/vx_data2d_nccf/nccf_file.cc b/src/libcode/vx_data2d_nccf/nccf_file.cc index 5514705a8c..9c562bb74f 100644 --- a/src/libcode/vx_data2d_nccf/nccf_file.cc +++ b/src/libcode/vx_data2d_nccf/nccf_file.cc @@ -1483,10 +1483,299 @@ void NcCfFile::get_grid_mapping_azimuthal_equidistant(const NcVar *grid_mapping_ void NcCfFile::get_grid_mapping_lambert_azimuthal_equal_area(const NcVar *grid_mapping_var) { static const string method_name = "NcCfFile::get_grid_mapping_lambert_azimuthal_equal_area()"; + double x_coord_to_m_cf = 1.0; + double y_coord_to_m_cf = 1.0; + + // Look for the x/y dimensions + + for (int dim_num = 0; dim_num < _numDims; ++dim_num) + { + // Get the standard name for the coordinate variable + + const NcVar coord_var = get_var(_ncFile, _dims[dim_num]->getName().c_str()); + if (IS_INVALID_NC(coord_var)) continue; + + const NcVarAtt *std_name_att = get_nc_att(&coord_var, (string)"standard_name"); + if (IS_INVALID_NC_P(std_name_att)) { + if (std_name_att) delete std_name_att; + continue; + } + + ConcatString dim_std_name; + if (!get_att_value_chars(std_name_att, dim_std_name)) { + if (std_name_att) delete std_name_att; + continue; + } + + if (std_name_att) { + delete std_name_att; + std_name_att = (NcVarAtt *)0; + } + + // See if this is an X or Y dimension + + if ( dim_std_name == x_dim_key_name ) + { + _xDim = _dims[dim_num]; + + x_dim_var_name = GET_NC_NAME_P(_xDim).c_str(); + for (int var_num = 0; var_num < Nvars; ++var_num) + { + if ( Var[var_num].name == x_dim_var_name) + { + _xCoordVar = Var[var_num].var; + break; + } + } + } + + if ( dim_std_name == y_dim_key_name) + { + _yDim = _dims[dim_num]; + + y_dim_var_name = GET_NC_NAME_P(_yDim).c_str(); + for (int var_num = 0; var_num < Nvars; ++var_num) + { + if (Var[var_num].name == y_dim_var_name) + { + _yCoordVar = Var[var_num].var; + break; + } + } + } + + } + + if (_xDim == 0) + { + mlog << Error << "\n" << method_name << " -> " + << "Didn't find X dimension (projection_x_coordinate) in netCDF file.\n\n"; + exit(1); + } + + if (_yDim == 0) + { + mlog << Error << "\n" << method_name << " -> " + << "Didn't find Y dimension (projection_y_coordinate) in netCDF file.\n\n"; + exit(1); + } + + if (_xCoordVar == 0) + { + mlog << Error << "\n" << method_name << " -> " + << "Didn't find X coord variable (" << GET_NC_NAME_P(_xDim) + << ") in netCDF file.\n\n"; + exit(1); + } + + if (_yCoordVar == 0) + { + mlog << Error << "\n" << method_name << " -> " + << "Didn't find Y coord variable (" << GET_NC_NAME_P(_yDim) + << ") in netCDF file.\n\n"; + exit(1); + } + + if (get_data_size(_xCoordVar) != (int) GET_NC_SIZE_P(_xDim) || + get_data_size(_yCoordVar) != (int) GET_NC_SIZE_P(_yDim)) + { + mlog << Error << "\n" << method_name << " -> " + << "Coordinate variables don't match dimension sizes in netCDF file.\n\n"; + exit(1); + } + + // Handle coordinate variable units + + ConcatString x_coord_units_name; + if (!get_var_units(_xCoordVar, x_coord_units_name)) { + mlog << Warning << "\n" << method_name << " -> " + << "Units not given for X coordinate variable -- " + << "assuming meters.\n\n"; + } + else { + if (0 == x_coord_units_name.length()) { + mlog << Warning << "\n" << method_name << " -> " + << "Cannot extract X coordinate units from netCDF file -- " + << "assuming meters.\n\n"; + } + else { + if ( x_coord_units_name == "m" || + x_coord_units_name == "meter" || + x_coord_units_name == "meters") x_coord_to_m_cf = 1.0; + else if (x_coord_units_name == "km") x_coord_to_m_cf = 1000.0; + else { + mlog << Error << "\n" << method_name << " -> " + << "The X coordinates must be in meters or kilometers for MET.\n\n"; + exit(1); + } + } + } + + ConcatString y_coord_units_name; + if (!get_var_units(_yCoordVar, y_coord_units_name)) { + mlog << Warning << "\n" << method_name << " -> " + << "Units not given for Y coordinate variable -- " + << "assuming meters.\n\n"; + } + else { + if (0 == y_coord_units_name.length()) { + mlog << Warning << "\n" << method_name << " -> " + << "Cannot extract Y coordinate units from netCDF file -- " + << "assuming meters.\n\n"; + } + else { + if ( y_coord_units_name == "m" || + y_coord_units_name == "meter" || + y_coord_units_name == "meters" ) y_coord_to_m_cf = 1.0; + else if (y_coord_units_name == "km") y_coord_to_m_cf = 1000.0; + else { + mlog << Error << "\n" << method_name << " -> " + << "The X coordinates must be in meters or kilometers for MET.\n\n"; + exit(1); + } + } + } + + // Figure out the dx/dy and x/y pin values from the dimension variables + + long x_counts = GET_NC_SIZE_P(_xDim); + double x_values[x_counts]; + + get_nc_data(_xCoordVar, x_values); + + long y_counts = GET_NC_SIZE_P(_yDim); + double y_values[y_counts]; + + get_nc_data(_yCoordVar, y_values); + + // Unit conversion + + for (int i = 0; i DELTA_TOLERANCE) + { + mlog << Error << "\n" << method_name << " -> " + << "MET can only process Lambert Azimuthal Equal Area files " + << "where the delta along the x-axis is constant (" + << curr_delta << " != " << dx_m_a << ")\n\n"; + exit(1); + } + } + + for (int i = 1; i < (int)y_counts; ++i) + { + double curr_delta = fabs(y_values[i] - y_values[i-1]); + if (fabs(curr_delta - dy_m_a) > DELTA_TOLERANCE) + { + mlog << Error << "\n" << method_name << " -> " + << "MET can only process Lambert Azimuthal Equal Area files " + << "where the delta along the y-axis is constant (" + << curr_delta << " != " << dy_m_a << ")\n\n"; + exit(1); + } + } + + // Fill in the data structure. Remember to negate the longitude + // values since MET uses the mathematical coordinate system centered on + // the center of the earth rather than the regular map coordinate system. + + LaeaNetcdfData data; + data.name = laea_proj_type; + + // longitude_of_projection_origin (convert degrees east to west) + + data.prime_meridian_lon = get_nc_var_att_double( + grid_mapping_var, "longitude_of_prime_meridian", false); + + if(is_bad_data(data.prime_meridian_lon)) data.prime_meridian_lon = 0.0; + else data.prime_meridian_lon *= -1.0; + + // semi_major_axis (convert m to km) + + data.semi_major_axis_km = get_nc_var_att_double( + grid_mapping_var, "semi_major_axis", false); + + if(is_bad_data(data.semi_major_axis_km)) data.semi_major_axis_km = EARTH_MAJOR_AXIS_km; + else data.semi_major_axis_km /= m_per_km; + + // semi_minor_axis (convert m to km) + + data.semi_minor_axis_km = get_nc_var_att_double( + grid_mapping_var, "semi_minor_axis", false); + + if(is_bad_data(data.semi_minor_axis_km)) data.semi_minor_axis_km = EARTH_MAJOR_AXIS_km; + else data.semi_minor_axis_km /= m_per_km; + + // latitude_of_projection_origin + + data.proj_origin_lat = get_nc_var_att_double( + grid_mapping_var, "latitude_of_projection_origin"); + + // longitude_of_projection_origin (convert degrees east to west) + + data.proj_origin_lon = -1.0 * get_nc_var_att_double( + grid_mapping_var, "longitude_of_projection_origin"); + + // false_easting + + double false_easting = get_nc_var_att_double( + grid_mapping_var, "false_easting", false); + + if(!is_bad_data(false_easting) && !is_eq(false_easting, 0.0)) + { + mlog << Error << "\n" << method_name << " -> " + << "MET cannot process Lambert Azimuthal Equal Area files " + << "with non-zero false_easting (" << false_easting + << ").\n\n"; + exit(1); + } + + // false_northing + + double false_northing = get_nc_var_att_double( + grid_mapping_var, "false_northing", false); + + if(!is_bad_data(false_northing) && !is_eq(false_northing, 0.0)) + { + mlog << Error << "\n" << method_name << " -> " + << "MET cannot process Lambert Azimuthal Equal Area files " + << "with non-zero false_northing (" << false_northing + << ").\n\n"; + exit(1); + } + + // Calculate the pin indices. The pin will be located at the grid's reference + // location since that's the only lat/lon location we know about. + + data.x_pin = -(x_values[0] / dx_m); + data.y_pin = -(y_values[0] / dy_m); + + data.dx_km = dx_m / m_per_km; + data.dy_km = dy_m / m_per_km; + data.nx = _xDim->getSize(); + data.ny = _yDim->getSize(); + + data.dump(); + + // Instantiate grid + + grid.set(data); + if (dy_m < 0) grid.set_swap_to_north(true); - mlog << Error << "\n" << method_name << " -> " - << "Lambert azimuthal equal area grid not handled in MET.\n\n"; - exit(1); } @@ -1643,7 +1932,6 @@ void NcCfFile::get_grid_mapping_lambert_conformal_conic(const NcVar *grid_mappin << "Units not given for X coordinate variable -- assuming meters.\n\n"; } else { - //const char *x_coord_units_name = x_coord_units_att->getValues(att->as_string(0); if (0 == x_coord_units_name.length()) { mlog << Warning << "\n" << method_name << " -> " << "Cannot extract X coordinate units from netCDF file -- " @@ -1668,7 +1956,6 @@ void NcCfFile::get_grid_mapping_lambert_conformal_conic(const NcVar *grid_mappin << "Units not given for Y coordinate variable -- assuming meters.\n\n"; } else { - //const char *y_coord_units_name = y_coord_units_att->getValues(att->as_string(0); if (0 == y_coord_units_name.length()) { mlog << Warning << "\n" << method_name << " -> " << "Cannot extract Y coordinate units from netCDF file -- " @@ -1692,13 +1979,11 @@ void NcCfFile::get_grid_mapping_lambert_conformal_conic(const NcVar *grid_mappin long x_counts = GET_NC_SIZE_P(_xDim); double x_values[x_counts]; - //_xCoordVar->get(x_values, &x_counts); get_nc_data(_xCoordVar, x_values); long y_counts = GET_NC_SIZE_P(_yDim); double y_values[y_counts]; - //_yCoordVar->get(y_values, &y_counts); get_nc_data(_yCoordVar, y_values); // Unit conversion @@ -1780,6 +2065,8 @@ void NcCfFile::get_grid_mapping_lambert_conformal_conic(const NcVar *grid_mappin data.ny = _yDim->getSize(); data.so2_angle = 0.0; + data.dump(); + grid.set(data); if (dy_m < 0) grid.set_swap_to_north(true); @@ -2277,6 +2564,8 @@ void NcCfFile::get_grid_mapping_polar_stereographic(const NcVar *grid_mapping_va data.false_east = false_east; data.false_north = false_north; + data.dump(); + grid.set(data); //Note: do not set grid.set_swap_to_north() @@ -2579,6 +2868,8 @@ void NcCfFile::get_grid_mapping_rotated_latitude_longitude(const NcVar *grid_map data.aux_rotation = 0; + data.dump(); + grid.set(data); grid.set_swap_to_north(swap_to_north); @@ -2889,6 +3180,8 @@ void NcCfFile::get_grid_mapping_geostationary( data.scene_id = scene_id_str; } + data.dump(); + // Note: Computing lat/lon was deferred because it took 1 minutes grid.set(data); @@ -3172,6 +3465,8 @@ void NcCfFile::get_grid_from_lat_lon_vars(NcVar *lat_var, NcVar *lon_var, lat_counts, lon_counts, swap_to_north); + data.dump(); + grid.set(data); // resets swap_to_north to false if (swap_to_north) grid.set_swap_to_north(true); } diff --git a/src/libcode/vx_grid/find_grid_by_name.cc b/src/libcode/vx_grid/find_grid_by_name.cc index 00365cb6f3..8197bbec49 100644 --- a/src/libcode/vx_grid/find_grid_by_name.cc +++ b/src/libcode/vx_grid/find_grid_by_name.cc @@ -60,8 +60,7 @@ if ( i.st ) { g.set( *(i.st) ); status = true; } if ( i.ll ) { g.set( *(i.ll) ); status = true; } if ( i.m ) { g.set( *(i.m) ); status = true; } if ( i.g ) { g.set( *(i.g) ); status = true; } -// if ( i.la ) { g.set( *(i.la) ); status = true; } -if ( i.lg ) { g.set( *(i.lg) ); status = true; } +if ( i.la ) { g.set( *(i.la) ); status = true; } return ( status ); @@ -205,22 +204,6 @@ for (j=0; jlat_LL == g2->lat_LL && - g1->lat_UL == g2->lat_UL && - g1->lat_LR == g2->lat_LR && - is_eq (rescale_lon(g1->lon_LL), - rescale_lon(g2->lon_LL), loose_tol) && - is_eq (rescale_lon(g1->lon_UL), - rescale_lon(g2->lon_UL), loose_tol) && - is_eq (rescale_lon(g1->lon_LR), - rescale_lon(g2->lon_LR), loose_tol) && - g1->nx == g2->nx && - g1->ny == g2->ny ) status = true; - -return ( status ); - -} -*/ - -//////////////////////////////////////////////////////////////////////// - - -bool is_eq(const LaeaGrib2Data * g1, const LaeaGrib2Data * g2) +bool is_eq(const LaeaData * g1, const LaeaData * g2) { diff --git a/src/libcode/vx_grid/grid_base.h b/src/libcode/vx_grid/grid_base.h index efc810b85e..81059c0353 100644 --- a/src/libcode/vx_grid/grid_base.h +++ b/src/libcode/vx_grid/grid_base.h @@ -70,8 +70,7 @@ class GridInfo { void set(const GaussianData &); void set(const GoesImagerData &); void set(const TcrmwData &); - // void set(const LaeaData &); - void set(const LaeaGrib2Data &); + void set(const LaeaData &); void set(const SemiLatLonData &); void create_grid(Grid &) const; @@ -88,8 +87,7 @@ class GridInfo { const GaussianData * g; // allocated const GoesImagerData * gi; // allocated const TcrmwData * tc; // allocated - // const LaeaData * la; // allocated - const LaeaGrib2Data * lg; // allocated + const LaeaData * la; // allocated const SemiLatLonData * sl; // allocated }; @@ -209,8 +207,8 @@ class Grid : public GridInterface { Grid(const GaussianData &); Grid(const GoesImagerData &); Grid(const TcrmwData &); - // Grid(const LaeaData &); - Grid(const LaeaGrib2Data &); + Grid(const LaeaData &); + Grid(const LaeaNetcdfData &); Grid(const SemiLatLonData &); virtual ~Grid(); Grid(const Grid &); @@ -231,8 +229,8 @@ class Grid : public GridInterface { void set (const GaussianData &); void set (const GoesImagerData &); void set (const TcrmwData &); - // void set (const LaeaData &); - void set (const LaeaGrib2Data &); + void set (const LaeaData &); + void set (const LaeaNetcdfData &); void set (const SemiLatLonData &); void set_swap_to_north(bool swap_to_north); @@ -293,8 +291,7 @@ extern bool is_eq(const RotatedLatLonData *, const RotatedLatLonData *); extern bool is_eq(const MercatorData *, const MercatorData *); extern bool is_eq(const GaussianData *, const GaussianData *); extern bool is_eq(const GoesImagerData *, const GoesImagerData *); -// extern bool is_eq(const LaeaData *, const LaeaData *); -extern bool is_eq(const LaeaGrib2Data *, const LaeaGrib2Data *); +extern bool is_eq(const LaeaData *, const LaeaData *); extern bool is_eq(const SemiLatLonData *, const SemiLatLonData *); diff --git a/src/libcode/vx_grid/laea_grid.cc b/src/libcode/vx_grid/laea_grid.cc index 914873662b..b9b84c1f48 100644 --- a/src/libcode/vx_grid/laea_grid.cc +++ b/src/libcode/vx_grid/laea_grid.cc @@ -68,143 +68,132 @@ clear(); //////////////////////////////////////////////////////////////////////// -/* + LaeaGrid::LaeaGrid(const LaeaData & data) { clear(); -memset(&Grib2Data, 0, sizeof(Grib2Data)); +memset(&Data, 0, sizeof(Data)); Data = data; +lat_LL = data.lat_first; +lon_LL = data.lon_first; + +lat_pole = data.standard_lat; +lon_pole = data.central_lon; + +Name = data.name; +SpheroidName = data.spheroid_name; +Data.spheroid_name = SpheroidName.c_str(); + Nx = data.nx; Ny = data.ny; - -Name = data.name; -SpheroidName = data.geoid; -if ( strcmp(data.geoid, "WGS_84") == 0 ) { +geoid.set_ab(data.equatorial_radius_km, data.polar_radius_km); - geoid.set_ab(6378.137, 6356.752); +geoid.set_name(data.spheroid_name); - geoid.set_name("WGS_84"); +aff.set_mb(1.0/(data.dx_km), 0.0, 0.0, 1.0/(data.dy_km), 0.0, 0.0); -} else { +double xx, yy; - mlog << Error << "\nLaeaGrid::LaeaGrid(const LaeaData &) -> " - << "unrecognized geoid ... \"" << data.geoid << "\"\n\n"; +latlon_to_xy(data.lat_first, data.lon_first, xx, yy); - exit ( 1 ); +aff.set_translation(-xx, -yy); -} + // + // done + // -calc_aff(); +return; } -*/ + //////////////////////////////////////////////////////////////////////// -LaeaGrid::LaeaGrid(const LaeaGrib2Data & grib2_data) +LaeaGrid::LaeaGrid(const LaeaNetcdfData & nc) { -clear(); - -memset(&Data, 0, sizeof(Data)); -Data = grib2_data; +double u, v, lat, lon; +const double tol = 1.0e-5; -lat_LL = grib2_data.lat_first; -lon_LL = grib2_data.lon_first; +clear(); -lat_pole = grib2_data.standard_lat; -lon_pole = grib2_data.central_lon; +Data.name = nc.name; -Name = grib2_data.name; -SpheroidName = grib2_data.spheroid_name; -Data.spheroid_name = SpheroidName.c_str(); +Data.radius_km = 0.0; -Nx = grib2_data.nx; -Ny = grib2_data.ny; +Data.is_sphere = false; -geoid.set_ab(grib2_data.equatorial_radius_km, grib2_data.polar_radius_km); +Data.equatorial_radius_km = nc.semi_major_axis_km; +Data.polar_radius_km = nc.semi_minor_axis_km; -geoid.set_name(grib2_data.spheroid_name); +Data.dx_km = nc.dx_km; +Data.dy_km = nc.dy_km; -// aff.set_mb(1.0, 0.0, 0.0, 1.0, 1.0, 1.0); -// aff.set_mb(2.0, 0.0, 0.0, 2.0, 1.0, 1.0); +Data.standard_lat = nc.proj_origin_lat; +Data.central_lon = nc.proj_origin_lon; -double s = 1.0; +Data.nx = nc.nx; +Data.ny = nc.ny; -aff.set_mb(s, 0.0, 0.0, s, 1.0, 1.0); +Nx = Data.nx; +Ny = Data.ny; -// aff.set_mb(grib2_data.dx_km, 0.0, 0.0, grib2_data.dy_km, 1.0, 1.0); -double xx, yy; +lat_pole = nc.proj_origin_lat; +lon_pole = nc.proj_origin_lon; -// latlon_to_xy(Data.lat_pole, Data.lon_pole, xx, yy); -// aff.set_pin(xx, yy, 0.5*(Nx - 1.0), 0.5*(Ny - 1.0)); +Name = Data.name; -latlon_to_xy(grib2_data.lat_first, grib2_data.lon_first, xx, yy); -aff.set_pin(xx, yy, 0.0, 0.0); +if ( fabs((nc.semi_major_axis_km - nc.semi_minor_axis_km)/(nc.semi_major_axis_km)) < tol ) { - // - // done - // + Data.radius_km = nc.semi_major_axis_km; -return; + Data.is_sphere = true; } +geoid.set_ab(Data.equatorial_radius_km, Data.polar_radius_km); -//////////////////////////////////////////////////////////////////////// +Data.spheroid_name = "Undefined"; -/* -void LaeaGrid::calc_aff() - -{ +geoid.set_name(Data.spheroid_name); -double u_LL, v_LL; -double u_LR, v_LR; -double u_UL, v_UL; +aff.set_mb(1.0/(Data.dx_km), 0.0, 0.0, 1.0/(Data.dy_km), 0.0, 0.0); -double x_LL, y_LL; -double x_LR, y_LR; -double x_UL, y_UL; +latlon_to_xy(nc.proj_origin_lat, nc.proj_origin_lon, u, v); +aff.set_translation(nc.x_pin - u, nc.y_pin - v); -snyder_latlon_to_xy(Data.lat_LL, Data.lon_LL, u_LL, v_LL); -snyder_latlon_to_xy(Data.lat_LR, Data.lon_LR, u_LR, v_LR); -snyder_latlon_to_xy(Data.lat_UL, Data.lon_UL, u_UL, v_UL); + //////////////////////// +xy_to_latlon(0.0, 0.0, lat, lon); -x_LL = 0.0; -y_LL = 0.0; +Data.lat_first = lat; +Data.lon_first = lon; -x_LR = Data.nx - 1.0; -y_LR = 0.0; +xy_to_latlon(Nx - 1.0, 0.0, lat, lon); -x_UL = 0.0; -y_UL = Data.ny - 1.0; +lat_LR = lat; +lon_LR = lon; +xy_to_latlon(0.0, Ny - 1.0, lat, lon); -aff.set_three_points( +lat_UL = lat; +lon_UL = lon; - u_LL, v_LL, x_LL, y_LL, - - u_LR, v_LR, x_LR, y_LR, - - u_UL, v_UL, x_UL, y_UL - -); - - -return; + // + // done + // } -*/ + //////////////////////////////////////////////////////////////////////// @@ -234,8 +223,7 @@ lon_LR = 0.0; lat_pole = 0.0; -memset(&Data, 0, sizeof(Data)); -// memset(&Grib2Data, 0, sizeof(Grib2Data)); +memset(&Data, 0, sizeof(Data)); return; @@ -251,13 +239,10 @@ void LaeaGrid::latlon_to_xy(double lat, double lon, double & x, double & y) cons double u, v; - snyder_latlon_to_xy(lat, lon, u, v); - uv_to_xy(u, v, x, y); - return; } @@ -289,7 +274,6 @@ lat1 = lat_pole; lambda0 = -lon_pole; - beta1 = geoid.beta(lat1); m1 = geoid.m_func(lat1); @@ -306,7 +290,6 @@ rho = sqrt( uu*uu + vv*vv ); // Eq 24-28, page 189 ce = 2.0*asind(rho/(2.0*Rq)); // Eq 24-29, page 189 - beta = asind( cosd(ce)*sind(beta1) + ((D*v)/rho)*sind(ce)*cosd(beta1) ); // Eq 24-30, page 189 s2 = sind(2.0*beta); @@ -317,25 +300,20 @@ num = u*sind(ce); denom = D*rho*cosd(beta1)*cosd(ce) - D*D*v*sind(beta1)*sind(ce); - -// lon = lambda0 + atand(num/denom); // maybe want atan2 here? // Eq 24-26, page 188 -lon = lambda0 + atan2d(num, denom); // maybe want atan2 here? +lon = lambda0 + atan2d(num, denom); // Eq 24-26, page 188 lon = -lon; - t2 = E2/3.0 + (31.0*E4)/180.0 + (517.0*E6)/5040.0; t4 = (23.0*E4)/360.0 + (251.0*E6)/3780.0; t6 = (761.0*E6)/45360.0; - cor = t2*s2 + t4*s4 + t6*s6; // Eq 3-18, page 189 lat = beta + cor*deg_per_rad; - // // done // @@ -355,13 +333,11 @@ double LaeaGrid::calc_area(int x, int y) const double u[4], v[4]; double sum; - xy_to_uv(x - 0.5, y - 0.5, u[0], v[0]); // lower left xy_to_uv(x + 0.5, y - 0.5, u[1], v[1]); // lower right xy_to_uv(x + 0.5, y + 0.5, u[2], v[2]); // upper right xy_to_uv(x - 0.5, y + 0.5, u[3], v[3]); // upper left - sum = uv_closedpolyline_area(u, v, 4); sum *= earth_radius_km*earth_radius_km; @@ -426,7 +402,7 @@ const char * LaeaGrid::projection_name() const { -return "Polar Laea"; +return "Laea"; } @@ -537,8 +513,6 @@ void LaeaGrid::dump(ostream & out, int depth) const Indent prefix(depth); - - out << prefix << "Name = "; out << prefix << "SpheroidName = "; @@ -572,14 +546,21 @@ ConcatString LaeaGrid::serialize(const char *sep) const { ConcatString a; +char junk[256]; -a << "Projection: Labmbert Azimuthal Equal Area" << sep; +a << "Projection: Lambert Azimuthal Equal Area" << sep; a << "Nx: " << Nx << sep; a << "Ny: " << Ny << sep; a << "SpheroidName: " << SpheroidName << sep; +snprintf(junk, sizeof(junk), "Lat_LL: %.3f", lat_LL); a << junk << sep; +snprintf(junk, sizeof(junk), "Lon_LL: %.3f", lon_LL); a << junk << sep; + +snprintf(junk, sizeof(junk), "Lat_Pole: %.3f", lat_pole); a << junk << sep; +snprintf(junk, sizeof(junk), "Lon_Pole: %.3f", lon_pole); a << junk << sep; + // // done // @@ -665,13 +646,11 @@ const double S = sind(lat); const double E = geoid.e(); const double es = E*S; - z = sqrt(1.0 - es*es); z = C/z; - return z; } @@ -694,56 +673,43 @@ double beta1, beta; double m1, lambda, lambda0, lat1, delta; A = geoid.a_km(); -// A = 1.0; - lambda = -lon; lambda0 = -lon_pole; lat1 = lat_pole; - delta = lambda - lambda0; - beta1 = geoid.beta(lat1); beta = geoid.beta(lat); - Qp = geoid.qp_direct(); Rq = A*sqrt(0.5*Qp); - m1 = snyder_m_func(lat1); - - B = 1.0 + sind(beta1)*sind(beta) + cosd(beta1)*cosd(beta)*cosd(delta); B = sqrt(2.0/B); B = Rq*B; - D = Rq*cosd(beta1); D = (A*m1)/D; - x_snyder = cosd(beta)*sind(delta); x_snyder *= B*D; - - y_snyder = cosd(beta1)*sind(beta) - sind(beta1)*cosd(beta)*cosd(delta); y_snyder *= B/D; - return; } @@ -812,7 +778,7 @@ exit ( 1 ); //////////////////////////////////////////////////////////////////////// -/* + Grid::Grid(const LaeaData & data) { @@ -822,11 +788,11 @@ init_from_scratch(); set(data); } -*/ + //////////////////////////////////////////////////////////////////////// -/* + void Grid::set(const LaeaData & data) { @@ -838,7 +804,7 @@ rep = new LaeaGrid (data); if ( !rep ) { mlog << Error << "\nGrid::set(const LaeaData &) -> " - << "memory allocation error\n\n"; + << "memory allocation error\n\n"; exit ( 1 ); @@ -847,18 +813,18 @@ if ( !rep ) { return; } -*/ + //////////////////////////////////////////////////////////////////////// -Grid::Grid(const LaeaGrib2Data & grib2_data) +Grid::Grid(const LaeaNetcdfData & data) { init_from_scratch(); -set(grib2_data); +set(data); } @@ -866,17 +832,17 @@ set(grib2_data); //////////////////////////////////////////////////////////////////////// -void Grid::set(const LaeaGrib2Data & grib2_data) +void Grid::set(const LaeaNetcdfData & data) { clear(); -rep = new LaeaGrid (grib2_data); +rep = new LaeaGrid (data); if ( !rep ) { - mlog << Error << "\nGrid::set(const LaeaGrib2Data &) -> " + mlog << Error << "\nGrid::set(const LaeaNetcdfData &) -> " << "memory allocation error\n\n"; exit ( 1 ); diff --git a/src/libcode/vx_grid/laea_grid.h b/src/libcode/vx_grid/laea_grid.h index 52347cf632..3a9c92751a 100644 --- a/src/libcode/vx_grid/laea_grid.h +++ b/src/libcode/vx_grid/laea_grid.h @@ -47,37 +47,30 @@ class LaeaGrid : public GridRep { LaeaGrid(); ~LaeaGrid(); - // LaeaGrid(const LaeaData &); - LaeaGrid(const LaeaGrib2Data &); + LaeaGrid(const LaeaData &); + LaeaGrid(const LaeaNetcdfData &); Spheroid geoid; - // LaeaData Data; - // LaeaGrib2Data Grib2Data; - LaeaGrib2Data Data; + LaeaData Data; ///////////////////////////////////// - // from the old Laea_Data class + double lat_LL; // lower left + double lon_LL; + double lat_UL; // upper left + double lon_UL; - double lat_LL; // lower left - double lon_LL; + double lat_LR; // lower right + double lon_LR; - double lat_UL; // upper left - double lon_UL; - - double lat_LR; // lower right - double lon_LR; - - double lat_pole; - double lon_pole; + double lat_pole; + double lon_pole; ///////////////////////////////////// - // bool UseGrib2Data; - double snyder_m_func(double lat) const; void snyder_latlon_to_xy(double lat, double lon, double & x_snyder, double & y_snyder) const; @@ -88,9 +81,6 @@ class LaeaGrid : public GridRep { void clear(); - - void calc_aff(); - void xy_to_uv(double x, double y, double & u, double & v) const; void uv_to_xy(double u, double v, double & x, double & y) const; diff --git a/src/libcode/vx_grid/laea_grid_defs.h b/src/libcode/vx_grid/laea_grid_defs.h index 925d977fa2..5382aa752b 100644 --- a/src/libcode/vx_grid/laea_grid_defs.h +++ b/src/libcode/vx_grid/laea_grid_defs.h @@ -16,23 +16,24 @@ //////////////////////////////////////////////////////////////////////// -/* -struct LaeaData { - const char * name; // not allocated - const char * geoid; // not allocated +struct LaeaNetcdfData { - double lat_LL; // lower left - double lon_LL; + const char * name; - double lat_UL; // upper left - double lon_UL; + double prime_meridian_lon; - double lat_LR; // lower right - double lon_LR; + double semi_major_axis_km; + double semi_minor_axis_km; - double lat_pole; - double lon_pole; + double proj_origin_lat; + double proj_origin_lon; + + double x_pin; + double y_pin; + + double dx_km; + double dy_km; int nx; int ny; @@ -40,7 +41,7 @@ struct LaeaData { void dump() const; }; -*/ + //////////////////////////////////////////////////////////////////////// @@ -56,7 +57,7 @@ struct LaeaData { // -struct LaeaGrib2Data { +struct LaeaData { const char * name; // not allocated const char * spheroid_name; // not allocated @@ -85,72 +86,6 @@ struct LaeaGrib2Data { }; -//////////////////////////////////////////////////////////////////////// - - - // - // Laea grid definitions - // - -/* -static const LaeaData laea_test_data = { - - - "laea_test", // name - - "WGS_84", // geoid - - 31.7462, // lat_LL // lower left - 10.4346, // lon_LL - - 67.0228, // lat_UL // upper left - 39.5358, // lon_UL - - 31.9877, // lat_LR // lower right - -29.421, // lon_LR - - 55.0, // lat_pole - -10.0, // lon_pole - - 1900, // nx - 2200, // ny - -}; -*/ - -static const LaeaGrib2Data laea_grib2_test_data = { - - - "laea_grib2_test", - - "Grib template 4", - - 0.0, // radius_km - - 0.5*6378.1370, // equatorial_radius_km - 0.5*6356.752314, // polar_radius_km - - 44.5172, // lat_first - 17.1171, // lon_first - - 54.900000, // standard_lat - 2.500000, // central_lon - - 2.000000, // dx_km - 2.000000, // dy_km - - 1042, // nx - 970, // ny - - // 970, // ny - // 1042, // nx - - false // is_sphere - - -}; - - //////////////////////////////////////////////////////////////////////// @@ -159,5 +94,3 @@ static const LaeaGrib2Data laea_grib2_test_data = { //////////////////////////////////////////////////////////////////////// - - diff --git a/src/libcode/vx_nc_util/grid_output.cc b/src/libcode/vx_nc_util/grid_output.cc index e0d2e57280..5d5cfbcfb8 100644 --- a/src/libcode/vx_nc_util/grid_output.cc +++ b/src/libcode/vx_nc_util/grid_output.cc @@ -36,8 +36,7 @@ static void rotated_latlon_grid_output (const GridInfo &, NcFile * ncfile); static void stereographic_grid_output (const GridInfo &, NcFile * ncfile); static void mercator_grid_output (const GridInfo &, NcFile * ncfile); static void gaussian_grid_output (const GridInfo &, NcFile * ncfile); -// static void laea_grid_output (const GridInfo &, NcFile * ncfile); -static void laea_grib2_grid_output (const GridInfo &, NcFile * ncfile); +static void laea_grid_output (const GridInfo &, NcFile * ncfile); static void semilatlon_grid_output (const GridInfo &, NcFile * ncfile, NcDim &, NcDim &); static void write_semilatlon_var (NcFile * ncfile, const char *, NcDim *, const NumArray &, const char *, @@ -67,8 +66,7 @@ else if ( info.ll ) latlon_grid_output (info, ncfile); else if ( info.rll ) rotated_latlon_grid_output (info, ncfile); else if ( info.m ) mercator_grid_output (info, ncfile); else if ( info.g ) gaussian_grid_output (info, ncfile); -// else if ( info.la ) laea_grid_output (info, ncfile); -else if ( info.lg ) laea_grib2_grid_output (info, ncfile); +else if ( info.la ) laea_grid_output (info, ncfile); else if ( info.sl ) semilatlon_grid_output (info, ncfile, lat_dim, lon_dim); else { @@ -646,7 +644,7 @@ return; //////////////////////////////////////////////////////////////////////// -/* + void laea_grid_output(const GridInfo & info, NcFile * ncfile) { @@ -657,126 +655,6 @@ const LaeaData & data = *(info.la); ncfile->putAtt("Projection", "Lambert Azimuthal Equal Area"); -if(data.geoid) ncfile->putAtt("geoid", data.geoid); - - // - // lat_LL - // - -snprintf(junk, sizeof(junk), "%f degrees_north", data.lat_LL); - -ncfile->putAtt("lat_LL", junk); - - // - // lon_LL - // - -t = data.lon_LL; - -if ( !west_longitude_positive ) t = -t; - -snprintf(junk, sizeof(junk), "%f degrees_east", t); - -ncfile->putAtt("lon_LL", junk); - - // - // lat_UL - // - -snprintf(junk, sizeof(junk), "%f degrees_north", data.lat_UL); - -ncfile->putAtt("lat_UL", junk); - - // - // lon_UL - // - -t = data.lon_UL; - -if ( !west_longitude_positive ) t = -t; - -snprintf(junk, sizeof(junk), "%f degrees_east", t); - -ncfile->putAtt("lon_UL", junk); - - // - // lat_LR - // - -snprintf(junk, sizeof(junk), "%f degrees_north", data.lat_LR); - -ncfile->putAtt("lat_LR", junk); - - // - // lon_LR - // - -t = data.lon_LR; - -if ( !west_longitude_positive ) t = -t; - -snprintf(junk, sizeof(junk), "%f degrees_east", t); - -ncfile->putAtt("lon_LR", junk); - - // - // lat_pole - // - -snprintf(junk, sizeof(junk), "%f degrees_north", data.lat_pole); - -ncfile->putAtt("lat_pole", junk); - - // - // lon_pole - // - -t = data.lon_pole; - -if ( !west_longitude_positive ) t = -t; - -snprintf(junk, sizeof(junk), "%f degrees_east", t); - -ncfile->putAtt("lon_pole", junk); - - // - // nx - // - -snprintf(junk, sizeof(junk), "%d", data.nx); - -ncfile->putAtt("nx", junk); - - // - // ny - // - -snprintf(junk, sizeof(junk), "%d", data.ny); - -ncfile->putAtt("ny", junk); - - // - // done - // - -return; - -} -*/ - -//////////////////////////////////////////////////////////////////////// - - -void laea_grib2_grid_output(const GridInfo & info, NcFile * ncfile) - -{ - -char junk[256]; -double t; -const LaeaGrib2Data & data = *(info.lg); - -ncfile->putAtt("Projection", "Grib2 Lambert Azimuthal Equal Area"); - ncfile->putAtt("spheroid_name", data.spheroid_name); // diff --git a/src/tools/other/pb2nc/pb2nc.cc b/src/tools/other/pb2nc/pb2nc.cc index 84f9a0fd9c..260062de28 100644 --- a/src/tools/other/pb2nc/pb2nc.cc +++ b/src/tools/other/pb2nc/pb2nc.cc @@ -54,9 +54,10 @@ // 015 02/10/18 Halley Gotway Add message_type_group_map. // 016 07/23/18 Halley Gotway Support masks defined by gen_vx_mask. // 017 07/06/22 Howard Soh METplus-Internal #19 Rename main to met_main -// 018 09/12/22 Prestopnik MET #2227 Remove namespace std and netCDF +// 018 09/12/22 Prestopnik, J. MET #2227 Remove namespace std and netCDF // from header files -// +// 019 07/21/23 Prestopnik, J. MET #2615 Add #include to compile +// successfully using gcc12 //////////////////////////////////////////////////////////////////////// @@ -72,6 +73,7 @@ #include #include +#include #include #include "main.h"