diff --git a/src/extpar_topo_to_buffer.f90 b/src/extpar_topo_to_buffer.f90 index 257fb832..ce724010 100644 --- a/src/extpar_topo_to_buffer.f90 +++ b/src/extpar_topo_to_buffer.f90 @@ -6,7 +6,7 @@ ! V1_0 2010/12/21 Hermann Asensio ! Initial release ! V1_1 2011/01/20 Hermann Asensio -! small bug fixes accroding to Fortran compiler warnings +! small bug fixes according to Fortran compiler warnings ! V1_2 2011/03/25 Hermann Asensio ! update to support ICON refinement grids ! V1_4 2011/04/21 Anne Roches @@ -44,6 +44,8 @@ !> \author Hermann Asensio PROGRAM extpar_topo_to_buffer + USE, INTRINSIC :: iso_c_binding + USE mo_logging USE info_extpar, ONLY: info_print USE mo_kind, ONLY: wp, i4 @@ -59,7 +61,8 @@ PROGRAM extpar_topo_to_buffer USE mo_cosmo_grid, ONLY: COSMO_grid - USE mo_icon_grid_data, ONLY: ICON_grid + USE mo_icon_grid_data, ONLY: ICON_grid, & + & icon_grid_region USE mo_io_units, ONLY: filename_max @@ -136,51 +139,92 @@ PROGRAM extpar_topo_to_buffer IMPLICIT NONE - CHARACTER (len=filename_max) :: netcdf_filename, & - & namelist_grid_def, & - & namelist_topo_data_input, & !< file with input namelist with GLOBE data information - & namelist_scale_sep_data_input, &!< file with input namelist with scale separated data information - & namelist_oro_smooth, & !< file with orography smoothing information (switches) - & namelist_lrad, & !< file with opo information (switches) - & topo_files(1:max_tiles), & !< filenames globe raw data - & sgsl_files(1:max_tiles), & !< filenames subgrid-slope - & orography_buffer_file, & !< name for orography buffer file - & orography_output_file, & !< name for orography output file - & sgsl_output_file, & !< name for sgsl output file - & raw_data_orography_path, & !< path to raw data - & raw_data_scale_sep_orography_path, & !< path to raw data - & scale_sep_files(1:max_tiles) !< filenames globe raw data + CHARACTER (len=filename_max) :: netcdf_filename, & + & namelist_grid_def, & + & namelist_topo_data_input, & !< file with input namelist with GLOBE data information + & namelist_scale_sep_data_input, &!< file with input namelist with scale separated data information + & namelist_oro_smooth, & !< file with orography smoothing information (switches) + & namelist_lrad, & !< file with opo information (switches) + & topo_files(1:max_tiles), & !< filenames globe raw data + & sgsl_files(1:max_tiles), & !< filenames subgrid-slope + & orography_buffer_file, & !< name for orography buffer file + & orography_output_file, & !< name for orography output file + & sgsl_output_file, & !< name for sgsl output file + & raw_data_orography_path, & !< path to raw data + & raw_data_scale_sep_orography_path, & !< path to raw data + & scale_sep_files(1:max_tiles) !< filenames globe raw data - REAL(KIND=wp) :: undefined !< value to indicate undefined grid elements + REAL(KIND=wp) :: undefined !< value to indicate undefined grid elements - INTEGER (KIND=i4) :: & - & k,ie,je,ke, & - & igrid_type, & !< target grid type, 1 for ICON, 2 for COSMO, 3 for GME grid - & ntiles_column, & !< number of tile columns in total domain - & ntiles_row, & !< number of tile rows in total domain - & ilow_pass_oro, & - & numfilt_oro, & - & ifill_valley, & - & ilow_pass_xso, & - & numfilt_xso - - INTEGER (KIND=i4), ALLOCATABLE :: topo_startrow(:), & !< startrow indeces for each GLOBE tile - & topo_endrow(:), & !< endrow indeces for each GLOBE tile - & topo_startcolumn(:), & !< starcolumn indeces for each GLOBE tile - & topo_endcolumn(:) !< endcolumn indeces for each GLOBE tile - - REAL (KIND=wp) :: eps_filter, & - & rfill_valley, & - & rxso_mask - - LOGICAL :: lsso_param, & - & lcompute_sgsl=.FALSE., & !compute subgrid slope - & lpreproc_oro = .FALSE., & !TRUE: preproc raw oro data FALSE: read directly from NetCDF - & lscale_separation=.FALSE., & - & lscale_file= .FALSE., & - & lsubtract_mean_slope = .FALSE., & - & lfilter_oro, & - & lxso_first + INTEGER (KIND=i4) :: & + & k,ie,je,ke, & + & igrid_type, & !< target grid type, 1 for ICON, 2 for COSMO, 3 for GME grid + & ntiles_column, & !< number of tile columns in total domain + & ntiles_row, & !< number of tile rows in total domain + & ilow_pass_oro, & + & numfilt_oro, & + & ifill_valley, & + & ilow_pass_xso, & + & numfilt_xso + + INTEGER (KIND=i4), ALLOCATABLE :: topo_startrow(:), & !< startrow indeces for each GLOBE tile + & topo_endrow(:), & !< endrow indeces for each GLOBE tile + & topo_startcolumn(:), & !< starcolumn indeces for each GLOBE tile + & topo_endcolumn(:) !< endcolumn indeces for each GLOBE tile + + INTEGER(c_int) :: num_cell_c, num_vertex_c, num_hori_c + INTEGER(c_int) :: grid_type_c + REAL (c_double) :: radius_c, ray_org_elev_c + INTEGER(c_int) :: refine_factor_c, itype_scaling_c + CHARACTER(KIND=c_char, len=2000) :: buffer_c + INTEGER(c_int) :: buffer_len_c + REAL(c_double), ALLOCATABLE :: clon_c(:), & + & clat_c(:), & + & hsurf_c(:), & + & vlon_c(:), & + & vlat_c(:), & + & horizon_topo_c(:, :), & + & skyview_topo_c(:) + INTEGER(c_int), ALLOCATABLE :: cells_of_vertex_c(:, :) + + REAL (KIND=wp) :: eps_filter, & + & rfill_valley, & + & rxso_mask + + LOGICAL :: lsso_param, & + & lcompute_sgsl=.FALSE., & !compute subgrid slope + & lpreproc_oro = .FALSE., & !TRUE: preproc raw oro data FALSE: read directly from NetCDF + & lscale_separation=.FALSE., & + & lscale_file= .FALSE., & + & lsubtract_mean_slope, & + & lfilter_oro, & + & lxso_first + + INTERFACE + SUBROUTINE horizon_svf_comp(clon_c, clat_c, hsurf_c, & + & vlon_c, vlat_c, & + & cells_of_vertex_c, & + & horizon_topo_c, skyview_topo_c, & + & num_cell_c, num_vertex_c, num_hori_c, & + & grid_type_c, radius_c, & + & ray_org_elev_c, refine_factor_c, & + & itype_scaling_c, buffer_c, buffer_len_c) bind(C, name="horizon_svf_comp") + USE iso_c_binding + IMPLICIT NONE + REAL(c_double), DIMENSION(*), INTENT(in) :: clon_c, clat_c + REAL(c_double), DIMENSION(*), INTENT(in) :: hsurf_c + REAL(c_double), DIMENSION(*), INTENT(in) :: vlon_c, vlat_c + INTEGER(c_int), DIMENSION(*), INTENT(in) :: cells_of_vertex_c + REAL(c_double), DIMENSION(*), INTENT(inout) :: horizon_topo_c + REAL(c_double), DIMENSION(*), INTENT(inout) :: skyview_topo_c + INTEGER(c_int), value :: num_cell_c, num_vertex_c, num_hori_c + INTEGER(c_int), value :: grid_type_c + REAL(c_double), value :: radius_c, ray_org_elev_c + INTEGER(c_int), value :: refine_factor_c, itype_scaling_c + CHARACTER(KIND=c_char), dimension(*) :: buffer_c + INTEGER(c_int) :: buffer_len_c + END SUBROUTINE horizon_svf_comp + END INTERFACE namelist_grid_def = 'INPUT_grid_org' namelist_scale_sep_data_input = 'INPUT_SCALE_SEP' @@ -518,8 +562,77 @@ PROGRAM extpar_topo_to_buffer CALL compute_lradtopo(nhori,tg,hh_topo,slope_asp_topo,slope_ang_topo, & & horizon_topo,skyview_topo) ELSEIF ( igrid_type == igrid_icon ) THEN - CALL lradtopo_icon(nhori, radius, min_circ_cov,tg, hh_topo, horizon_topo, & - & skyview_topo, max_missing, itype_scaling) + + ! ----------------------------------------------------------------------- + ! Old Fortran implementation + ! ----------------------------------------------------------------------- + + ! CALL lradtopo_icon(nhori, radius, min_circ_cov,tg, hh_topo, horizon_topo, & + ! & skyview_topo, max_missing, itype_scaling) + + ! ----------------------------------------------------------------------- + ! New C++ ray-tracing based implementation + ! ----------------------------------------------------------------------- + + ! Cast non-contiguous arrays to C types + ALLOCATE(clon_c(icon_grid_region%ncells)) + ALLOCATE(clat_c(icon_grid_region%ncells)) + DO k = 1, icon_grid_region%ncells + clon_c(k) = REAL(icon_grid_region%cells%center(k)%lon, KIND=c_double) + clat_c(k) = REAL(icon_grid_region%cells%center(k)%lat, KIND=c_double) + END DO + ALLOCATE(vlon_c(icon_grid_region%nverts)) + ALLOCATE(vlat_c(icon_grid_region%nverts)) + DO k = 1, icon_grid_region%nverts + vlon_c(k) = REAL(icon_grid_region%verts%vertex(k)%lon, KIND=c_double) + vlat_c(k) = REAL(icon_grid_region%verts%vertex(k)%lat, KIND=c_double) + END DO + + ! Cast contiguous array to C type + ALLOCATE(hsurf_c(icon_grid_region%ncells)) + hsurf_c = REAL(hh_topo(:,1,1), KIND=c_double) + + ! Adjust two-dimensional index arary: type and starting index + ALLOCATE(cells_of_vertex_c(icon_grid_region%nverts, 6)) + cells_of_vertex_c = INT(icon_grid_region%verts%cell_index - 1, & + & KIND=c_int) + + ! Cast scalar values to C types + num_cell_c = INT(icon_grid_region%ncells, KIND=c_int) + num_vertex_c = INT(icon_grid_region%nverts, KIND=c_int) + num_hori_c = INT(nhori, KIND=c_int) + radius_c = REAL(radius, KIND=c_double) + itype_scaling_c = INT(itype_scaling, KIND=c_int) + + ! Constant settings + grid_type_c = 1 ! triangle mesh buidling from ICON grid (0 or 1) + ray_org_elev_c = 0.2 ! elevation of ray origin above ground level [m] + refine_factor_c = 10 ! number of sub-sampling within azimuth sector + + ! Allocate output arrays + ALLOCATE(horizon_topo_c(icon_grid_region%ncells, nhori)) + ALLOCATE(skyview_topo_c(icon_grid_region%ncells)) + + buffer_len_c = len(buffer_c) + CALL horizon_svf_comp(clon_c, clat_c, hsurf_c, & + & vlon_c, vlat_c, & + & cells_of_vertex_c, & + & horizon_topo_c, skyview_topo_c, & + & num_cell_c, num_vertex_c, num_hori_c, & + & grid_type_c, radius_c, & + & ray_org_elev_c, refine_factor_c, & + & itype_scaling_c, buffer_c, buffer_len_c) + WRITE(message_text,*) buffer_c(3:buffer_len_c) + CALL logging%info(message_text) + + ! Cast output to Fortran types + horizon_topo(:,1,1,:) = REAL(horizon_topo_c) + skyview_topo(:,1,1) = REAL(skyview_topo_c) + + DEALLOCATE(clon_c, clat_c, hsurf_c, vlon_c, vlat_c) + DEALLOCATE(cells_of_vertex_c) + DEALLOCATE(horizon_topo_c, skyview_topo_c) + ENDIF ENDIF diff --git a/src/mo_lradtopo_horayzon.cpp b/src/mo_lradtopo_horayzon.cpp new file mode 100755 index 00000000..25c51636 --- /dev/null +++ b/src/mo_lradtopo_horayzon.cpp @@ -0,0 +1,786 @@ +// C++ program to compute topographic horizon and sky view factor + +#define _USE_MATH_DEFINES +#include +#include +#include // Fortran interface +#include // Fortran interface +#include +#include +#include +#include +#include +#include +#include +#include + +// Namespace +#if defined(RTC_NAMESPACE_USE) + RTC_NAMESPACE_USE +#endif + +//----------------------------------------------------------------------------- +// Definition of geometries +//----------------------------------------------------------------------------- + +// Point in 3D space +struct geom_point{ + double x, y, z; + // also used for geographic coordinates: lon (x), lat (y), elevation (z) +}; + +// Vector in 3D space +struct geom_vector{ + double x, y, z; +}; + +// Vertex (for Embree) +struct Vertex{ + float x, y, z; +}; + +// Triangle specified by vertex indices (for Embree) +struct Triangle{ + int v0, v1, v2; +}; +// Indices should be 32-bit unsigned integers according to the Embree +// documentation. However, until 2'147'483'647, the binary representation +// between signed/unsigned integers is identical. + +//----------------------------------------------------------------------------- +// Functions (not dependent on Embree) +//----------------------------------------------------------------------------- + +/** + * @brief Converts degree to radian. + * @param ang Input angle [deg]. + * @return Output angle [rad]. + */ +inline double deg2rad(double ang) { + return ((ang / 180.0) * M_PI); +} + +/** + * @brief Converts radian to degree. + * @param ang Input angle [rad]. + * @return Output angle [deg]. + */ +inline double rad2deg(double ang) { + return ((ang / M_PI) * 180.0); +} + +/** + * @brief Computes the dot product between two vectors. + * @param a Vector a. + * @param b Vector b. + * @return Resulting dot product. + */ +inline double dot_product(geom_vector a, geom_vector b) { + return (a.x * b.x + a.y * b.y + a.z * b.z); +} + +/** + * @brief Computes cross dot product between two vectors. + * @param a Vector a. + * @param b Vector b. + * @return Resulting cross product. + */ +inline geom_vector cross_product(geom_vector a, geom_vector b) { + geom_vector c = {a.y * b.z - a.z * b.y, + a.z * b.x - a.x * b.z, + a.x * b.y - a.y * b.x}; + return c; +} + +/** + * @brief Computes the unit vector (normalised vector) of a vector in-place. + * @param a Vector a. + */ +void unit_vector(geom_vector& a) { + double vector_mag = sqrt(a.x * a.x + a.y * a.y + a.z * a.z); + a.x /= vector_mag; + a.y /= vector_mag; + a.z /= vector_mag; +} + +/** + * @brief Rotates vector v around unit vector k with a given angle. + * + * This function rotates vector v around a unit vector k with a given angle + * according to the Rodrigues' rotation formula. For performance reasons, + * trigonometric function have to be pre-computed. + * + * @param v Vector that should be rotated. + * @param k Unit vector specifying the rotation axis. + * @param ang_rot_sin Sine of the rotation angle. + * @param ang_rot_cos Cosine of the rotation angle. + * @return Rotated vector. + */ +inline geom_vector vector_rotation(geom_vector v, geom_vector k, + double ang_rot_sin, double ang_rot_cos) { + geom_vector v_rot; + double term = dot_product(k, v) * (1.0 - ang_rot_cos); + v_rot.x = v.x * ang_rot_cos + (k.y * v.z - k.z * v.y) * ang_rot_sin + + k.x * term; + v_rot.y = v.y * ang_rot_cos + (k.z * v.x - k.x * v.z) * ang_rot_sin + + k.y * term; + v_rot.z = v.z * ang_rot_cos + (k.x * v.y - k.y * v.x) * ang_rot_sin + + k.z * term; + return v_rot; +} + +/** + * @brief Returns indices that would sort an array in ascending order. + * @param values Input values. + * @return Indices that would sort the array. + */ +std::vector sort_index(std::vector& values){ + std::vector index(values.size()); + for (size_t i = 0 ; i < index.size() ; i++) { + index[i] = i; + } + std::sort(index.begin(), index.end(), [&](const int& a, const int& b){ + return (values[a] < values[b]); + }); + return index; +} + +/** + * @brief Transforms geographic to ECEF coordinates in-place. + * + * This function transforms geographic longitude/latitude to earth-centered, + * earth-fixed (ECEF) coordinates. A spherical Earth is assumed. + * + * @param points Points (lon, lat elevation) in geographic coordinates + * [rad, rad, m]. + * @param rad_earth Radius of Earth [m]. + */ +void lonlat2ecef(std::vector& points, double rad_earth){ + for (size_t i = 0; i < points.size(); i++){ + double sin_lon = sin(points[i].x); + double cos_lon = cos(points[i].x); + double sin_lat = sin(points[i].y); + double cos_lat = cos(points[i].y); + double elevation = points[i].z; + points[i].x = (rad_earth + elevation) * cos_lat * cos_lon; + points[i].y = (rad_earth + elevation) * cos_lat * sin_lon; + points[i].z = (rad_earth + elevation) * sin_lat; + } +} + +/** + * @brief Transforms points from ECEF to ENU coordinates in-place. + * @param points Points (x, y, z) in ECEF coordinates [m]. + * @param lon_orig Longitude of ENU coordinate system origin [rad]. + * @param lat_orig Latitude of ENU coordinate system origin [rad]. + * @param rad_earth Radius of Earth [m]. + */ +void ecef2enu_point(std::vector& points, double lon_orig, + double lat_orig, double rad_earth){ + double sin_lon = sin(lon_orig); + double cos_lon = cos(lon_orig); + double sin_lat = sin(lat_orig); + double cos_lat = cos(lat_orig); + double x_ecef_orig = rad_earth * cos(lat_orig) * cos(lon_orig); + double y_ecef_orig = rad_earth * cos(lat_orig) * sin(lon_orig); + double z_ecef_orig = rad_earth * sin(lat_orig); + double x_enu, y_enu, z_enu; + for (size_t i = 0; i < points.size(); i++){ + x_enu = - sin_lon * (points[i].x - x_ecef_orig) + + cos_lon * (points[i].y - y_ecef_orig); + y_enu = - sin_lat * cos_lon * (points[i].x - x_ecef_orig) + - sin_lat * sin_lon * (points[i].y - y_ecef_orig) + + cos_lat * (points[i].z - z_ecef_orig); + z_enu = + cos_lat * cos_lon * (points[i].x - x_ecef_orig) + + cos_lat * sin_lon * (points[i].y - y_ecef_orig) + + sin_lat * (points[i].z - z_ecef_orig); + points[i].x = x_enu; + points[i].y = y_enu; + points[i].z = z_enu; + } +} + +/** + * @brief Builds the triangle mesh from the ICON grid. + * + * This function builds the triangle mesh from the ICON grid cell circumcenters + * (and vertices). Two options are available: + * - 0: Build triangle mesh solely from ICON grid cell circumcenters + * (non-unique triangulation of hexa- and pentagons; relatively long + * triangle edges can cause artefacts in horizon computation) + * - 1: Build triangle mesh from ICON grid cell circumcenters and vertices + * (elevation at vertices is computed as mean from adjacent cell + * circumcenters; triangulation is unique and artefacts are reduced) + * + * @param clon Longitude of ICON grid cell circumcenters [rad]. + * @param clat Latitudes of ICON grid cell circumcenters [rad]. + * @param hsurf Elevation of ICON grid cell circumcenters [m]. + * @param vlon Longitude of ICON grid cell vertices [rad]. + * @param vlat Latitudes of ICON grid cell vertices [rad]. + * @param cells_of_vertex Indices of ICON cells adjacent to ICON vertices. + * @param num_cell Number of ICON grid cells. + * @param num_vertex Number of ICON grid vertices. + * @param grid_type Grid type option for building mesh. + * @param vertices Vertices of build triangle mesh. + * @param vertex_of_triangle Indices of triangles' vertices. + * @return Number of triangles in build mesh. + */ +int build_triangle_mesh(double* clon, double* clat, double* hsurf, + double* vlon, double* vlat, int* cells_of_vertex, + int num_cell, int num_vertex, int grid_type, + std::vector& vertices, std::vector& vertex_of_triangle){ + int ind_cell; + int ind_cov; + for (int i = 0; i < num_cell; i++){ + vertices[i].x = clon[i]; + vertices[i].y = clat[i]; + vertices[i].z = hsurf[i]; + } + if (grid_type == 0) { + std::cout << "Build triangle mesh solely from ICON grid cell" + << " circumcenters\n (non-unique triangulation)" << std::endl; + int ind_1, ind_2; + for (int ind_vertex = 0; ind_vertex < num_vertex; ind_vertex++){ + std::vector angles; + angles.reserve(6); + for (int j = 0; j < 6; j++){ + ind_cell = cells_of_vertex[num_vertex * j + ind_vertex]; + if (ind_cell != -2) { + double angle = atan2(clon[ind_cell] - vlon[ind_vertex], + clat[ind_cell] - vlat[ind_vertex]); + // clockwise angle from positive latitude-axis (y-axis) + if (angle < 0.0) { + angle += 2.0 * M_PI; + } + angles.push_back(angle); + } + } + if (angles.size() >= 3){ + // at least 3 vertices are needed to create one or multiple + // triangles(s) from the polygon + std::vector ind_sort = sort_index(angles); + ind_1 = 1; + ind_2 = 2; + for (size_t j = 0; j < (angles.size() - 2); j++){ + ind_cov = num_vertex * ind_sort[0] + ind_vertex; + vertex_of_triangle.push_back(cells_of_vertex[ind_cov]); + ind_cov = num_vertex * ind_sort[ind_1] + ind_vertex; + vertex_of_triangle.push_back(cells_of_vertex[ind_cov]); + ind_1 ++; + ind_cov = num_vertex * ind_sort[ind_2] + ind_vertex; + vertex_of_triangle.push_back(cells_of_vertex[ind_cov]); + ind_2 ++; + // add indices of triangle's vertices in clockwise order + } + } + } + } else { + std::cout << "Build triangle mesh from ICON grid cell circumcenters" + << " and vertices\n (unique triangulation)" << std::endl; + int ind_add = num_cell; + int ind[7] = {0, 1, 2, 3, 4, 5, 0}; + for (int ind_vertex = 0; ind_vertex < num_vertex; ind_vertex++){ + std::vector angles; + angles.reserve(6); + double hsurf_mean = 0.0; + for (int j = 0; j < 6; j++){ + ind_cell = cells_of_vertex[num_vertex * j + ind_vertex]; + if (ind_cell != -2) { + double angle = atan2(clon[ind_cell] - vlon[ind_vertex], + clat[ind_cell] - vlat[ind_vertex]); + // clockwise angle from positive latitude-axis (y-axis) + if (angle < 0.0) { + angle += 2.0 * M_PI; + } + angles.push_back(angle); + hsurf_mean += hsurf[ind_cell]; + } + } + if (angles.size() == 6){ + vertices.push_back({vlon[ind_vertex], vlat[ind_vertex], + hsurf_mean / 6.0}); + std::vector ind_sort = sort_index(angles); + for (int j = 0; j < 6; j++){ + ind_cov = num_vertex * ind_sort[ind[j]] + ind_vertex; + vertex_of_triangle.push_back(cells_of_vertex[ind_cov]); + ind_cov = num_vertex * ind_sort[ind[j + 1]] + ind_vertex; + vertex_of_triangle.push_back(cells_of_vertex[ind_cov]); + vertex_of_triangle.push_back(ind_add); + } + ind_add += 1; + } + } + } + int num_triangle = vertex_of_triangle.size() / 3; + std::cout << "Number of triangles in mesh: " << num_triangle << std::endl; + return num_triangle; +} + +/** + * @brief Computes the sky view factor for a horizontally aligned plane. + * + * This function computes the sky view factor (SVF) for a horizontally aligned + * plane. Three methods are available: + * - Visible sky fraction / pure geometric sky view factor + * - Sky view factor / geometric scaled with sin(horizon) + * - Sky view factor additionally scaled with sin(horizon) / + * geometric scaled with sin(horizon)**2 + * + * @param horizon_cell Horizon array [rad]. + * @param horizon_cell_len Length of the horizon array. + * @return Sky view factor [-]. + */ +double (*function_pointer)(double* horizon_cell, int horizon_cell_len); +double pure_geometric_svf(double* horizon_cell, int horizon_cell_len){ + double svf = 0.0; + for(int i = 0; i < horizon_cell_len; i++){ + svf += (1.0 - sin(horizon_cell[i])); + } + svf /= (double)horizon_cell_len; + return svf; +} +double geometric_svf_scaled_1(double* horizon_cell, int horizon_cell_len){ + double svf = 0.0; + for(int i = 0; i < horizon_cell_len; i++){ + svf += (1.0 - (sin(horizon_cell[i]) * sin(horizon_cell[i]))); + } + svf /= (double)horizon_cell_len; + return svf; +} +double geometric_svf_scaled_2(double* horizon_cell, int horizon_cell_len){ + double svf = 0.0; + for(int i = 0; i < horizon_cell_len; i++){ + svf += (1.0 - (sin(horizon_cell[i]) * sin(horizon_cell[i]) + * sin(horizon_cell[i]))); + } + svf /= (double)horizon_cell_len; + return svf; +} + +//----------------------------------------------------------------------------- +// Functions (Embree related) +//----------------------------------------------------------------------------- + +/** + * @brief Error function for device initialiser. + * @param userPtr + * @param error + * @param str + */ +void errorFunction(void* userPtr, enum RTCError error, const char* str) { + printf("error %d: %s\n", error, str); +} + +/** + * @brief Initialises device and registers error handler + * @return Device instance. + */ +RTCDevice initializeDevice() { + RTCDevice device = rtcNewDevice(NULL); + if (!device) { + printf("error %d: cannot create device\n", rtcGetDeviceError(NULL)); + } + rtcSetDeviceErrorFunction(device, errorFunction, NULL); + return device; +} + +/** + * @brief Initialises the Embree scene. + * @param device Initialised device. + * @param vertex_of_triangle Indices of the triangle vertices. + * @param num_triangle Number of triangles. + * @param vertices Vertices of the triangles [m]. + * @return Embree scene. + */ +RTCScene initializeScene(RTCDevice device, int* vertex_of_triangle, + int num_triangle, std::vector& vertices){ + + RTCScene scene = rtcNewScene(device); + rtcSetSceneFlags(scene, RTC_SCENE_FLAG_ROBUST); + RTCGeometry geom = rtcNewGeometry(device, RTC_GEOMETRY_TYPE_TRIANGLE); + + // Vertices + Vertex* vertices_embree = (Vertex*) rtcSetNewGeometryBuffer(geom, + RTC_BUFFER_TYPE_VERTEX, 0, RTC_FORMAT_FLOAT3, sizeof(Vertex), + vertices.size()); + for (size_t i = 0; i < vertices.size(); i++) { + vertices_embree[i].x = (float)vertices[i].x; + vertices_embree[i].y = (float)vertices[i].y; + vertices_embree[i].z = (float)vertices[i].z; + } + + // Cell (triangle) indices to vertices + rtcSetSharedGeometryBuffer(geom, RTC_BUFFER_TYPE_INDEX, 0, + RTC_FORMAT_UINT3, vertex_of_triangle, 0, 3*sizeof(int), num_triangle); + + auto start = std::chrono::high_resolution_clock::now(); + + // Commit geometry and scene + rtcCommitGeometry(geom); + rtcAttachGeometry(scene, geom); + rtcReleaseGeometry(geom); + rtcCommitScene(scene); + + auto end = std::chrono::high_resolution_clock::now(); + std::chrono::duration time = end - start; + std::cout << std::setprecision(2) << std::fixed; + std::cout << "Building bounding volume hierarchy (BVH): " << time.count() + << " s" << std::endl; + + return scene; + +} + +/** + * @brief Ray casting with occlusion testing (hit / no hit). + * @param scene Embree scene. + * @param ox x-coordinate of the ray origin [m]. + * @param oy y-coordinate of the ray origin [m]. + * @param oz z-coordinate of the ray origin [m]. + * @param dx x-component of the ray direction [m]. + * @param dy y-component of the ray direction [m]. + * @param dz z-component of the ray direction [m]. + * @param dist_search Search distance for potential collision [m]. + * @return Collision status (true: hit, false: no hit). + */ +bool castRay_occluded1(RTCScene scene, float ox, float oy, float oz, float dx, + float dy, float dz, float dist_search){ + struct RTCRay ray; + ray.org_x = ox; + ray.org_y = oy; + ray.org_z = oz; + ray.dir_x = dx; + ray.dir_y = dy; + ray.dir_z = dz; + ray.tnear = 0.0; + ray.tfar = dist_search; + ray.mask = 1; + rtcOccluded1(scene, &ray); // intersect ray with scene + return (ray.tfar < 0.0); +} + +/** + * @brief Computes the terrain horizon for a specific point. + * + * This function computes the terrain horizon for a specific point on the + * triangle mesh. It iteratively samples a certain azimuth direction with rays + * until the horizon is found. For all but the first azimuth direction, the + * elevation angle for the search is initialised with a value equal to the + * horizon from the previous azimuth direction +/- the horizon accuracy value. + * + * @param ray_org_x x-coordinate of the ray origin [m]. + * @param ray_org_y y-coordinate of the ray origin [m]. + * @param ray_org_z z-coordinate of the ray origin [m]. + * @param hori_acc Horizon accuracy [rad]. + * @param dist_search Search distance for potential collision [m]. + * @param elev_ang_thresh Threshold angle for sampling in negative elevation + * angle direction [rad]. + * @param scene Embree scene. + * @param num_rays Number of rays casted. + * @param horizon_cell Horizon array [rad]. + * @param horizon_cell_len Length of the horizon array. + * @param azim_shift Azimuth shift for the first azimuth sector [rad]. + * @param sphere_normal Sphere normal at the point location [m]. + * @param north_direction North direction at the point location [m]. + * @param azim_sin Sine of the azimuth angle spacing. + * @param azim_cos Cosine of the azimuth angle spacing. + * @param elev_sin_2ha Sine of the double elevation angle spacing. + * @param elev_cos_2ha Cosine of the double elevation angle spacing. + */ +void terrain_horizon(float ray_org_x, float ray_org_y, float ray_org_z, + double hori_acc, float dist_search, double elev_ang_thresh, + RTCScene scene, size_t &num_rays, + double* horizon_cell, int horizon_cell_len, + double azim_shift, + geom_vector sphere_normal, geom_vector north_direction, + double azim_sin, double azim_cos, + double elev_sin_2ha, double elev_cos_2ha){ + + // Initial ray direction + geom_vector ray_dir; + ray_dir.x = north_direction.x; + ray_dir.y = north_direction.y; + ray_dir.z = north_direction.z; + + // Shift azimuth angle in case of 'refine_factor' > 1 so that first + // azimuth sector is centred around 0.0 deg (North) + ray_dir = vector_rotation(ray_dir, sphere_normal, sin(-azim_shift), + cos(-azim_shift)); + + // Sample along azimuth + double elev_ang = 0.0; + for (int i = 0; i < horizon_cell_len; i++){ + + // Rotation axis + geom_vector rot_axis = cross_product(ray_dir, sphere_normal); + unit_vector(rot_axis); + // not necessarily a unit vector because vectors are mostly not + // perpendicular + + // Find terrain horizon by iterative ray sampling + bool hit = castRay_occluded1(scene, ray_org_x, ray_org_y, + ray_org_z, (float)ray_dir.x, (float)ray_dir.y, (float)ray_dir.z, + dist_search); + num_rays += 1; + if (hit) { // terrain hit -> increase elevation angle + while (hit){ + elev_ang += (2.0 * hori_acc); + ray_dir = vector_rotation(ray_dir, rot_axis, elev_sin_2ha, + elev_cos_2ha); + hit = castRay_occluded1(scene, ray_org_x, ray_org_y, + ray_org_z, (float)ray_dir.x, (float)ray_dir.y, + (float)ray_dir.z, dist_search); + num_rays += 1; + } + horizon_cell[i] = elev_ang - hori_acc; + } else { // terrain not hit -> decrease elevation angle + while ((!hit) && (elev_ang > elev_ang_thresh)){ + elev_ang -= (2.0 * hori_acc); + ray_dir = vector_rotation(ray_dir, rot_axis, -elev_sin_2ha, + elev_cos_2ha); // sin(-x) == -sin(x), cos(x) == cos(-x) + hit = castRay_occluded1(scene, ray_org_x, ray_org_y, + ray_org_z, (float)ray_dir.x, (float)ray_dir.y, + (float)ray_dir.z, dist_search); + num_rays += 1; + } + horizon_cell[i] = elev_ang + hori_acc; + } + + // Azimuthal rotation of ray direction (clockwise; first to east) + ray_dir = vector_rotation(ray_dir, sphere_normal, -azim_sin, + azim_cos); // sin(-x) == -sin(x), cos(x) == cos(-x) + + } + +} + +//----------------------------------------------------------------------------- +// Main function +//----------------------------------------------------------------------------- + +extern "C" { // Fortran interface +void horizon_svf_comp(double* clon, double* clat, double* hsurf, + double* vlon, double* vlat, + int* cells_of_vertex, + double* horizon, double* skyview, + int num_cell, int num_vertex, int azim_num, + int grid_type, double dist_search_dp, + double ray_org_elev, int refine_factor, + // int svf_type){ // Fortran interface + int svf_type, char* buffer, int* buffer_len){ // Fortran interface + + // Redirect std::cout (Fortran interface) + std::stringstream string_stream; + auto* old_buf = std::cout.rdbuf(string_stream.rdbuf()); + + // Fixed settings + double hori_acc = deg2rad(0.25); // horizon accuracy [rad] + double elev_ang_thresh = deg2rad(-85.0); + // threshold for sampling in negative elevation angle direction [rad] + // - relevant for 'void sampling directions' at edge of mesh + // - necessary requirement: (elev_ang_thresh - (2.0 * hori_acc)) > -90.0 + + // Constants + double rad_earth = 6371229.0; // ICON/COSMO earth radius [m] + + // Type casting + float dist_search = (float)dist_search_dp; + + std::cout << "------------------------------------------------------------" + << "-------------------" << std::endl; + std::cout << "Horizon and SVF computation with Intel Embree (v1.1)" + << std::endl; + std::cout << "------------------------------------------------------------" + << "-------------------" << std::endl; + + // Build triangle mesh from ICON grid + auto start_mesh = std::chrono::high_resolution_clock::now(); + std::vector vertices(num_cell); + std::vector vertex_of_triangle; + int num_triangle = build_triangle_mesh(clon, clat, hsurf, + vlon, vlat, cells_of_vertex, + num_cell, num_vertex, grid_type, + vertices, vertex_of_triangle); + auto end_mesh = std::chrono::high_resolution_clock::now(); + std::chrono::duration time_mesh = end_mesh - start_mesh; + std::cout << std::setprecision(2) << std::fixed; + std::cout << "Triangle mesh building: " << time_mesh.count() << " s" + << std::endl; + + // In-place transformation from geographic to ECEF coordinates + lonlat2ecef(vertices, rad_earth); + + // Earth center and North Pole in ECEF coordinates + std::vector earth_centre(1); + earth_centre[0].x = 0.0; + earth_centre[0].y = 0.0; + earth_centre[0].z = 0.0; + std::vector north_pole(1); + north_pole[0].x = 0.0; + north_pole[0].y = 0.0; + north_pole[0].z = rad_earth; + + // Origin of ENU coordinate system + double x_orig = 0.0; + double y_orig = 0.0; + double z_orig = 0.0; + for (int i = 0; i < num_cell; i++){ + x_orig += vertices[i].x; + y_orig += vertices[i].y; + z_orig += vertices[i].z; + } + double radius = sqrt(x_orig * x_orig + y_orig * y_orig + z_orig * z_orig); + double lon_orig = atan2(y_orig, x_orig); + double lat_orig = asin(z_orig / radius); + // works correctly for ICON domains containing the North/South Pole and/or + // crossing the +/- 180 deg meridian + + // In-place transformation from ECEF to ENU coordinates + std::cout << std::setprecision(4) << std::fixed; + std::cout << "Origin of ENU coordinate system: " << rad2deg(lat_orig) + << " deg lat, " << rad2deg(lon_orig) << " deg lon" << std::endl; + ecef2enu_point(vertices, lon_orig, lat_orig, rad_earth); + ecef2enu_point(earth_centre, lon_orig, lat_orig, rad_earth); + ecef2enu_point(north_pole, lon_orig, lat_orig, rad_earth); + + // Build bounding volume hierarchy (BVH) + RTCDevice device = initializeDevice(); + RTCScene scene = initializeScene(device, vertex_of_triangle.data(), + num_triangle, vertices); + + // Evaluated trigonometric functions for rotation along azimuth/elevation + // angle + int horizon_cell_len = azim_num * refine_factor; + double azim_sin = sin(deg2rad(360.0) / (double)horizon_cell_len); + double azim_cos = cos(deg2rad(360.0) / (double)horizon_cell_len); + double elev_sin_2ha = sin(2.0 * hori_acc); + double elev_cos_2ha = cos(2.0 * hori_acc); + // Note: sin(-x) == -sin(x), cos(x) == cos(-x) + + // Compute shift for azimuth angle so that first azimuth sector is + // centred around 0.0 deg (North) in case of 'refine_factor' > 1 + double azim_shift; + if (refine_factor == 1) { + azim_shift = 0.0; + } else { + azim_shift = -(deg2rad(360.0) / (2.0 * azim_num)) + + (deg2rad(360.0) / (2.0 * (double)horizon_cell_len)); + } + + // Select algorithm for sky view factor computation + std::cout << "Sky View Factor computation algorithm: "; + if (svf_type == 0) { + std::cout << "pure geometric SVF" << std::endl; + function_pointer = pure_geometric_svf; + } else if (svf_type == 1) { + std::cout << "geometric scaled with sin(horizon)" << std::endl; + function_pointer = geometric_svf_scaled_1; + } else if (svf_type == 2) { + std::cout << "geometric scaled with sin(horizon)**2" << std::endl; + function_pointer = geometric_svf_scaled_2; + } + + auto start_ray = std::chrono::high_resolution_clock::now(); + size_t num_rays = 0; + + num_rays += tbb::parallel_reduce( + tbb::blocked_range(0, num_cell), 0.0, + [&](tbb::blocked_range r, size_t num_rays) { // parallel + + // for (size_t i = 0; i < (size_t)num_cell; i++){ // serial + for (size_t i=r.begin(); i()); // parallel + + auto end_ray = std::chrono::high_resolution_clock::now(); + std::chrono::duration time_ray = end_ray - start_ray; + std::cout << std::setprecision(2) << std::fixed; + std::cout << "Ray tracing: " << time_ray.count() << " s" << std::endl; + + // Print number of rays needed for location and azimuth direction + std::cout << "Number of rays shot: " << num_rays << std::endl; + double ratio = (double)num_rays / (double)(num_cell * azim_num); + std::cout << std::setprecision(2) << std::fixed; + std::cout << "Average number of rays per cell and azimuth sector: " + << ratio << std::endl; + + std::cout << "------------------------------------------------------------" + << "-------------------" << std::endl; + + // Restore original std::cout and copy output to buffer (Fortran interface) + std::cout.rdbuf(old_buf); + std::string output = string_stream.str(); + int copy_len = std::min(*buffer_len, (int)output.size()); + std::memcpy(buffer, output.c_str(), copy_len); + *buffer_len = copy_len; + +} +} // Fortran interface