From 6e9709c5826c82e830a6c8ee9e5132f4893d13a3 Mon Sep 17 00:00:00 2001 From: Anthony Islas Date: Sat, 29 Mar 2025 14:34:31 -0700 Subject: [PATCH 01/20] First draft, open geo nc file and index file for topo --- util/CMakeLists.txt | 29 ++- util/src/CMakeLists.txt | 31 ++- util/src/compute_gwdo.F | 420 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 478 insertions(+), 2 deletions(-) create mode 100644 util/src/compute_gwdo.F diff --git a/util/CMakeLists.txt b/util/CMakeLists.txt index 684ecae3..ba2cbbae 100644 --- a/util/CMakeLists.txt +++ b/util/CMakeLists.txt @@ -11,6 +11,7 @@ add_executable( calc_ecmwf_p ) add_executable( mod_levs ) if ( DEFINED WRF_DIR ) add_executable( int2nc ) + add_executable( compute_gwdo ) endif() add_executable( height_ukmo ) @@ -56,6 +57,27 @@ if ( DEFINED WRF_DIR ) PRIVATE ${netCDF_INCLUDE_DIRS} ) + target_link_libraries( + compute_gwdo + PRIVATE + WRF::WRF_Core + netCDF::netcdff + ) + set_target_properties( + compute_gwdo + PROPERTIES + # Just dump everything in here + Fortran_MODULE_DIRECTORY ${CMAKE_INSTALL_PREFIX}/modules/util/compute_gwdo + Fortran_FORMAT FREE + ) + install( + TARGETS compute_gwdo + EXPORT ${EXPORT_NAME}Targets + RUNTIME DESTINATION bin/ + ARCHIVE DESTINATION lib/ + LIBRARY DESTINATION lib/ + ) + endif() # Every single target @@ -66,7 +88,7 @@ foreach( TARGET ${ALL_UTIL_TARGETS} util_common ) ${TARGET} PROPERTIES # Just dump everything in here - Fortran_MODULE_DIRECTORY ${CMAKE_INSTALL_PREFIX}/modules/util/ + Fortran_MODULE_DIRECTORY ${CMAKE_INSTALL_PREFIX}/modules/util/${TARGET}/ Fortran_FORMAT FREE ) @@ -77,6 +99,11 @@ foreach( TARGET ${ALL_UTIL_TARGETS} util_common ) PRIVATE util_common ) + target_include_directories( + ${TARGET} + PRIVATE + $ + ) endif() # Add these to the export targets diff --git a/util/src/CMakeLists.txt b/util/src/CMakeLists.txt index d357c8d9..f162f0e5 100644 --- a/util/src/CMakeLists.txt +++ b/util/src/CMakeLists.txt @@ -69,4 +69,33 @@ target_sources( height_ukmo PRIVATE height_ukmo.F - ) \ No newline at end of file + ) + +target_sources( + compute_gwdo + PRIVATE + compute_gwdo.F + # module_input_module.F90 + # ${PROJECT_SOURCE_DIR}/geogrid/src/cio.c + # ${PROJECT_SOURCE_DIR}/geogrid/src/wrf_debug.F + # # ${PROJECT_SOURCE_DIR}/geogrid/src/bitarray_module.F + # ${PROJECT_SOURCE_DIR}/geogrid/src/constants_module.F + # # ${PROJECT_SOURCE_DIR}/geogrid/src/module_stringutil.F + # # ${PROJECT_SOURCE_DIR}/geogrid/src/geogrid.F + # ${PROJECT_SOURCE_DIR}/geogrid/src/gridinfo_module.F + # # ${PROJECT_SOURCE_DIR}/geogrid/src/hash_module.F + # # ${PROJECT_SOURCE_DIR}/geogrid/src/interp_module.F + # # ${PROJECT_SOURCE_DIR}/geogrid/src/list_module.F + # # ${PROJECT_SOURCE_DIR}/geogrid/src/llxy_module.F + # ${PROJECT_SOURCE_DIR}/geogrid/src/misc_definitions_module.F + # ${PROJECT_SOURCE_DIR}/geogrid/src/module_debug.F + # # ${PROJECT_SOURCE_DIR}/geogrid/src/module_map_utils.F + # # ${PROJECT_SOURCE_DIR}/geogrid/src/output_module.F + # ${PROJECT_SOURCE_DIR}/geogrid/src/parallel_module.F + # # ${PROJECT_SOURCE_DIR}/geogrid/src/process_tile_module.F + # # ${PROJECT_SOURCE_DIR}/geogrid/src/proc_point_module.F + # ${PROJECT_SOURCE_DIR}/geogrid/src/queue_module.F + # # ${PROJECT_SOURCE_DIR}/geogrid/src/read_geogrid.c + # # ${PROJECT_SOURCE_DIR}/geogrid/src/smooth_module.F + # # ${PROJECT_SOURCE_DIR}/geogrid/src/source_data_module.F + ) diff --git a/util/src/compute_gwdo.F b/util/src/compute_gwdo.F new file mode 100644 index 00000000..9903d6e0 --- /dev/null +++ b/util/src/compute_gwdo.F @@ -0,0 +1,420 @@ +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! MODULE GRIDINFO_MODULE +! +! This module handles (i.e., acquires, stores, and makes available) all data +! describing the model domains to be processed. +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +module gridinfo_module + + ! Parameters + integer, parameter :: MAX_DOMAINS = 21 + + ! Variables + integer :: iproj_type, n_domains, io_form_output, dyn_opt + integer, dimension(MAX_DOMAINS) :: parent_grid_ratio, parent_id, ixdim, jydim + integer, dimension(MAX_DOMAINS) :: subgrid_ratio_x, subgrid_ratio_y + real :: known_lat, known_lon, pole_lat, pole_lon, stand_lon, truelat1, truelat2, & + known_x, known_y, dxkm, dykm, phi, lambda, ref_lat, ref_lon, ref_x, ref_y, & + dlatdeg, dlondeg + real, dimension(MAX_DOMAINS) :: parent_ll_x, parent_ll_y, parent_ur_x, parent_ur_y + character (len=256) :: geog_data_path, opt_output_from_geogrid_path, opt_geogrid_tbl_path + + character (len=128), dimension(MAX_DOMAINS) :: geog_data_res + character (len=1) :: gridtype + logical :: do_tiled_output + logical, dimension(MAX_DOMAINS) :: grid_is_active + integer :: debug_level + +contains + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Name: get_grid_params + ! + ! Purpose: This subroutine retrieves all parameters regarding the model domains + ! to be processed by geogrid.exe. This includes map parameters, domain + ! size and location, and nest information. + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine get_grid_params() + + implicit none + + ! Local variables + integer :: i, j, max_dom, funit, io_form_geogrid, interval_seconds + real :: dx, dy, rparent_gridpts + integer, dimension(MAX_DOMAINS) :: i_parent_start, j_parent_start, & + s_we, e_we, s_sn, e_sn, & + start_year, start_month, start_day, start_hour, start_minute, start_second, & + end_year, end_month, end_day, end_hour, end_minute, end_second + character (len=128) :: map_proj + character (len=128), dimension(MAX_DOMAINS) :: start_date, end_date + character (len=3) :: wrf_core + logical :: is_used, nest_outside + logical, dimension(MAX_DOMAINS) :: active_grid + logical :: nocolons + + namelist /share/ wrf_core, max_dom, start_date, end_date, & + start_year, end_year, start_month, end_month, & + start_day, end_day, start_hour, end_hour, & + start_minute, end_minute, start_second, end_second, & + interval_seconds, & + io_form_geogrid, opt_output_from_geogrid_path, & + debug_level, active_grid, & + subgrid_ratio_x, subgrid_ratio_y, & + nocolons + namelist /geogrid/ parent_id, parent_grid_ratio, & + i_parent_start, j_parent_start, s_we, e_we, s_sn, e_sn, & + map_proj, ref_x, ref_y, ref_lat, ref_lon, & + pole_lat, pole_lon, truelat1, truelat2, stand_lon, & + dx, dy, geog_data_res, geog_data_path, opt_geogrid_tbl_path + + ! Set defaults for namelist variables + debug_level = 0 + io_form_geogrid = 2 + wrf_core = 'ARW' + max_dom = 1 + geog_data_path = 'NOT_SPECIFIED' + ref_x = 1.e20 + ref_y = 1.e20 + ref_lat = 1.e20 + ref_lon = 1.e20 + dx = 1.e20 + dy = 1.e20 + map_proj = 'Lambert' + pole_lat = 90.0 + pole_lon = 0.0 + truelat1 = 1.e20 + truelat2 = 1.e20 + stand_lon = 1.e20 + do i=1,MAX_DOMAINS + geog_data_res(i) = 'default' + parent_id(i) = 1 + parent_grid_ratio(i) = 1e9 + s_we(i) = 1 + e_we(i) = 1e9 + s_sn(i) = 1 + e_sn(i) = 1e9 + start_year(i) = 0 + start_month(i) = 0 + start_day(i) = 0 + start_hour(i) = 0 + start_minute(i) = 0 + start_second(i) = 0 + end_year(i) = 0 + end_month(i) = 0 + end_day(i) = 0 + end_hour(i) = 0 + end_minute(i) = 0 + end_second(i) = 0 + start_date(i) = '0000-00-00_00:00:00' + end_date(i) = '0000-00-00_00:00:00' + active_grid(i) = .true. + subgrid_ratio_x(i) = 1 + subgrid_ratio_y(i) = 1 + end do + opt_output_from_geogrid_path = './' + opt_geogrid_tbl_path = 'geogrid/' + interval_seconds = 1e9 + nocolons = .false. + + ! Read parameters from Fortran namelist + do funit=10,100 + inquire(unit=funit, opened=is_used) + if (.not. is_used) exit + end do + open(funit,file='namelist.wps',status='old',form='formatted',err=1000) + read(funit,share) + read(funit,geogrid) + close(funit) + + n_domains = max_dom + + return + + 1000 write( *, * ) 'Error opening file namelist.wps' + + end subroutine get_grid_params + +end module gridinfo_module + +module nc_helpers +contains + subroutine nf90_check_error( stat ) + use netcdf + implicit none + integer, intent( in ) :: stat + if ( stat /= nf90_noerr ) then + ! We had an error + write( 0, * ) __FILE__, ":", __LINE__, " : ", trim( nf90_strerror( stat ) ) + stop 1 + end if + end subroutine nf90_check_error +end module nc_helpers + +module source_data_module + type :: source_index_raw + character (len=32) :: projection = "required" + character (len=32) :: source_type = "required" + logical :: signed = .false. + character (len=32) :: units = "required" + character (len=64) :: description = "required" + real :: dx = -1.0 !< required + real :: dy = -1.0 !< required + real :: known_x = 1.0 + real :: known_y = 1.0 + real :: known_lat = -1.0 !< required + real :: known_lon = -1.0 !< required + real :: stdlon = -2.0 !< not required + real :: truelat1 = -2.0 !< not required + real :: truelat2 = -2.0 !< not required + integer :: wordsize = -1 !< required + integer :: tile_x = -1 !< required + integer :: tile_y = -1 !< required + integer :: tile_z = -2 !< not required + integer :: tile_z_start = -2 !< not required + integer :: tile_z_end = -2 !< not required + integer :: category_min = -2 !< not required + integer :: category_max = -2 !< not required + integer :: tile_bdr = 0 + real :: missing_value = -2.0 !< not required + real :: scale_factor = 1.0 + character (len=32) :: row_order = "bottom_top" + character (len=32) :: endian = "big" + integer :: iswater = 16 + integer :: islake = -1 !< i.e. no separate inland water category (does not mean required) + integer :: isice = 24 + integer :: isurban = 1 + integer :: isoilwater = 14 + character (len=32) :: mminlu = "USGS" + integer :: filename_digits = 5 + end type source_index_raw + +contains + + subroutine read_index_file( source_path, index_data ) + use, intrinsic :: iso_fortran_env + + implicit none + + character (len=*), intent( in ) :: source_path + character (len=256) :: index_path, line_buffer, key, var + logical :: is_used, kv_exist + integer :: iostatus, funit + + type( source_index_raw ), intent( inout ) :: index_data + + write( index_path, "(a)" ) trim(source_path) // 'index' + write( *, * ) "Reading 'index' file at ", index_path + + ! + ! Open the index file for the data source for this field, and read in metadata specs + ! + do funit=10,100 + inquire(unit=funit, opened=is_used) + if (.not. is_used) exit + end do + open(funit,file=trim(index_path),form='formatted',status='old',iostat=iostatus) + + ! TODO add optional to control when not found + if (iostatus /= 0 ) then + write( *, * ) 'Could not open ', trim(index_path) + end if + + do + read( funit, '(a)', iostat=iostatus ) line_buffer + if ( iostatus /= 0 ) then + if ( iostatus == iostat_end ) exit + end if + + ! Process line + ! write( *, * ) "Line : ", trim(line_buffer) + kv_exist = get_kv_pair( line_buffer, key, var ) + if ( kv_exist ) then + ! write( *, * ) "Found kv ", key, " = ", var + if ( key == "projection" ) index_data % projection = var + if ( key == "source_type" ) index_data % source_type = var + if ( key == "signed" .and. var == "yes" ) index_data % signed = .true. + if ( key == "units" ) index_data % units = var + if ( key == "description" ) index_data % description = var + if ( key == "dx" ) read( var, * ) index_data % dx + if ( key == "dy" ) read( var, * ) index_data % dy + if ( key == "known_x" ) read( var, * ) index_data % known_x + if ( key == "known_y" ) read( var, * ) index_data % known_y + if ( key == "known_lat" ) read( var, * ) index_data % known_lat + if ( key == "known_lon" ) read( var, * ) index_data % known_lon + if ( key == "stdlon" ) read( var, * ) index_data % stdlon + if ( key == "truelat1" ) read( var, * ) index_data % truelat1 + if ( key == "truelat2" ) read( var, * ) index_data % truelat2 + if ( key == "wordsize" ) read( var, * ) index_data % wordsize + if ( key == "tile_x" ) read( var, * ) index_data % tile_x + if ( key == "tile_y" ) read( var, * ) index_data % tile_y + if ( key == "tile_z" ) read( var, * ) index_data % tile_z + if ( key == "tile_z_start" ) read( var, * ) index_data % tile_z_start + if ( key == "tile_z_end" ) read( var, * ) index_data % tile_z_end + if ( key == "category_min" ) read( var, * ) index_data % category_min + if ( key == "category_max" ) read( var, * ) index_data % category_max + if ( key == "tile_bdr" ) read( var, * ) index_data % tile_bdr + if ( key == "missing_value" ) read( var, * ) index_data % missing_value + if ( key == "scale_factor" ) read( var, * ) index_data % scale_factor + if ( key == "row_order" ) index_data % row_order = var + if ( key == "endian" ) index_data % endian = var + if ( key == "iswater" ) read( var, * ) index_data % iswater + if ( key == "islake" ) read( var, * ) index_data % islake + if ( key == "isice" ) read( var, * ) index_data % isice + if ( key == "isurban" ) read( var, * ) index_data % isurban + if ( key == "isoilwater" ) read( var, * ) index_data % isoilwater + if ( key == "mminlu" ) index_data % mminlu = var + if ( key == "filename_digits" ) read( var, * ) index_data % filename_digits + end if + end do + end subroutine + + logical function get_kv_pair( line, key, var ) result( found ) + implicit none + character (len=*), intent( in ) :: line + character (len=*), intent( inout ) :: key, var + integer :: idx, idx_comment + ! logical, intent( out ) :: found + + found = .false. + + idx = index( line, "=" ) + idx_comment = index( line, "#" ) + + if ( idx == 0 ) return + + if ( idx_comment /= 0 .and. idx_comment < idx ) then + ! We have a comment before our = and that is not good + write( *, * ) "Key-value parse error, comment before = " + return + end if + + if ( idx_comment == 0 ) idx_comment = len( line ) + + ! We are fine to parse + key = trim(line(1:idx-1)) + var = trim(line(idx+1:idx_comment)) + found = .true. + + end function get_kv_pair + +end module source_data_module + +program compute_gwdo + use netcdf + use nc_helpers + use source_data_module + use gridinfo_module + + implicit none + character (len=128), parameter :: geo_em_prefix = "geo_em.d" + character (len=128) :: geo_em_file, tmp + integer :: i, j, k, l + integer :: ierr, ncid, ndims, dim1, dim2, dim3 + integer, dimension(nf90_max_var_dims) :: dimids + integer :: start_i, end_i, start_j, end_j, start_k, end_k + integer :: mapfacxid, mapfacyid + integer :: latvarid, lonvarid + integer :: convarid, varvarid + integer, dimension(4) :: olvarid, oavarid + + ! I will probably need more + real, dimension(:,:,:), allocatable :: con_data + real, dimension(:,:,:), allocatable :: lats, lons + + real :: cur_lat, cur_lon + type( source_index_raw ) :: index_data + + + + + call get_grid_params() + + write( *, * ) "Found geogrid : ", n_domains + write( *, * ) "Using data from : ", geog_data_path + + write( tmp, "(a)" ) trim(geog_data_path) // "topo_gmted2010_30s/" + call read_index_file( tmp, index_data ) + + write( *, * ) "Source data dx / dy : ", index_data % dx, ", ", index_data % dy + + do i = 1, n_domains + write( geo_em_file, "(A8,I0.2,A3)" ) trim( geo_em_prefix ), i, ".nc" + write( *, * ) "Processing ", geo_em_file + ierr = nf90_open( geo_em_file, NF90_WRITE, ncid ) + call nf90_check_error( ierr ) + + ierr = nf90_inq_varid( ncid, "XLAT_M", latvarid ) + call nf90_check_error( ierr ) + + ierr = nf90_inq_varid( ncid, "XLONG_M", lonvarid ) + call nf90_check_error( ierr ) + + ierr = nf90_inq_varid( ncid, "MAPFAC_MX", mapfacxid ) + call nf90_check_error( ierr ) + + ierr = nf90_inq_varid( ncid, "MAPFAC_MY", mapfacyid ) + call nf90_check_error( ierr ) + + ierr = nf90_inquire_variable( ncid, lonvarid, ndims = ndims, dimids = dimids ) + call nf90_check_error( ierr ) + + write( *, * ) "Variable CLONG has ", ndims, " dims" + + ierr = nf90_inq_varid( ncid, "CON", convarid ) + call nf90_check_error( ierr ) + + ierr = nf90_inq_varid( ncid, "VAR", varvarid ) + call nf90_check_error( ierr ) + + do j = 1, 4 + write( tmp, "(A2,I1)" ) "OL", j + ierr = nf90_inq_varid( ncid, tmp, olvarid(j) ) + write( tmp, "(A2,I1)" ) "OA", j + ierr = nf90_inq_varid( ncid, tmp, oavarid(j) ) + end do + + ierr = nf90_inquire_dimension( ncid, dimids(1), len = dim1 ) + ierr = nf90_inquire_dimension( ncid, dimids(2), len = dim2 ) + ierr = nf90_inquire_dimension( ncid, dimids(3), len = dim3 ) + + write( *, * ) "Dimensions sizes are ", dim1, ", ", dim2, ", ", dim3 + + allocate( con_data( dim1, dim2, dim3 ) ) + allocate( lats( dim1, dim2, dim3 ) ) + allocate( lons( dim1, dim2, dim3 ) ) + + ierr = nf90_get_var( ncid, latvarid, lats ) + ierr = nf90_get_var( ncid, lonvarid, lons ) + + ! Just stuff for now + do j = 1, dim1 + do k = 1, dim2 + do l = 1, dim3 + con_data(j,k,l) = 3.1415926535 + end do + end do + end do + + ! Just stuff for now + do j = 1, dim1 + do k = 1, dim2 + do l = 1, dim3 + cur_lat = lats(j,k,l) + cur_lon = lons(j,k,l) + ! write( *, * ) "Current LL : ", cur_lat, ", ", cur_lon + end do + end do + end do + + ierr = nf90_put_var( ncid, convarid, con_data ) + call nf90_check_error( ierr ) + + + ierr = nf90_close( ncid ) + call nf90_check_error( ierr ) + + end do + + return +end program compute_gwdo + From 6ffb45c6d3ac75c0b8b1843d0c4413a9abe98066 Mon Sep 17 00:00:00 2001 From: Anthony Islas Date: Mon, 31 Mar 2025 12:30:18 -0700 Subject: [PATCH 02/20] Getting source file location for particular latlon in data --- util/src/CMakeLists.txt | 41 +++---- util/src/compute_gwdo.F | 246 ++++++++++++++++------------------------ 2 files changed, 121 insertions(+), 166 deletions(-) diff --git a/util/src/CMakeLists.txt b/util/src/CMakeLists.txt index f162f0e5..691aeff8 100644 --- a/util/src/CMakeLists.txt +++ b/util/src/CMakeLists.txt @@ -76,26 +76,27 @@ target_sources( PRIVATE compute_gwdo.F # module_input_module.F90 - # ${PROJECT_SOURCE_DIR}/geogrid/src/cio.c - # ${PROJECT_SOURCE_DIR}/geogrid/src/wrf_debug.F - # # ${PROJECT_SOURCE_DIR}/geogrid/src/bitarray_module.F - # ${PROJECT_SOURCE_DIR}/geogrid/src/constants_module.F - # # ${PROJECT_SOURCE_DIR}/geogrid/src/module_stringutil.F - # # ${PROJECT_SOURCE_DIR}/geogrid/src/geogrid.F - # ${PROJECT_SOURCE_DIR}/geogrid/src/gridinfo_module.F - # # ${PROJECT_SOURCE_DIR}/geogrid/src/hash_module.F - # # ${PROJECT_SOURCE_DIR}/geogrid/src/interp_module.F - # # ${PROJECT_SOURCE_DIR}/geogrid/src/list_module.F - # # ${PROJECT_SOURCE_DIR}/geogrid/src/llxy_module.F - # ${PROJECT_SOURCE_DIR}/geogrid/src/misc_definitions_module.F - # ${PROJECT_SOURCE_DIR}/geogrid/src/module_debug.F - # # ${PROJECT_SOURCE_DIR}/geogrid/src/module_map_utils.F - # # ${PROJECT_SOURCE_DIR}/geogrid/src/output_module.F + ${PROJECT_SOURCE_DIR}/geogrid/src/cio.c + ${PROJECT_SOURCE_DIR}/geogrid/src/wrf_debug.F + # ${PROJECT_SOURCE_DIR}/geogrid/src/bitarray_module.F + ${PROJECT_SOURCE_DIR}/geogrid/src/constants_module.F + # ${PROJECT_SOURCE_DIR}/geogrid/src/module_stringutil.F + # ${PROJECT_SOURCE_DIR}/geogrid/src/geogrid.F + ${PROJECT_SOURCE_DIR}/geogrid/src/gridinfo_module.F + # ${PROJECT_SOURCE_DIR}/geogrid/src/hash_module.F + # ${PROJECT_SOURCE_DIR}/geogrid/src/interp_module.F + ${PROJECT_SOURCE_DIR}/geogrid/src/list_module.F + ${PROJECT_SOURCE_DIR}/geogrid/src/llxy_module.F + ${PROJECT_SOURCE_DIR}/geogrid/src/misc_definitions_module.F + ${PROJECT_SOURCE_DIR}/geogrid/src/module_debug.F + ${PROJECT_SOURCE_DIR}/geogrid/src/module_map_utils.F + # ${PROJECT_SOURCE_DIR}/geogrid/src/module_map_utils.F + # ${PROJECT_SOURCE_DIR}/geogrid/src/output_module.F # ${PROJECT_SOURCE_DIR}/geogrid/src/parallel_module.F - # # ${PROJECT_SOURCE_DIR}/geogrid/src/process_tile_module.F - # # ${PROJECT_SOURCE_DIR}/geogrid/src/proc_point_module.F + # ${PROJECT_SOURCE_DIR}/geogrid/src/process_tile_module.F + # ${PROJECT_SOURCE_DIR}/geogrid/src/proc_point_module.F # ${PROJECT_SOURCE_DIR}/geogrid/src/queue_module.F - # # ${PROJECT_SOURCE_DIR}/geogrid/src/read_geogrid.c - # # ${PROJECT_SOURCE_DIR}/geogrid/src/smooth_module.F - # # ${PROJECT_SOURCE_DIR}/geogrid/src/source_data_module.F + # ${PROJECT_SOURCE_DIR}/geogrid/src/read_geogrid.c + # ${PROJECT_SOURCE_DIR}/geogrid/src/smooth_module.F + # ${PROJECT_SOURCE_DIR}/geogrid/src/source_data_module.F ) diff --git a/util/src/compute_gwdo.F b/util/src/compute_gwdo.F index 9903d6e0..7de2bd24 100644 --- a/util/src/compute_gwdo.F +++ b/util/src/compute_gwdo.F @@ -1,141 +1,3 @@ -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! MODULE GRIDINFO_MODULE -! -! This module handles (i.e., acquires, stores, and makes available) all data -! describing the model domains to be processed. -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -module gridinfo_module - - ! Parameters - integer, parameter :: MAX_DOMAINS = 21 - - ! Variables - integer :: iproj_type, n_domains, io_form_output, dyn_opt - integer, dimension(MAX_DOMAINS) :: parent_grid_ratio, parent_id, ixdim, jydim - integer, dimension(MAX_DOMAINS) :: subgrid_ratio_x, subgrid_ratio_y - real :: known_lat, known_lon, pole_lat, pole_lon, stand_lon, truelat1, truelat2, & - known_x, known_y, dxkm, dykm, phi, lambda, ref_lat, ref_lon, ref_x, ref_y, & - dlatdeg, dlondeg - real, dimension(MAX_DOMAINS) :: parent_ll_x, parent_ll_y, parent_ur_x, parent_ur_y - character (len=256) :: geog_data_path, opt_output_from_geogrid_path, opt_geogrid_tbl_path - - character (len=128), dimension(MAX_DOMAINS) :: geog_data_res - character (len=1) :: gridtype - logical :: do_tiled_output - logical, dimension(MAX_DOMAINS) :: grid_is_active - integer :: debug_level - -contains - - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! Name: get_grid_params - ! - ! Purpose: This subroutine retrieves all parameters regarding the model domains - ! to be processed by geogrid.exe. This includes map parameters, domain - ! size and location, and nest information. - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine get_grid_params() - - implicit none - - ! Local variables - integer :: i, j, max_dom, funit, io_form_geogrid, interval_seconds - real :: dx, dy, rparent_gridpts - integer, dimension(MAX_DOMAINS) :: i_parent_start, j_parent_start, & - s_we, e_we, s_sn, e_sn, & - start_year, start_month, start_day, start_hour, start_minute, start_second, & - end_year, end_month, end_day, end_hour, end_minute, end_second - character (len=128) :: map_proj - character (len=128), dimension(MAX_DOMAINS) :: start_date, end_date - character (len=3) :: wrf_core - logical :: is_used, nest_outside - logical, dimension(MAX_DOMAINS) :: active_grid - logical :: nocolons - - namelist /share/ wrf_core, max_dom, start_date, end_date, & - start_year, end_year, start_month, end_month, & - start_day, end_day, start_hour, end_hour, & - start_minute, end_minute, start_second, end_second, & - interval_seconds, & - io_form_geogrid, opt_output_from_geogrid_path, & - debug_level, active_grid, & - subgrid_ratio_x, subgrid_ratio_y, & - nocolons - namelist /geogrid/ parent_id, parent_grid_ratio, & - i_parent_start, j_parent_start, s_we, e_we, s_sn, e_sn, & - map_proj, ref_x, ref_y, ref_lat, ref_lon, & - pole_lat, pole_lon, truelat1, truelat2, stand_lon, & - dx, dy, geog_data_res, geog_data_path, opt_geogrid_tbl_path - - ! Set defaults for namelist variables - debug_level = 0 - io_form_geogrid = 2 - wrf_core = 'ARW' - max_dom = 1 - geog_data_path = 'NOT_SPECIFIED' - ref_x = 1.e20 - ref_y = 1.e20 - ref_lat = 1.e20 - ref_lon = 1.e20 - dx = 1.e20 - dy = 1.e20 - map_proj = 'Lambert' - pole_lat = 90.0 - pole_lon = 0.0 - truelat1 = 1.e20 - truelat2 = 1.e20 - stand_lon = 1.e20 - do i=1,MAX_DOMAINS - geog_data_res(i) = 'default' - parent_id(i) = 1 - parent_grid_ratio(i) = 1e9 - s_we(i) = 1 - e_we(i) = 1e9 - s_sn(i) = 1 - e_sn(i) = 1e9 - start_year(i) = 0 - start_month(i) = 0 - start_day(i) = 0 - start_hour(i) = 0 - start_minute(i) = 0 - start_second(i) = 0 - end_year(i) = 0 - end_month(i) = 0 - end_day(i) = 0 - end_hour(i) = 0 - end_minute(i) = 0 - end_second(i) = 0 - start_date(i) = '0000-00-00_00:00:00' - end_date(i) = '0000-00-00_00:00:00' - active_grid(i) = .true. - subgrid_ratio_x(i) = 1 - subgrid_ratio_y(i) = 1 - end do - opt_output_from_geogrid_path = './' - opt_geogrid_tbl_path = 'geogrid/' - interval_seconds = 1e9 - nocolons = .false. - - ! Read parameters from Fortran namelist - do funit=10,100 - inquire(unit=funit, opened=is_used) - if (.not. is_used) exit - end do - open(funit,file='namelist.wps',status='old',form='formatted',err=1000) - read(funit,share) - read(funit,geogrid) - close(funit) - - n_domains = max_dom - - return - - 1000 write( *, * ) 'Error opening file namelist.wps' - - end subroutine get_grid_params - -end module gridinfo_module - module nc_helpers contains subroutine nf90_check_error( stat ) @@ -152,6 +14,7 @@ end module nc_helpers module source_data_module type :: source_index_raw + character (len=256):: source_path = "required" character (len=32) :: projection = "required" character (len=32) :: source_type = "required" logical :: signed = .false. @@ -202,8 +65,9 @@ subroutine read_index_file( source_path, index_data ) type( source_index_raw ), intent( inout ) :: index_data + write( index_data % source_path, "(a)" ) trim(source_path) write( index_path, "(a)" ) trim(source_path) // 'index' - write( *, * ) "Reading 'index' file at ", index_path + write( *, * ) "Reading 'index' file at ", index_data % source_path ! ! Open the index file for the data source for this field, and read in metadata specs @@ -216,7 +80,7 @@ subroutine read_index_file( source_path, index_data ) ! TODO add optional to control when not found if (iostatus /= 0 ) then - write( *, * ) 'Could not open ', trim(index_path) + write( *, * ) 'Could not open ', trim(index_data % source_path) end if do @@ -225,11 +89,10 @@ subroutine read_index_file( source_path, index_data ) if ( iostatus == iostat_end ) exit end if - ! Process line - ! write( *, * ) "Line : ", trim(line_buffer) + ! Process line iff a key-value pair was found kv_exist = get_kv_pair( line_buffer, key, var ) if ( kv_exist ) then - ! write( *, * ) "Found kv ", key, " = ", var + ! consider multiline else if even if not neater if ( key == "projection" ) index_data % projection = var if ( key == "source_type" ) index_data % source_type = var if ( key == "signed" .and. var == "yes" ) index_data % signed = .true. @@ -291,12 +154,86 @@ logical function get_kv_pair( line, key, var ) result( found ) if ( idx_comment == 0 ) idx_comment = len( line ) ! We are fine to parse - key = trim(line(1:idx-1)) - var = trim(line(idx+1:idx_comment)) + key = trim(adjustl(line(1:idx-1))) + var = trim(adjustl(line(idx+1:idx_comment))) found = .true. end function get_kv_pair + integer function get_proj_code( index_data ) result( proj_code ) + use misc_definitions_module + implicit none + type( source_index_raw ), intent( in ) :: index_data + if ( index_data % projection == "lambert" ) then; proj_code = PROJ_LC + else if ( index_data % projection == "polar_wgs84" ) then; proj_code = PROJ_PS_WGS84 + else if ( index_data % projection == "albers_nad83" ) then; proj_code = PROJ_ALBERS_NAD83 + else if ( index_data % projection == "polar" ) then; proj_code = PROJ_PS + else if ( index_data % projection == "mercator" ) then; proj_code = PROJ_MERC + else if ( index_data % projection == "regular_ll" ) then; proj_code = PROJ_LATLON + else; write( *, * ) "ERROR: Unknown projection '", index_data % projection, "'" + end if + end function get_proj_code + + + ! Taken from the souce_data_module of geogrid, but simplified to take stateful index data + subroutine get_tile_fname(index_data, lat, lon, tile_fname, istatus) + + use llxy_module + use gridinfo_module + + implicit none + + type( source_index_raw ), intent( in ) :: index_data + real, intent(in) :: lat, lon + character (len=*), intent(out) :: tile_fname + integer, intent(out) :: istatus + + ! Local variables + integer :: current_domain, idx + real :: rx, ry + + istatus = 0 + write(tile_fname, '(a)') 'TILE.OUTSIDE.DOMAIN' + + ! do idx=1,num_entries + ! if ((index(source_fieldname(idx),trim(fieldname)) /= 0) .and. & + ! (len_trim(source_fieldname(idx)) == len_trim(fieldname))) then + ! if (ilevel == source_priority(idx)) then + ! istatus = 0 + ! exit + ! end if + ! end if + ! end do + + ! if (istatus /= 0) return + + current_domain = iget_selected_domain() + call select_domain(SOURCE_PROJ) + call lltoxy(lat, lon, rx, ry, M) + call select_domain(current_domain) + + ! rx = real(source_tile_x(idx)) * real(floor((rx-0.5*source_dx(idx))/ real(source_tile_x(idx)))) + 1.0 + ! ry = real(source_tile_y(idx)) * real(floor((ry-0.5*source_dy(idx))/ real(source_tile_y(idx)))) + 1.0 + rx = real( index_data % tile_x ) * real(floor((rx-0.5) / real( index_data % tile_x ))) + 1.0 + ry = real( index_data % tile_y ) * real(floor((ry-0.5) / real( index_data % tile_y ))) + 1.0 + + if (rx > 0. .and. ry > 0.) then + select case(index_data % filename_digits) + case (5) + write(tile_fname, '(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim( index_data % source_path ), & + nint(rx),'-',nint(rx)+index_data % tile_x -1,'.',nint(ry),'-',nint(ry)+index_data % tile_y -1 + case (6) + write(tile_fname, '(a,i6.6,a1,i6.6,a1,i6.6,a1,i6.6)') trim( index_data % source_path ), & + nint(rx),'-',nint(rx)+index_data % tile_x -1,'.',nint(ry),'-',nint(ry)+index_data % tile_y -1 + case default + call mprintf(.true., ERROR, 'In GEOGRID.TBL, for index file of data at '// & + 'entry %i, filename_digits must be either 5 or 6.', i1=idx) + istatus = 1 + return + end select + end if + + end subroutine get_tile_fname end module source_data_module program compute_gwdo @@ -304,6 +241,7 @@ program compute_gwdo use nc_helpers use source_data_module use gridinfo_module + use llxy_module implicit none character (len=128), parameter :: geo_em_prefix = "geo_em.d" @@ -335,6 +273,16 @@ program compute_gwdo write( tmp, "(a)" ) trim(geog_data_path) // "topo_gmted2010_30s/" call read_index_file( tmp, index_data ) + ! For now set the source data projection here + call push_source_projection( get_proj_code( index_data ), index_data % stdlon, & + index_data % truelat1, index_data % truelat2, & + index_data % dx, index_data % dy, & + index_data % dy, index_data % dx, & + index_data % known_x, index_data % known_y, & + index_data % known_lat, index_data % known_lon ) + + + write( *, * ) "Source data dx / dy : ", index_data % dx, ", ", index_data % dy do i = 1, n_domains @@ -386,6 +334,9 @@ program compute_gwdo ierr = nf90_get_var( ncid, latvarid, lats ) ierr = nf90_get_var( ncid, lonvarid, lons ) + ! Find correct file and location within file to grab + + ! Just stuff for now do j = 1, dim1 do k = 1, dim2 @@ -401,7 +352,10 @@ program compute_gwdo do l = 1, dim3 cur_lat = lats(j,k,l) cur_lon = lons(j,k,l) - ! write( *, * ) "Current LL : ", cur_lat, ", ", cur_lon + call get_tile_fname( index_data, cur_lat, cur_lon, tmp, ierr ) + if ( j == 10 .and. k == 10 ) then + write( *, * ) "Filename for [", i, ", ", j, ", ", k, "] at LL [", cur_lat, ", ", cur_lon, "] :", tmp + end if end do end do end do From 6aebd97177b7fe661d4255e01d6103793f6996bb Mon Sep 17 00:00:00 2001 From: Anthony Islas Date: Tue, 15 Apr 2025 12:56:03 -0700 Subject: [PATCH 03/20] Source data reading of geogrid format --- util/SourceData.py | 254 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 254 insertions(+) create mode 100644 util/SourceData.py diff --git a/util/SourceData.py b/util/SourceData.py new file mode 100644 index 00000000..90d803ce --- /dev/null +++ b/util/SourceData.py @@ -0,0 +1,254 @@ +import numpy as np +import cartopy + +def get_kv_pair( line, delimiter="=", comment="#" ): + non_comment = line.split( comment, maxsplit=1 )[0] + + if delimiter in non_comment: + key, _, value = non_comment.partition( delimiter ) + return key.strip(), value.strip() + else: + return None, None + +class IndexData: + + def __init__( self, path ): + self.source_path_ = path + # Default values of index file from : + # https://www2.mmm.ucar.edu/wrf/users/wrf_users_guide/build/html/wps.html#index-options + self.projection_ = "required" + self.source_type_ = "required" + self.signed_ = False + self.units_ = "required" + self.description_ = "required" + self.dx_ = -1.0 # required + self.dy_ = -1.0 # required + self.known_x_ = 1.0 + self.known_y_ = 1.0 + self.known_lat_ = -1.0 # required + self.known_lon_ = -1.0 # required + self.stdlon_ = -2.0 # not required + self.truelat1_ = -2.0 # not required + self.truelat2_ = -2.0 # not required + self.wordsize_ = -1 # required + self.tile_x_ = -1 # required + self.tile_y_ = -1 # required + self.tile_z_ = -2 # not required + self.tile_z_start_ = -2 # not required + self.tile_z_end_ = -2 # not required + self.category_min_ = -2 # not required + self.category_max_ = -2 # not required + self.tile_bdr_ = 0 + self.missing_value_ = -2.0 # not required + self.scale_factor_ = 1.0 + self.row_order_ = "bottom_top" + self.endian_ = "big" + self.iswater_ = 16 + self.islake_ = -1 # i.e. no separate inland water category (does not mean required) + self.isice_ = 24 + self.isurban_ = 1 + self.isoilwater_ = 14 + self.mminlu_ = "USGS" + self.filename_digits_ = 5 + + self.read_file() + + def read_file( self ): + with open( self.source_path_ + "/index" ) as f: + for line in f: + key, value = get_kv_pair( line ) + if key is not None: + if key == "projection" : self.projection_ = value + elif key == "source_type" : self.source_type_ = value + elif key == "signed" : self.signed_ = ( value == "yes" ) + elif key == "units" : self.units_ = value + elif key == "description" : self.description_ = value + elif key == "dx" : self.dx_ = float( value ) + elif key == "dy" : self.dy_ = float( value ) + elif key == "known_x" : self.known_x_ = float( value ) + elif key == "known_y" : self.known_y_ = float( value ) + elif key == "known_lat" : self.known_lat_ = float( value ) + elif key == "known_lon" : self.known_lon_ = float( value ) + elif key == "stdlon" : self.stdlon_ = float( value ) + elif key == "truelat1" : self.truelat1_ = float( value ) + elif key == "truelat2" : self.truelat2_ = float( value ) + elif key == "wordsize" : self.wordsize_ = int( value ) + elif key == "tile_x" : self.tile_x_ = int( value ) + elif key == "tile_y" : self.tile_y_ = int( value ) + elif key == "tile_z" : self.tile_z_ = int( value ) + elif key == "tile_z_start" : self.tile_z_start_ = int( value ) + elif key == "tile_z_end" : self.tile_z_end_ = int( value ) + elif key == "category_min" : self.category_min_ = int( value ) + elif key == "category_max" : self.category_max_ = int( value ) + elif key == "tile_bdr" : self.tile_bdr_ = int( value ) + elif key == "missing_value" : self.missing_value_ = float( value ) + elif key == "scale_factor" : self.scale_factor_ = float( value ) + elif key == "row_order" : self.row_order_ = value + elif key == "endian" : self.endian_ = value + elif key == "iswater" : self.iswater_ = int( value ) + elif key == "islake" : self.islake_ = int( value ) + elif key == "isice" : self.isice_ = int( value ) + elif key == "isurban" : self.isurban_ = int( value ) + elif key == "isoilwater" : self.isoilwater_ = int( value ) + elif key == "mminlu" : self.mminlu_ = value + elif key == "filename_digits" : self.filename_digits_ = int( value ) + + def get_dtype( self ): + endian = ">" if self.endian_ == "big" else "<" + signed = "i" if self.signed_ else "u" + return f"{endian}{signed}{self.wordsize_}" + + +class SourceData: + + def __init__( self, name, path ): + self.earth_radius_ = 6371229.0 + self.name_ = name + self.source_path_ = path + self.index_ = IndexData( self.source_path_ ) + + self.projection_ = None + self.npts_x_ = int( 360.0 / self.index_.dx_ ) + self.pts_per_deg_ = int( 1.0 / self.index_.dx_ ) + self.subgrid_m_dx_ = 2.0 * np.pi * self.earth_radius_ / self.npts_x_ + + self.data_ = None + + + # These are all probably wrong + if self.index_.projection_ == "lambert": + self.projection_ = cartopy.crs.LambertConformal( + central_longitude =self.index_.known_lon_, + central_latitude =self.index_.known_lat_, + standard_parallels =( self.index_.truelat1_, self.index_.truelat2_ ) + ) + elif self.index_.projection_ == "polar_wgs84" or self.index_.projection_ == "polar": + self.projection_ = cartopy.crs.Stereographic( + central_longitude =self.index_.known_lon_, + central_latitude =self.index_.known_lat_, + true_scale_latitude =self.index_.truelat1_ + ) + elif self.index_.projection_ == "albers_nad83": + self.projection_ = cartopy.crs.AlbersEqualArea( + central_longitude =self.index_.known_lon_, + central_latitude =self.index_.known_lat_, + standard_parallels =( self.index_.truelat1_, self.index_.truelat2_ ) + ) + elif self.index_.projection_ == "mercator": + self.projection_ = cartopy.crs.Mercator( + central_longitude =self.index_.known_lon_, + # central_latitude =self.index_.known_lat_, + latitude_true_scale =self.index_.truelat1_ + ) + elif self.index_.projection_ == "regular_ll": + self.projection_ = cartopy.crs.LambertCylindrical( central_longitude =self.index_.known_lon_ ) + + + def read_geogrid( self, file ): + rawdata = np.fromfile( + file, + dtype=self.index_.get_dtype() + ) + data = rawdata.astype( np.float32, casting="unsafe" ) * self.index_.scale_factor_ + z_dim = 1 + if self.index_.tile_z_start_ > 0 and self.index_.tile_z_end_ > 0: + z_dim = self.index_.tile_z_end_ - self.index_.tile_z_start_ + elif self.index_.tile_z_ > 0: + z_dim = self.index_.tile_z_ + + data = data.reshape( + ( + z_dim, + self.index_.tile_y_ + 2 * self.index_.tile_bdr_, + self.index_.tile_x_ + 2 * self.index_.tile_bdr_ + ) + ) + if self.index_.row_order_ != "top_bottom": + data = np.flip( data, axis=1 ) + + return data + + def latlon_to_ij( self, lat, lon ): + # Ignore staggering for now + i = 0.0 + j = 0.0 + if self.index_.projection_ == "regular_ll": + delta_lat = lat - self.index_.known_lat_ + delta_lon = lon - self.index_.known_lon_ + + i = ( delta_lon / self.index_.dx_ + self.index_.known_x_ ) % int( 360.0 / self.index_.dx_ ) + j = ( delta_lat / self.index_.dy_ + self.index_.known_y_ ) + + return i, j + + def ij_to_latlon( self, i, j ): + lat = 0.0 + lon = 0.0 + if self.index_.projection_ == "regular_ll": + lon = ( i - self.index_.known_x_ ) * self.index_.dx_ + self.index_.known_lon_ + lat = ( j - self.index_.known_y_ ) * self.index_.dy_ + self.index_.known_lat_ + + if lon > 180.0: + lon -= 360.0 + + return lat, lon + + def get_tile_extent( self, starti, startj ): + start_lat, start_lon = self.ij_to_latlon( starti - self.index_.tile_bdr_, startj - self.index_.tile_bdr_ ) + stop_lat, stop_lon = self.ij_to_latlon( starti + self.index_.tile_x_ - 1 + self.index_.tile_bdr_, startj + self.index_.tile_y_ - 1 + self.index_.tile_bdr_) + print( ( starti - self.index_.tile_bdr_, startj - self.index_.tile_bdr_ ) ) + print( ( starti + self.index_.tile_x_ - 1 + self.index_.tile_bdr_, startj + self.index_.tile_y_ - 1 + self.index_.tile_bdr_ ) ) + return ( start_lon, stop_lon, start_lat, stop_lat ) + + def get_tile_name( self, lat, lon ): + i, j = self.latlon_to_ij( lat, lon ) + + tile_i = self.index_.tile_x_ * int( int( i ) / self.index_.tile_x_ ) + 1 + tile_j = self.index_.tile_y_ * int( int( j ) / self.index_.tile_y_ ) + 1 + + path = f"{self.source_path_}/{{xstart:0{self.index_.filename_digits_}d}}-{{xstop:0{self.index_.filename_digits_}d}}.{{ystart:0{self.index_.filename_digits_}d}}-{{ystop:0{self.index_.filename_digits_}d}}" + return path.format( + xstart=tile_i, + xstop =tile_i+self.index_.tile_x_-1, + ystart=tile_j, + ystop =tile_j+self.index_.tile_y_-1 + ), tile_i, tile_j + + + def get_box( self, lat, lon, size ): + # X span in points + nx = 0 + if ( np.cos( np.deg2rad( lat ) ) > ( 2.0 * self.pts_per_deg_ * size * 180.0 ) / ( self.npts_x_ * np.pi * self.earth_radius_ ) ): + nx = np.ceil( ( 180.0 * size * self.pts_per_deg_ ) / ( np.pi * self.earth_radius_ * np.cos( np.deg2rad( lat ) ) ) ) + else: + nx = int( self.npts_x_ / 2 ) + + # Y span in points + ny = np.ceil( ( 180.0 * size * self.pts_per_deg_ ) / ( np.pi * self.earth_radius_ ) ) + + + #################################################################################################################### + #################################################################################################################### + #################################################################################################################### + ## Just load the tile in for now + self.data_ = self.read_geogrid( self.get_tile_name( lat, lon )[0] ) + truei, truej = self.latlon_to_ij( lat, lon ) + # Relative i, j in data + reli = truei % self.index_.tile_x_ + relj = truej % self.index_.tile_y_ + + print( reli ) + print( relj ) + print( nx ) + print( ny ) + print( (int(relj-ny/2), int(relj+ny/2)) ) + print( (int(reli-nx/2), int(reli+nx/2)) ) + print( self.data_.shape ) + # Assume it contains all that we need for now + box = self.data_[ 0, int(relj-ny/2):int(relj+ny/2), int(reli-nx/2):int(reli+nx/2) ] + + #################################################################################################################### + #################################################################################################################### + #################################################################################################################### + + return box From 1c79c685fb4ae4ca7b00cd53b0cb0028a5d04943 Mon Sep 17 00:00:00 2001 From: Anthony Islas Date: Tue, 15 Apr 2025 12:56:25 -0700 Subject: [PATCH 04/20] Orographic statistics for GWDO --- util/OrographyStats.py | 86 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 86 insertions(+) create mode 100644 util/OrographyStats.py diff --git a/util/OrographyStats.py b/util/OrographyStats.py new file mode 100644 index 00000000..b39db5ef --- /dev/null +++ b/util/OrographyStats.py @@ -0,0 +1,86 @@ +import numpy as np + + +class OrographyStats: + + def __init__( self, box ): + self.box_ = box + + self.mean_ = np.mean( self.box_ ) #< mean + self.std_ = np.std ( self.box_, mean=self.mean_ ) #< var (actulally stddev) + + #< Critical height used in calculation of orographic effective length + self.hc_ = 1116.2 - 0.878 * self.std_ + + # This currently does not use the landuse to average over only land and zero + # out on mostly water + if self.std_ < 1.0: #< con (convexity) + self.con_ = 0.0 + else: + var4 = np.sum( [ ( h - self.mean_ )**4 for h in np.nditer( box ) ] ) + self.con_ = var4 / ( self.box_.size * self.std_**4 ) + + self.oa_ = np.zeros( (4) ) #< oa (orographic asymmetry) + self.calc_oa() + + self.ol_ = np.zeros( (4) ) # ol (orographic effective length) + self.calc_ol() + + def __repr__( self ): + return f"Mean : {self.mean_} Std : {self.std_} CON : {self.con_} OA : {self.oa_} OL : {self.ol_}" + + def calc_oa( self ): + # Note that right now the assumption is that the box is laid out in col major order + # with [y,x] indices + + # oa1 is the orographic asymmetry in the West direction + nu = np.sum( self.box_[:,:int(self.box_.shape[1]/2)] > self.mean_ ) + nd = np.sum( self.box_[:,int(self.box_.shape[1]/2):] > self.mean_ ) + self.oa_[0] = ( nu -nd ) / ( nu + nd ) if ( ( nu + nd ) > 0 ) else 0.0 + + # oa2 is the orographic asymmetry in the South direction + nu = np.sum( self.box_[:int(self.box_.shape[0]/2),:] > self.mean_ ) + nd = np.sum( self.box_[int(self.box_.shape[0]/2):,:] > self.mean_ ) + self.oa_[1] = ( nu -nd ) / ( nu + nd ) if ( ( nu + nd ) > 0 ) else 0.0 + + # Pre-compute the geometric diagonal of the box + slope = self.box_.shape[1] / self.box_.shape[0] + j, i = np.indices( self.box_.shape ) + upstream = i <= ( j * slope ) + + # oa3 is the orographic asymmetry in the South-West direction + nu = np.sum( self.box_[upstream] > self.mean_ ) + nd = np.sum( self.box_[~upstream] > self.mean_ ) + self.oa_[2] = ( nu -nd ) / ( nu + nd ) if ( ( nu + nd ) > 0 ) else 0.0 + + # oa4 is the orographic asymmetry in the North-West direction + upstream = np.flip( upstream, axis=0 ) + nu = np.sum( self.box_[upstream] > self.mean_ ) + nd = np.sum( self.box_[~upstream] > self.mean_ ) + self.oa_[3] = ( nu -nd ) / ( nu + nd ) if ( ( nu + nd ) > 0 ) else 0.0 + + def calc_ol( self ): + + # ol1 is the effective orographic length in the West direction + interior = self.box_[int(np.floor(self.box_.shape[0]*.25)):int(np.ceil(self.box_.shape[0]*.75)),:] + self.ol_[0] = np.sum( interior > self.hc_ ) / interior.size + + # ol2 is the effective orographic length in the South direction + interior = self.box_[:,int(np.floor(self.box_.shape[1]*.25)):int(np.ceil(self.box_.shape[1]*.75))] + self.ol_[1] = np.sum( interior > self.hc_ ) / interior.size + + # The prescribed methodology uses 4 quadrants to get the diagonals and effectively + # test half of the box, however this does not actually test the interior half + # of the area of the box in the wind direction... + + # ol3 is the effective orographic length in the South-West direction + interiorA = self.box_[int(self.box_.shape[0]/2):,int(self.box_.shape[1]/2):] # first half of x second half of y + interiorB = self.box_[:int(self.box_.shape[0]/2),:int(self.box_.shape[1]/2)] # second half of x first half of y + + self.ol_[2] = ( np.sum( interiorA > self.hc_ ) + np.sum( interiorB > self.hc_ ) ) / ( interiorA.size + interiorB.size ) + + # ol4 is the effective orographic length in the North-West direction + interiorA = self.box_[:int(self.box_.shape[0]/2),:int(self.box_.shape[1]/2)] # second half of x first half of y + interiorB = self.box_[int(self.box_.shape[0]/2):,int(self.box_.shape[1]/2):] # first half of x second half of y + + self.ol_[3] = ( np.sum( interiorA > self.hc_ ) + np.sum( interiorB > self.hc_ ) ) / ( interiorA.size + interiorB.size ) From 9b38f823d79ee9813303f0008292e5164b0ab040 Mon Sep 17 00:00:00 2001 From: Anthony Islas Date: Wed, 16 Apr 2025 18:03:46 -0700 Subject: [PATCH 05/20] Tile loading helper --- util/TileData.py | 105 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 105 insertions(+) create mode 100644 util/TileData.py diff --git a/util/TileData.py b/util/TileData.py new file mode 100644 index 00000000..44aef49f --- /dev/null +++ b/util/TileData.py @@ -0,0 +1,105 @@ +import numpy as np + +class TileData: + + def __init__( self, nx, ny, tile_x, tile_y, load_func ): + # https://docs.python.org/3/faq/programming.html#how-do-i-create-a-multidimensional-list + # Beware multiplying lists + self.tiles_ = [ [None] * nx for y in range( ny ) ] + self.nx_ = nx + self.ny_ = ny + self.tile_x_ = tile_x + self.tile_y_ = tile_y + self.load_func_ = load_func + + def print_present_tile_grid( self ): + for l in self.tiles_: + r = [ 1 if c is not None else 0 for c in l ] + print( r ) + + def load_tile( self, tile_i, tile_j ): + # Convert tile_i and tile_j to full ij index for the start of the tile + # This assumes geogrid start of 1,1 + self.add_tile_by_index( tile_i, tile_j, self.load_func_( tile_i * self.tile_x_, tile_j * self.tile_y_ ) ) + + + def get_tile_index( self, i, j ): + tile_i = int( i / self.tile_x_ ) + tile_j = int( j / self.tile_y_ ) + return tile_i, tile_j + + def add_tile_by_ij( self, i, j, data ): + # i and j are raw indices + tile_i, tile_j = self.get_tile_index( i, j ) + self.add_tile_by_tile_index( tile_i,tile_j, data ) + + def add_tile_by_index( self, tile_i, tile_j, data ): + self.tiles_[tile_j][tile_i] = data + + def get_tiles( self, i, j, tile_indices=np.array([[0,0]]) ): + tiles = [] + tile_i, tile_j = self.get_tile_index( i, j ) + + for ij in tiles: + tiles.append( self.get_tile( tile_j + ij[1], tile_i + ij[0] ) ) + + return tiles + + def get_tile_by_index( self, tile_i, tile_j ): + """ + Go into the tiles and grab the corresponding tile, periodicity is resolved here + """ + # Lat-lon periodicity + true_tile_i = tile_i % self.nx_ + true_tile_j = tile_j + rotate = False + + if tile_j < 0 or tile_j > self.ny_: + # Flip at the poles so return a view of the data rotated 180 + true_tile_i = ( tile_i + self.nx_ / 2 ) % self.nx_ + true_tile_j = self.ny_ - ( tile_j % ( self.ny_ + 1 ) ) + rotate = True + + data = self.tiles_[ true_tile_j ][ true_tile_i ] + if data is None: + print( f"Tile at {true_tile_i} {true_tile_j} is not loaded yet. Loading..." ) + self.load_tile( true_tile_i, true_tile_j ) + # Grab the data again + data = self.tiles_[ true_tile_j ][ true_tile_i ] + else: + print( f"Tile at {true_tile_i} {true_tile_j} is already loaded" ) + + + if rotate: + data = np.rot90( data, 2 ) + + return data + + def get_box( self, indices ): + print( np.einsum( "ijk, i -> ijk", indices, [ 1/self.tile_y_, 1/self.tile_x_ ] ) ) + # First transform indices to regions of distinct tiles by dividing the x and y + # indices by the respective tile sizes + as_tile_idx = np.floor( np.einsum( "ijk, i -> ijk", indices, [ 1/self.tile_y_, 1/self.tile_x_ ] ) ).astype( int ) + # Create our unique set of tiles needed and load them into an analogous structure + uniq_tiles = np.unique( np.unique( as_tile_idx, axis=1 ), axis=2 ).astype( int ) + uniq_tiles_list = list( map( tuple, uniq_tiles.transpose( 1, 2, 0 ).reshape( int( uniq_tiles.size / 2 ), 2 ) ) ) + print( uniq_tiles_list ) + tiles = {} + for tile in uniq_tiles_list: + print( f"Getting data for tile {tile}" ) + tiles[tile] = self.get_tile_by_index( tile[1], tile[0] ) + + box = np.zeros( indices[0].shape, dtype=next(iter(tiles.values())).dtype ) + # Now we have our tile set easily accessible based on the 1-to-1 mapping of indices to tile idx + + for j in range( indices.shape[1] ): + for i in range( indices.shape[2] ): + reli = indices[1,j,i] % self.tile_x_ + relj = indices[0,j,i] % self.tile_y_ + # print( f"i : {i} j : {j}") + # print( f"reli : {reli} relj : {relj}" ) + # print( tiles[tuple(as_tile_idx[:,j,i])] ) + box[j,i] = tiles[tuple(as_tile_idx[:,j,i])][relj,reli] + + self.print_present_tile_grid() + return box From dad08519bc087b74cd239cbcef99a92f22004a82 Mon Sep 17 00:00:00 2001 From: Anthony Islas Date: Wed, 16 Apr 2025 18:04:08 -0700 Subject: [PATCH 06/20] Switch to using tile loader --- util/SourceData.py | 86 ++++++++++++++++++++++++++++++++++++---------- 1 file changed, 67 insertions(+), 19 deletions(-) diff --git a/util/SourceData.py b/util/SourceData.py index 90d803ce..9b149939 100644 --- a/util/SourceData.py +++ b/util/SourceData.py @@ -1,6 +1,8 @@ import numpy as np import cartopy +from TileData import TileData + def get_kv_pair( line, delimiter="=", comment="#" ): non_comment = line.split( comment, maxsplit=1 )[0] @@ -108,6 +110,7 @@ def __init__( self, name, path ): self.index_ = IndexData( self.source_path_ ) self.projection_ = None + self.tile_data_ = None self.npts_x_ = int( 360.0 / self.index_.dx_ ) self.pts_per_deg_ = int( 1.0 / self.index_.dx_ ) self.subgrid_m_dx_ = 2.0 * np.pi * self.earth_radius_ / self.npts_x_ @@ -142,9 +145,44 @@ def __init__( self, name, path ): ) elif self.index_.projection_ == "regular_ll": self.projection_ = cartopy.crs.LambertCylindrical( central_longitude =self.index_.known_lon_ ) + self.tile_data_ = TileData( + int( int( 360.0 / self.index_.dx_ ) / self.index_.tile_x_ ), + int( int( 180.0 / self.index_.dy_ ) / self.index_.tile_y_ ), + self.index_.tile_x_, + self.index_.tile_y_, + load_func=lambda i, j: + self.read_geogrid( + self.get_tile_name_ij( i, j )[0] )[ # Limit the view to just the tile data for now, remove border + 0, + self.index_.tile_bdr_:self.index_.tile_y_+self.index_.tile_bdr_, + self.index_.tile_bdr_:self.index_.tile_x_+self.index_.tile_bdr_ + ] + ) def read_geogrid( self, file ): + """ + Read in the geogrid raw data using the format provided here: + https://www2.mmm.ucar.edu/wrf/users/wrf_users_guide/build/html/wps.html#writing-static-data-to-the-geogrid-binary-format + + Note that data passed out is assumed to be oriented as follows: + SW + <- -x +x -> (increasing longitude) + ^ +---+---+---+---+---+---+ + | | | | | | | | + -y +---+---+---+---+---+---+ + | | | | | | | + +---+---+---+---+---+---+ + | | | | | | | + +---+---+---+---+---+---+ + | | | | | | | + +y +---+---+---+---+---+---+ + | | | | | | | | + v +---+---+---+---+---+---+ + (increasing latitude) + NE + """ + print( "Reading " + file ) rawdata = np.fromfile( file, dtype=self.index_.get_dtype() @@ -163,7 +201,7 @@ def read_geogrid( self, file ): self.index_.tile_x_ + 2 * self.index_.tile_bdr_ ) ) - if self.index_.row_order_ != "top_bottom": + if self.index_.row_order_ == "top_bottom": data = np.flip( data, axis=1 ) return data @@ -200,9 +238,11 @@ def get_tile_extent( self, starti, startj ): print( ( starti + self.index_.tile_x_ - 1 + self.index_.tile_bdr_, startj + self.index_.tile_y_ - 1 + self.index_.tile_bdr_ ) ) return ( start_lon, stop_lon, start_lat, stop_lat ) - def get_tile_name( self, lat, lon ): + def get_tile_name_ll( self, lat, lon ): i, j = self.latlon_to_ij( lat, lon ) + return self.get_tile_name_ij( i, j ) + def get_tile_name_ij( self, i, j ): tile_i = self.index_.tile_x_ * int( int( i ) / self.index_.tile_x_ ) + 1 tile_j = self.index_.tile_y_ * int( int( j ) / self.index_.tile_y_ ) + 1 @@ -219,34 +259,42 @@ def get_box( self, lat, lon, size ): # X span in points nx = 0 if ( np.cos( np.deg2rad( lat ) ) > ( 2.0 * self.pts_per_deg_ * size * 180.0 ) / ( self.npts_x_ * np.pi * self.earth_radius_ ) ): - nx = np.ceil( ( 180.0 * size * self.pts_per_deg_ ) / ( np.pi * self.earth_radius_ * np.cos( np.deg2rad( lat ) ) ) ) + nx = int( np.ceil( ( 180.0 * size * self.pts_per_deg_ ) / ( np.pi * self.earth_radius_ * np.cos( np.deg2rad( lat ) ) ) ) ) else: nx = int( self.npts_x_ / 2 ) # Y span in points - ny = np.ceil( ( 180.0 * size * self.pts_per_deg_ ) / ( np.pi * self.earth_radius_ ) ) + ny = int( np.ceil( ( 180.0 * size * self.pts_per_deg_ ) / ( np.pi * self.earth_radius_ ) ) ) #################################################################################################################### #################################################################################################################### #################################################################################################################### - ## Just load the tile in for now - self.data_ = self.read_geogrid( self.get_tile_name( lat, lon )[0] ) - truei, truej = self.latlon_to_ij( lat, lon ) - # Relative i, j in data - reli = truei % self.index_.tile_x_ - relj = truej % self.index_.tile_y_ + # ## Just load the tile in for now + # self.data_ = self.read_geogrid( self.get_tile_name( lat, lon )[0] ) + true_i, true_j = self.latlon_to_ij( lat, lon ) + # # Relative i, j in data + # reli = truei % self.index_.tile_x_ + # relj = truej % self.index_.tile_y_ + + # print( reli ) + # print( relj ) + # print( nx ) + # print( ny ) + + # print( (int(relj-ny/2), int(relj+ny/2)) ) + # print( (int(reli-nx/2), int(reli+nx/2)) ) + # print( self.data_.shape ) + # # Assume it contains all that we need for now + # box = self.data_[ 0, int(relj-ny/2):int(relj+ny/2), int(reli-nx/2):int(reli+nx/2) ] - print( reli ) - print( relj ) - print( nx ) - print( ny ) - print( (int(relj-ny/2), int(relj+ny/2)) ) - print( (int(reli-nx/2), int(reli+nx/2)) ) - print( self.data_.shape ) - # Assume it contains all that we need for now - box = self.data_[ 0, int(relj-ny/2):int(relj+ny/2), int(reli-nx/2):int(reli+nx/2) ] + # Generate the indices for this box regardless of tile periodicity, let the tile data handle that + indices = np.indices( ( ny, nx ) ) + indices[0] += int( true_j - ny / 2 ) + indices[1] += int( true_i - nx / 2 ) + # print( indices ) + box = self.tile_data_.get_box( indices ) #################################################################################################################### #################################################################################################################### #################################################################################################################### From c4cb4c16d977bb03bc570487bbbdd65c605056ed Mon Sep 17 00:00:00 2001 From: Anthony Islas Date: Thu, 17 Apr 2025 17:24:40 -0700 Subject: [PATCH 07/20] Account for SW start of data (0,0) --- util/OrographyStats.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/util/OrographyStats.py b/util/OrographyStats.py index b39db5ef..2c1b1ae9 100644 --- a/util/OrographyStats.py +++ b/util/OrographyStats.py @@ -39,14 +39,14 @@ def calc_oa( self ): self.oa_[0] = ( nu -nd ) / ( nu + nd ) if ( ( nu + nd ) > 0 ) else 0.0 # oa2 is the orographic asymmetry in the South direction - nu = np.sum( self.box_[:int(self.box_.shape[0]/2),:] > self.mean_ ) - nd = np.sum( self.box_[int(self.box_.shape[0]/2):,:] > self.mean_ ) + nu = np.sum( self.box_[int(self.box_.shape[0]/2):,:] > self.mean_ ) + nd = np.sum( self.box_[:int(self.box_.shape[0]/2),:] > self.mean_ ) self.oa_[1] = ( nu -nd ) / ( nu + nd ) if ( ( nu + nd ) > 0 ) else 0.0 # Pre-compute the geometric diagonal of the box slope = self.box_.shape[1] / self.box_.shape[0] j, i = np.indices( self.box_.shape ) - upstream = i <= ( j * slope ) + upstream = np.flip( i <= ( j * slope ), axis=0 ) # oa3 is the orographic asymmetry in the South-West direction nu = np.sum( self.box_[upstream] > self.mean_ ) @@ -74,13 +74,13 @@ def calc_ol( self ): # of the area of the box in the wind direction... # ol3 is the effective orographic length in the South-West direction - interiorA = self.box_[int(self.box_.shape[0]/2):,int(self.box_.shape[1]/2):] # first half of x second half of y - interiorB = self.box_[:int(self.box_.shape[0]/2),:int(self.box_.shape[1]/2)] # second half of x first half of y + interiorA = self.box_[int(self.box_.shape[0]/2):,int(self.box_.shape[1]/2):] # first half of x first half of y + interiorB = self.box_[:int(self.box_.shape[0]/2),:int(self.box_.shape[1]/2)] # second half of x second half of y self.ol_[2] = ( np.sum( interiorA > self.hc_ ) + np.sum( interiorB > self.hc_ ) ) / ( interiorA.size + interiorB.size ) # ol4 is the effective orographic length in the North-West direction - interiorA = self.box_[:int(self.box_.shape[0]/2),:int(self.box_.shape[1]/2)] # second half of x first half of y - interiorB = self.box_[int(self.box_.shape[0]/2):,int(self.box_.shape[1]/2):] # first half of x second half of y + interiorA = self.box_[int(self.box_.shape[0]/2):,:int(self.box_.shape[1]/2)] # first half of x second half of y + interiorB = self.box_[:int(self.box_.shape[0]/2),int(self.box_.shape[1]/2):] # second half of x first half of y self.ol_[3] = ( np.sum( interiorA > self.hc_ ) + np.sum( interiorB > self.hc_ ) ) / ( interiorA.size + interiorB.size ) From 3e9f90bd6f3f3b35bd0270fddb9c0714788cd232 Mon Sep 17 00:00:00 2001 From: Anthony Islas Date: Thu, 17 Apr 2025 17:25:58 -0700 Subject: [PATCH 08/20] Use faster evaluation of tile points --- util/TileData.py | 41 ++++++++++++++++++++++++++++------------- 1 file changed, 28 insertions(+), 13 deletions(-) diff --git a/util/TileData.py b/util/TileData.py index 44aef49f..1af9b98f 100644 --- a/util/TileData.py +++ b/util/TileData.py @@ -56,7 +56,7 @@ def get_tile_by_index( self, tile_i, tile_j ): if tile_j < 0 or tile_j > self.ny_: # Flip at the poles so return a view of the data rotated 180 - true_tile_i = ( tile_i + self.nx_ / 2 ) % self.nx_ + true_tile_i = int( tile_i + self.nx_ / 2 ) % self.nx_ true_tile_j = self.ny_ - ( tile_j % ( self.ny_ + 1 ) ) rotate = True @@ -71,35 +71,50 @@ def get_tile_by_index( self, tile_i, tile_j ): if rotate: + print( f"Tile at {tile_i} {tile_j} must be rotated" ) data = np.rot90( data, 2 ) return data def get_box( self, indices ): - print( np.einsum( "ijk, i -> ijk", indices, [ 1/self.tile_y_, 1/self.tile_x_ ] ) ) + # First transform indices to regions of distinct tiles by dividing the x and y # indices by the respective tile sizes as_tile_idx = np.floor( np.einsum( "ijk, i -> ijk", indices, [ 1/self.tile_y_, 1/self.tile_x_ ] ) ).astype( int ) + # Create our unique set of tiles needed and load them into an analogous structure - uniq_tiles = np.unique( np.unique( as_tile_idx, axis=1 ), axis=2 ).astype( int ) + # Naive creation by checking for uniqueness + # uniq_tiles = np.unique( np.unique( as_tile_idx, axis=1 ), axis=2 ).astype( int ) + # We know they will increase monotonically so just generate the range + uniq_tiles = np.indices( ( as_tile_idx[0,-1,0] - as_tile_idx[0,0,0] + 1, as_tile_idx[1,0,-1] - as_tile_idx[1,0,0] + 1 ) ) + uniq_tiles[0] += as_tile_idx[0,0,0] + uniq_tiles[1] += as_tile_idx[1,0,0] + uniq_tiles_list = list( map( tuple, uniq_tiles.transpose( 1, 2, 0 ).reshape( int( uniq_tiles.size / 2 ), 2 ) ) ) - print( uniq_tiles_list ) + tiles = {} for tile in uniq_tiles_list: - print( f"Getting data for tile {tile}" ) tiles[tile] = self.get_tile_by_index( tile[1], tile[0] ) box = np.zeros( indices[0].shape, dtype=next(iter(tiles.values())).dtype ) # Now we have our tile set easily accessible based on the 1-to-1 mapping of indices to tile idx - for j in range( indices.shape[1] ): - for i in range( indices.shape[2] ): - reli = indices[1,j,i] % self.tile_x_ - relj = indices[0,j,i] % self.tile_y_ - # print( f"i : {i} j : {j}") - # print( f"reli : {reli} relj : {relj}" ) - # print( tiles[tuple(as_tile_idx[:,j,i])] ) - box[j,i] = tiles[tuple(as_tile_idx[:,j,i])][relj,reli] + for tile in uniq_tiles_list: + print( f"Processing tile {tile}") + # Get the box indices corresponding to this tile + box_indices = np.logical_and( as_tile_idx[0] == tile[0], as_tile_idx[1] == tile[1] ) + + # Bulk assign to box this region of uniq tile by getting the relative index + # into that tile from original indices + box[box_indices] = tiles[tile][indices[0,box_indices] % self.tile_y_, indices[1,box_indices] % self.tile_x_] + + # Naively assign one-by-one + # for j in range( indices.shape[1] ): + # for i in range( indices.shape[2] ): + # reli = indices[1,j,i] % self.tile_x_ + # relj = indices[0,j,i] % self.tile_y_ + # box[j,i] = tiles[tuple(as_tile_idx[:,j,i])][relj,reli] self.print_present_tile_grid() + return box From abcef4b11342c9b76ee0133d79ade40747562f9c Mon Sep 17 00:00:00 2001 From: Anthony Islas Date: Thu, 17 Apr 2025 17:26:06 -0700 Subject: [PATCH 09/20] Clean up --- util/SourceData.py | 17 ----------------- 1 file changed, 17 deletions(-) diff --git a/util/SourceData.py b/util/SourceData.py index 9b149939..d6a45c41 100644 --- a/util/SourceData.py +++ b/util/SourceData.py @@ -271,28 +271,11 @@ def get_box( self, lat, lon, size ): #################################################################################################################### #################################################################################################################### # ## Just load the tile in for now - # self.data_ = self.read_geogrid( self.get_tile_name( lat, lon )[0] ) true_i, true_j = self.latlon_to_ij( lat, lon ) - # # Relative i, j in data - # reli = truei % self.index_.tile_x_ - # relj = truej % self.index_.tile_y_ - - # print( reli ) - # print( relj ) - # print( nx ) - # print( ny ) - - # print( (int(relj-ny/2), int(relj+ny/2)) ) - # print( (int(reli-nx/2), int(reli+nx/2)) ) - # print( self.data_.shape ) - # # Assume it contains all that we need for now - # box = self.data_[ 0, int(relj-ny/2):int(relj+ny/2), int(reli-nx/2):int(reli+nx/2) ] - # Generate the indices for this box regardless of tile periodicity, let the tile data handle that indices = np.indices( ( ny, nx ) ) indices[0] += int( true_j - ny / 2 ) indices[1] += int( true_i - nx / 2 ) - # print( indices ) box = self.tile_data_.get_box( indices ) #################################################################################################################### From 7c0fa7587113bcd52211d9b8e61a2376e0bf9aee Mon Sep 17 00:00:00 2001 From: Anthony Islas Date: Fri, 18 Apr 2025 11:20:41 -0700 Subject: [PATCH 10/20] Add rudimentary cache clearing --- util/TileData.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/util/TileData.py b/util/TileData.py index 1af9b98f..8b444f70 100644 --- a/util/TileData.py +++ b/util/TileData.py @@ -1,3 +1,4 @@ +from collections import deque import numpy as np class TileData: @@ -11,6 +12,8 @@ def __init__( self, nx, ny, tile_x, tile_y, load_func ): self.tile_x_ = tile_x self.tile_y_ = tile_y self.load_func_ = load_func + self.tile_cache_ = deque() + self.max_cache_ = 4 def print_present_tile_grid( self ): for l in self.tiles_: @@ -18,10 +21,22 @@ def print_present_tile_grid( self ): print( r ) def load_tile( self, tile_i, tile_j ): + if len( self.tile_cache_ ) == self.max_cache_: + # make space for incoming tile + tile = self.tile_cache_.popleft() + print( f"Removing tile {tile} from loaded tiles cache" ) + + data = self.tiles_[tile[0]][tile[1]] + self.tiles_[tile[0]][tile[1]] = None + del data + + # Convert tile_i and tile_j to full ij index for the start of the tile # This assumes geogrid start of 1,1 self.add_tile_by_index( tile_i, tile_j, self.load_func_( tile_i * self.tile_x_, tile_j * self.tile_y_ ) ) + self.tile_cache_.append( ( tile_j, tile_i ) ) + def get_tile_index( self, i, j ): tile_i = int( i / self.tile_x_ ) From 8637148dc6ac96b608a7e13dcea8500c0fb5977b Mon Sep 17 00:00:00 2001 From: Anthony Islas Date: Fri, 18 Apr 2025 13:26:17 -0700 Subject: [PATCH 11/20] Suppress debug print statements --- util/SourceData.py | 2 +- util/TileData.py | 13 +++++++------ 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/util/SourceData.py b/util/SourceData.py index d6a45c41..1816495a 100644 --- a/util/SourceData.py +++ b/util/SourceData.py @@ -182,7 +182,7 @@ def read_geogrid( self, file ): (increasing latitude) NE """ - print( "Reading " + file ) + # print( "Reading " + file ) rawdata = np.fromfile( file, dtype=self.index_.get_dtype() diff --git a/util/TileData.py b/util/TileData.py index 8b444f70..06137415 100644 --- a/util/TileData.py +++ b/util/TileData.py @@ -14,6 +14,7 @@ def __init__( self, nx, ny, tile_x, tile_y, load_func ): self.load_func_ = load_func self.tile_cache_ = deque() self.max_cache_ = 4 + self.debug_ = False def print_present_tile_grid( self ): for l in self.tiles_: @@ -24,7 +25,7 @@ def load_tile( self, tile_i, tile_j ): if len( self.tile_cache_ ) == self.max_cache_: # make space for incoming tile tile = self.tile_cache_.popleft() - print( f"Removing tile {tile} from loaded tiles cache" ) + if self.debug_: print( f"Removing tile {tile} from loaded tiles cache" ) data = self.tiles_[tile[0]][tile[1]] self.tiles_[tile[0]][tile[1]] = None @@ -77,16 +78,16 @@ def get_tile_by_index( self, tile_i, tile_j ): data = self.tiles_[ true_tile_j ][ true_tile_i ] if data is None: - print( f"Tile at {true_tile_i} {true_tile_j} is not loaded yet. Loading..." ) + if self.debug_: print( f"Tile at {true_tile_i} {true_tile_j} is not loaded yet. Loading..." ) self.load_tile( true_tile_i, true_tile_j ) # Grab the data again data = self.tiles_[ true_tile_j ][ true_tile_i ] else: - print( f"Tile at {true_tile_i} {true_tile_j} is already loaded" ) + if self.debug_: print( f"Tile at {true_tile_i} {true_tile_j} is already loaded" ) if rotate: - print( f"Tile at {tile_i} {tile_j} must be rotated" ) + if self.debug_: print( f"Tile at {tile_i} {tile_j} must be rotated" ) data = np.rot90( data, 2 ) return data @@ -115,7 +116,7 @@ def get_box( self, indices ): # Now we have our tile set easily accessible based on the 1-to-1 mapping of indices to tile idx for tile in uniq_tiles_list: - print( f"Processing tile {tile}") + if self.debug_: print( f"Processing tile {tile}") # Get the box indices corresponding to this tile box_indices = np.logical_and( as_tile_idx[0] == tile[0], as_tile_idx[1] == tile[1] ) @@ -130,6 +131,6 @@ def get_box( self, indices ): # relj = indices[0,j,i] % self.tile_y_ # box[j,i] = tiles[tuple(as_tile_idx[:,j,i])][relj,reli] - self.print_present_tile_grid() + if self.debug_: self.print_present_tile_grid() return box From 0853cd39185365f0a33905beed1d62be982f76ae Mon Sep 17 00:00:00 2001 From: Anthony Islas Date: Fri, 18 Apr 2025 14:08:45 -0700 Subject: [PATCH 12/20] Simple driver to process geogrid netCDF output for GWDO --- util/compute_gwdo.py | 65 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 65 insertions(+) create mode 100755 util/compute_gwdo.py diff --git a/util/compute_gwdo.py b/util/compute_gwdo.py new file mode 100755 index 00000000..6616c117 --- /dev/null +++ b/util/compute_gwdo.py @@ -0,0 +1,65 @@ +#!/usr/bin/env python3 + +import sys + +import numpy as np +from netCDF4 import Dataset + +from SourceData import SourceData +from OrographyStats import OrographyStats + + +topo_data = "/home/aislas/wrf-model/DATA/WPS_GEOG/topo_gmted2010_30s" +topo_source = SourceData( "gmted2010_30s", topo_data ) + +# For now just get cmdline arg +geo_files = sys.argv[1:] +box_size = 30000 + +for geo in geo_files: + print( f"Processing {geo}" ) + geo_data = Dataset( geo, "r+" ) + + xlat = geo_data.variables[ "XLAT_M" ] + xlon = geo_data.variables[ "XLONG_M" ] + + # Fields to overwrite + con = geo_data.variables[ "CON" ] + var = geo_data.variables[ "VAR" ] + + oa1 = geo_data.variables[ "OA1" ] + oa2 = geo_data.variables[ "OA2" ] + oa3 = geo_data.variables[ "OA3" ] + oa4 = geo_data.variables[ "OA4" ] + + ol1 = geo_data.variables[ "OL1" ] + ol2 = geo_data.variables[ "OL2" ] + ol3 = geo_data.variables[ "OL3" ] + ol4 = geo_data.variables[ "OL4" ] + + ns_size = xlat.shape[1] + we_size = xlat.shape[2] + + for j in range( ns_size ): + for i in range( we_size ): + lat = xlat[0,j,i] + lon = xlon[0,j,i] + box = topo_source.get_box( lat, lon, box_size ) + + oro_stats = OrographyStats( box ) + con[0,j,i] = oro_stats.con_ + var[0,j,i] = oro_stats.std_ + + oa1[0,j,i] = oro_stats.oa_[0] + oa2[0,j,i] = oro_stats.oa_[1] + oa3[0,j,i] = oro_stats.oa_[2] + oa4[0,j,i] = oro_stats.oa_[3] + + ol1[0,j,i] = oro_stats.ol_[0] + ol2[0,j,i] = oro_stats.ol_[1] + ol3[0,j,i] = oro_stats.ol_[2] + ol4[0,j,i] = oro_stats.ol_[3] + + geo_data.close() + +print( "Done!" ) From a05691138d4149bc523507dd007247d15a38daf7 Mon Sep 17 00:00:00 2001 From: Anthony Islas Date: Mon, 5 May 2025 16:35:47 -0700 Subject: [PATCH 13/20] Fix oa2 calculation for nd and nu --- util/OrographyStats.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/util/OrographyStats.py b/util/OrographyStats.py index 2c1b1ae9..b51b5d13 100644 --- a/util/OrographyStats.py +++ b/util/OrographyStats.py @@ -39,8 +39,8 @@ def calc_oa( self ): self.oa_[0] = ( nu -nd ) / ( nu + nd ) if ( ( nu + nd ) > 0 ) else 0.0 # oa2 is the orographic asymmetry in the South direction - nu = np.sum( self.box_[int(self.box_.shape[0]/2):,:] > self.mean_ ) - nd = np.sum( self.box_[:int(self.box_.shape[0]/2),:] > self.mean_ ) + nu = np.sum( self.box_[:int(self.box_.shape[0]/2),:] > self.mean_ ) + nd = np.sum( self.box_[int(self.box_.shape[0]/2):,:] > self.mean_ ) self.oa_[1] = ( nu -nd ) / ( nu + nd ) if ( ( nu + nd ) > 0 ) else 0.0 # Pre-compute the geometric diagonal of the box From 1149ce6889753b1f8b2a931eafd78c9f17610597 Mon Sep 17 00:00:00 2001 From: Anthony Islas Date: Mon, 5 May 2025 16:36:18 -0700 Subject: [PATCH 14/20] Flip variable calcs to match comments, output is the same --- util/OrographyStats.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/util/OrographyStats.py b/util/OrographyStats.py index b51b5d13..c22fbb52 100644 --- a/util/OrographyStats.py +++ b/util/OrographyStats.py @@ -74,8 +74,8 @@ def calc_ol( self ): # of the area of the box in the wind direction... # ol3 is the effective orographic length in the South-West direction - interiorA = self.box_[int(self.box_.shape[0]/2):,int(self.box_.shape[1]/2):] # first half of x first half of y - interiorB = self.box_[:int(self.box_.shape[0]/2),:int(self.box_.shape[1]/2)] # second half of x second half of y + interiorA = self.box_[:int(self.box_.shape[0]/2),:int(self.box_.shape[1]/2)] # first half of x first half of y + interiorB = self.box_[int(self.box_.shape[0]/2):,int(self.box_.shape[1]/2):] # second half of x second half of y self.ol_[2] = ( np.sum( interiorA > self.hc_ ) + np.sum( interiorB > self.hc_ ) ) / ( interiorA.size + interiorB.size ) From 402f9bd4692f511d966089dc0a22efaccfd2d1a0 Mon Sep 17 00:00:00 2001 From: Anthony Islas Date: Mon, 13 Oct 2025 15:53:20 -0700 Subject: [PATCH 15/20] Clean up code and add parsing from namelist --- util/OrographyStats.py | 81 +++++----- util/SourceData.py | 349 +++++++++++++++++++++-------------------- util/TileData.py | 125 ++++++++------- util/compute_gwdo.py | 214 ++++++++++++++++++------- 4 files changed, 441 insertions(+), 328 deletions(-) diff --git a/util/OrographyStats.py b/util/OrographyStats.py index c22fbb52..8402c615 100644 --- a/util/OrographyStats.py +++ b/util/OrographyStats.py @@ -2,85 +2,88 @@ class OrographyStats: - def __init__( self, box ): - self.box_ = box + self.box = box - self.mean_ = np.mean( self.box_ ) #< mean - self.std_ = np.std ( self.box_, mean=self.mean_ ) #< var (actulally stddev) + self.mean = np.mean( self.box ) + # var (actulally stddev) + self.std = np.std( self.box, mean=self.mean ) + self.max = np.max( box ) - #< Critical height used in calculation of orographic effective length - self.hc_ = 1116.2 - 0.878 * self.std_ + # Critical height used in calculation of orographic effective length + self.hc = 1116.2 - 0.878 * self.std # This currently does not use the landuse to average over only land and zero # out on mostly water - if self.std_ < 1.0: #< con (convexity) - self.con_ = 0.0 + # con (convexity) + if self.std < 1.0: + self.con = 0.0 else: - var4 = np.sum( [ ( h - self.mean_ )**4 for h in np.nditer( box ) ] ) - self.con_ = var4 / ( self.box_.size * self.std_**4 ) - - self.oa_ = np.zeros( (4) ) #< oa (orographic asymmetry) + var4 = np.sum( [ ( h - self.mean )**4 for h in np.nditer( box ) ] ) + self.con = var4 / ( self.box.size * self.std**4 ) + + # oa (orographic asymmetry) + self.oa = np.zeros( (4) ) self.calc_oa() - self.ol_ = np.zeros( (4) ) # ol (orographic effective length) + # ol (orographic effective length) + self.ol = np.zeros( (4) ) self.calc_ol() def __repr__( self ): - return f"Mean : {self.mean_} Std : {self.std_} CON : {self.con_} OA : {self.oa_} OL : {self.ol_}" + return f"Mean : {self.mean} Std : {self.std} Max: {self.max} CON : {self.con} OA : {self.oa} OL : {self.ol}" def calc_oa( self ): # Note that right now the assumption is that the box is laid out in col major order # with [y,x] indices # oa1 is the orographic asymmetry in the West direction - nu = np.sum( self.box_[:,:int(self.box_.shape[1]/2)] > self.mean_ ) - nd = np.sum( self.box_[:,int(self.box_.shape[1]/2):] > self.mean_ ) - self.oa_[0] = ( nu -nd ) / ( nu + nd ) if ( ( nu + nd ) > 0 ) else 0.0 + nu = np.sum( self.box[:, :int(self.box.shape[1] / 2)] > self.mean ) + nd = np.sum( self.box[:, int(self.box.shape[1] / 2):] > self.mean ) + self.oa[0] = ( nu - nd ) / ( nu + nd ) if ( ( nu + nd ) > 0 ) else 0.0 # oa2 is the orographic asymmetry in the South direction - nu = np.sum( self.box_[:int(self.box_.shape[0]/2),:] > self.mean_ ) - nd = np.sum( self.box_[int(self.box_.shape[0]/2):,:] > self.mean_ ) - self.oa_[1] = ( nu -nd ) / ( nu + nd ) if ( ( nu + nd ) > 0 ) else 0.0 + nu = np.sum( self.box[:int(self.box.shape[0] / 2), :] > self.mean ) + nd = np.sum( self.box[int(self.box.shape[0] / 2):, :] > self.mean ) + self.oa[1] = ( nu - nd ) / ( nu + nd ) if ( ( nu + nd ) > 0 ) else 0.0 # Pre-compute the geometric diagonal of the box - slope = self.box_.shape[1] / self.box_.shape[0] - j, i = np.indices( self.box_.shape ) + slope = self.box.shape[1] / self.box.shape[0] + j, i = np.indices( self.box.shape ) upstream = np.flip( i <= ( j * slope ), axis=0 ) # oa3 is the orographic asymmetry in the South-West direction - nu = np.sum( self.box_[upstream] > self.mean_ ) - nd = np.sum( self.box_[~upstream] > self.mean_ ) - self.oa_[2] = ( nu -nd ) / ( nu + nd ) if ( ( nu + nd ) > 0 ) else 0.0 + nu = np.sum( self.box[upstream] > self.mean ) + nd = np.sum( self.box[~upstream] > self.mean ) + self.oa[2] = ( nu - nd ) / ( nu + nd ) if ( ( nu + nd ) > 0 ) else 0.0 # oa4 is the orographic asymmetry in the North-West direction upstream = np.flip( upstream, axis=0 ) - nu = np.sum( self.box_[upstream] > self.mean_ ) - nd = np.sum( self.box_[~upstream] > self.mean_ ) - self.oa_[3] = ( nu -nd ) / ( nu + nd ) if ( ( nu + nd ) > 0 ) else 0.0 + nu = np.sum( self.box[upstream] > self.mean ) + nd = np.sum( self.box[~upstream] > self.mean ) + self.oa[3] = ( nu - nd ) / ( nu + nd ) if ( ( nu + nd ) > 0 ) else 0.0 def calc_ol( self ): - # ol1 is the effective orographic length in the West direction - interior = self.box_[int(np.floor(self.box_.shape[0]*.25)):int(np.ceil(self.box_.shape[0]*.75)),:] - self.ol_[0] = np.sum( interior > self.hc_ ) / interior.size + interior = self.box[int(np.floor(self.box.shape[0] * .25)):int(np.ceil(self.box.shape[0] * .75)), :] + self.ol[0] = np.sum( interior > self.hc ) / interior.size # ol2 is the effective orographic length in the South direction - interior = self.box_[:,int(np.floor(self.box_.shape[1]*.25)):int(np.ceil(self.box_.shape[1]*.75))] - self.ol_[1] = np.sum( interior > self.hc_ ) / interior.size + interior = self.box[:, int(np.floor(self.box.shape[1] * .25)):int(np.ceil(self.box.shape[1] * .75))] + self.ol[1] = np.sum( interior > self.hc ) / interior.size # The prescribed methodology uses 4 quadrants to get the diagonals and effectively # test half of the box, however this does not actually test the interior half # of the area of the box in the wind direction... # ol3 is the effective orographic length in the South-West direction - interiorA = self.box_[:int(self.box_.shape[0]/2),:int(self.box_.shape[1]/2)] # first half of x first half of y - interiorB = self.box_[int(self.box_.shape[0]/2):,int(self.box_.shape[1]/2):] # second half of x second half of y + interiorA = self.box[:int(self.box.shape[0] / 2), :int(self.box.shape[1] / 2)] # first half of x first half of y + interiorB = self.box[int(self.box.shape[0] / 2):, int(self.box.shape[1] / 2):] # second half of x second half of y - self.ol_[2] = ( np.sum( interiorA > self.hc_ ) + np.sum( interiorB > self.hc_ ) ) / ( interiorA.size + interiorB.size ) + self.ol[2] = ( np.sum( interiorA > self.hc ) + np.sum( interiorB > self.hc ) ) / ( interiorA.size + interiorB.size ) # ol4 is the effective orographic length in the North-West direction - interiorA = self.box_[int(self.box_.shape[0]/2):,:int(self.box_.shape[1]/2)] # first half of x second half of y - interiorB = self.box_[:int(self.box_.shape[0]/2),int(self.box_.shape[1]/2):] # second half of x first half of y + interiorA = self.box[int(self.box.shape[0] / 2):, :int(self.box.shape[1] / 2)] # first half of x second half of y + interiorB = self.box[:int(self.box.shape[0] / 2), int(self.box.shape[1] / 2):] # second half of x first half of y - self.ol_[3] = ( np.sum( interiorA > self.hc_ ) + np.sum( interiorB > self.hc_ ) ) / ( interiorA.size + interiorB.size ) + self.ol[3] = ( np.sum( interiorA > self.hc ) + np.sum( interiorB > self.hc ) ) / ( interiorA.size + interiorB.size ) diff --git a/util/SourceData.py b/util/SourceData.py index 1816495a..8f4a38b2 100644 --- a/util/SourceData.py +++ b/util/SourceData.py @@ -3,6 +3,7 @@ from TileData import TileData + def get_kv_pair( line, delimiter="=", comment="#" ): non_comment = line.split( comment, maxsplit=1 )[0] @@ -12,153 +13,149 @@ def get_kv_pair( line, delimiter="=", comment="#" ): else: return None, None -class IndexData: +class IndexData: def __init__( self, path ): - self.source_path_ = path + self.source_path = path # Default values of index file from : # https://www2.mmm.ucar.edu/wrf/users/wrf_users_guide/build/html/wps.html#index-options - self.projection_ = "required" - self.source_type_ = "required" - self.signed_ = False - self.units_ = "required" - self.description_ = "required" - self.dx_ = -1.0 # required - self.dy_ = -1.0 # required - self.known_x_ = 1.0 - self.known_y_ = 1.0 - self.known_lat_ = -1.0 # required - self.known_lon_ = -1.0 # required - self.stdlon_ = -2.0 # not required - self.truelat1_ = -2.0 # not required - self.truelat2_ = -2.0 # not required - self.wordsize_ = -1 # required - self.tile_x_ = -1 # required - self.tile_y_ = -1 # required - self.tile_z_ = -2 # not required - self.tile_z_start_ = -2 # not required - self.tile_z_end_ = -2 # not required - self.category_min_ = -2 # not required - self.category_max_ = -2 # not required - self.tile_bdr_ = 0 - self.missing_value_ = -2.0 # not required - self.scale_factor_ = 1.0 - self.row_order_ = "bottom_top" - self.endian_ = "big" - self.iswater_ = 16 - self.islake_ = -1 # i.e. no separate inland water category (does not mean required) - self.isice_ = 24 - self.isurban_ = 1 - self.isoilwater_ = 14 - self.mminlu_ = "USGS" - self.filename_digits_ = 5 + self.projection = "required" + self.source_type = "required" + self.signed = False + self.units = "required" + self.description = "required" + self.dx = -1.0 # required + self.dy = -1.0 # required + self.known_x = 1.0 + self.known_y = 1.0 + self.known_lat = -1.0 # required + self.known_lon = -1.0 # required + self.stdlon = -2.0 # not required + self.truelat1 = -2.0 # not required + self.truelat2 = -2.0 # not required + self.wordsize = -1 # required + self.tile_x = -1 # required + self.tile_y = -1 # required + self.tile_z = -2 # not required + self.tile_z_start = -2 # not required + self.tile_z_end = -2 # not required + self.category_min = -2 # not required + self.category_max = -2 # not required + self.tile_bdr = 0 + self.missing_value = -2.0 # not required + self.scale_factor = 1.0 + self.row_order = "bottom_top" + self.endian = "big" + self.iswater = 16 + self.islake = -1 # i.e. no separate inland water category (does not mean required) + self.isice = 24 + self.isurban = 1 + self.isoilwater = 14 + self.mminlu = "USGS" + self.filename_digits = 5 self.read_file() def read_file( self ): - with open( self.source_path_ + "/index" ) as f: + with open( self.source_path + "/index" ) as f: for line in f: key, value = get_kv_pair( line ) if key is not None: - if key == "projection" : self.projection_ = value - elif key == "source_type" : self.source_type_ = value - elif key == "signed" : self.signed_ = ( value == "yes" ) - elif key == "units" : self.units_ = value - elif key == "description" : self.description_ = value - elif key == "dx" : self.dx_ = float( value ) - elif key == "dy" : self.dy_ = float( value ) - elif key == "known_x" : self.known_x_ = float( value ) - elif key == "known_y" : self.known_y_ = float( value ) - elif key == "known_lat" : self.known_lat_ = float( value ) - elif key == "known_lon" : self.known_lon_ = float( value ) - elif key == "stdlon" : self.stdlon_ = float( value ) - elif key == "truelat1" : self.truelat1_ = float( value ) - elif key == "truelat2" : self.truelat2_ = float( value ) - elif key == "wordsize" : self.wordsize_ = int( value ) - elif key == "tile_x" : self.tile_x_ = int( value ) - elif key == "tile_y" : self.tile_y_ = int( value ) - elif key == "tile_z" : self.tile_z_ = int( value ) - elif key == "tile_z_start" : self.tile_z_start_ = int( value ) - elif key == "tile_z_end" : self.tile_z_end_ = int( value ) - elif key == "category_min" : self.category_min_ = int( value ) - elif key == "category_max" : self.category_max_ = int( value ) - elif key == "tile_bdr" : self.tile_bdr_ = int( value ) - elif key == "missing_value" : self.missing_value_ = float( value ) - elif key == "scale_factor" : self.scale_factor_ = float( value ) - elif key == "row_order" : self.row_order_ = value - elif key == "endian" : self.endian_ = value - elif key == "iswater" : self.iswater_ = int( value ) - elif key == "islake" : self.islake_ = int( value ) - elif key == "isice" : self.isice_ = int( value ) - elif key == "isurban" : self.isurban_ = int( value ) - elif key == "isoilwater" : self.isoilwater_ = int( value ) - elif key == "mminlu" : self.mminlu_ = value - elif key == "filename_digits" : self.filename_digits_ = int( value ) - + if key == "projection": self.projection = value + elif key == "source_type": self.source_type = value + elif key == "signed": self.signed = ( value == "yes" ) + elif key == "units": self.units = value + elif key == "description": self.description = value + elif key == "dx": self.dx = float( value ) + elif key == "dy": self.dy = float( value ) + elif key == "known_x": self.known_x = float( value ) + elif key == "known_y": self.known_y = float( value ) + elif key == "known_lat": self.known_lat = float( value ) + elif key == "known_lon": self.known_lon = float( value ) + elif key == "stdlon": self.stdlon = float( value ) + elif key == "truelat1": self.truelat1 = float( value ) + elif key == "truelat2": self.truelat2 = float( value ) + elif key == "wordsize": self.wordsize = int( value ) + elif key == "tile_x": self.tile_x = int( value ) + elif key == "tile_y": self.tile_y = int( value ) + elif key == "tile_z": self.tile_z = int( value ) + elif key == "tile_z_start": self.tile_z_start = int( value ) + elif key == "tile_z_end": self.tile_z_end = int( value ) + elif key == "category_min": self.category_min = int( value ) + elif key == "category_max": self.category_max = int( value ) + elif key == "tile_bdr": self.tile_bdr = int( value ) + elif key == "missing_value": self.missing_value = float( value ) + elif key == "scale_factor": self.scale_factor = float( value ) + elif key == "row_order": self.row_order = value + elif key == "endian": self.endian = value + elif key == "iswater": self.iswater = int( value ) + elif key == "islake": self.islake = int( value ) + elif key == "isice": self.isice = int( value ) + elif key == "isurban": self.isurban = int( value ) + elif key == "isoilwater": self.isoilwater = int( value ) + elif key == "mminlu": self.mminlu = value + elif key == "filename_digits": self.filename_digits = int( value ) + def get_dtype( self ): - endian = ">" if self.endian_ == "big" else "<" - signed = "i" if self.signed_ else "u" - return f"{endian}{signed}{self.wordsize_}" + endian = ">" if self.endian == "big" else "<" + signed = "i" if self.signed else "u" + return f"{endian}{signed}{self.wordsize}" class SourceData: - def __init__( self, name, path ): - self.earth_radius_ = 6371229.0 - self.name_ = name - self.source_path_ = path - self.index_ = IndexData( self.source_path_ ) - - self.projection_ = None - self.tile_data_ = None - self.npts_x_ = int( 360.0 / self.index_.dx_ ) - self.pts_per_deg_ = int( 1.0 / self.index_.dx_ ) - self.subgrid_m_dx_ = 2.0 * np.pi * self.earth_radius_ / self.npts_x_ - - self.data_ = None - - - # These are all probably wrong - if self.index_.projection_ == "lambert": - self.projection_ = cartopy.crs.LambertConformal( - central_longitude =self.index_.known_lon_, - central_latitude =self.index_.known_lat_, - standard_parallels =( self.index_.truelat1_, self.index_.truelat2_ ) - ) - elif self.index_.projection_ == "polar_wgs84" or self.index_.projection_ == "polar": - self.projection_ = cartopy.crs.Stereographic( - central_longitude =self.index_.known_lon_, - central_latitude =self.index_.known_lat_, - true_scale_latitude =self.index_.truelat1_ - ) - elif self.index_.projection_ == "albers_nad83": - self.projection_ = cartopy.crs.AlbersEqualArea( - central_longitude =self.index_.known_lon_, - central_latitude =self.index_.known_lat_, - standard_parallels =( self.index_.truelat1_, self.index_.truelat2_ ) - ) - elif self.index_.projection_ == "mercator": - self.projection_ = cartopy.crs.Mercator( - central_longitude =self.index_.known_lon_, - # central_latitude =self.index_.known_lat_, - latitude_true_scale =self.index_.truelat1_ - ) - elif self.index_.projection_ == "regular_ll": - self.projection_ = cartopy.crs.LambertCylindrical( central_longitude =self.index_.known_lon_ ) - self.tile_data_ = TileData( - int( int( 360.0 / self.index_.dx_ ) / self.index_.tile_x_ ), - int( int( 180.0 / self.index_.dy_ ) / self.index_.tile_y_ ), - self.index_.tile_x_, - self.index_.tile_y_, - load_func=lambda i, j: - self.read_geogrid( - self.get_tile_name_ij( i, j )[0] )[ # Limit the view to just the tile data for now, remove border - 0, - self.index_.tile_bdr_:self.index_.tile_y_+self.index_.tile_bdr_, - self.index_.tile_bdr_:self.index_.tile_x_+self.index_.tile_bdr_ - ] - ) - + self._earth_radius = 6371229.0 + self._name = name + self._source_path = path + self._index = IndexData( self._source_path ) + + self._projection = None + self._npts_x = int( 360.0 / self._index.dx ) + self._pts_per_deg = int( 1.0 / self._index.dx ) + self._subgrid_m_dx = 2.0 * np.pi * self._earth_radius / self._npts_x + + # # These are all probably wrong + # if self._index.projection_ == "lambert": + # self._projection = cartopy.crs.LambertConformal( + # central_longitude =self._index.known_lon_, + # central_latitude =self._index.known_lat_, + # standard_parallels =( self._index.truelat1_, self._index.truelat2_ ) + # ) + # elif self._index.projection_ == "polar_wgs84" or self._index.projection_ == "polar": + # self._projection = cartopy.crs.Stereographic( + # central_longitude =self._index.known_lon_, + # central_latitude =self._index.known_lat_, + # true_scale_latitude =self._index.truelat1_ + # ) + # elif self._index.projection_ == "albers_nad83": + # self._projection = cartopy.crs.AlbersEqualArea( + # central_longitude =self._index.known_lon_, + # central_latitude =self._index.known_lat_, + # standard_parallels =( self._index.truelat1_, self._index.truelat2_ ) + # ) + # elif self._index.projection_ == "mercator": + # self._projection = cartopy.crs.Mercator( + # central_longitude =self._index.known_lon_, + # # central_latitude =self._index.known_lat_, + # latitude_true_scale =self._index.truelat1_ + # ) + # elif self._index.projection_ == "regular_ll": + # self._projection = cartopy.crs.LambertCylindrical( central_longitude =self._index.known_lon_ ) + + self._tile_data = TileData( + int( int( 360.0 / self._index.dx ) / self._index.tile_x ), + int( int( 180.0 / self._index.dy ) / self._index.tile_y ), + self._index.tile_x, + self._index.tile_y, + load_func=lambda i, j: + self.read_geogrid( + # Limit the view to just the tile data for now, remove border + self.get_tile_name_ij( i, j )[0] )[ + 0, + self._index.tile_bdr:self._index.tile_y + self._index.tile_bdr, + self._index.tile_bdr:self._index.tile_x + self._index.tile_bdr + ] + ) def read_geogrid( self, file ): """ @@ -183,48 +180,48 @@ def read_geogrid( self, file ): NE """ # print( "Reading " + file ) - rawdata = np.fromfile( + rawdata = np.fromfile( file, - dtype=self.index_.get_dtype() + dtype=self._index.get_dtype() ) - data = rawdata.astype( np.float32, casting="unsafe" ) * self.index_.scale_factor_ + data = rawdata.astype( np.float32, casting="unsafe" ) * self._index.scale_factor z_dim = 1 - if self.index_.tile_z_start_ > 0 and self.index_.tile_z_end_ > 0: - z_dim = self.index_.tile_z_end_ - self.index_.tile_z_start_ - elif self.index_.tile_z_ > 0: - z_dim = self.index_.tile_z_ + if self._index.tile_z_start > 0 and self._index.tile_z_end > 0: + z_dim = self._index.tile_z_end - self._index.tile_z_start + elif self._index.tile_z > 0: + z_dim = self._index.tile_z data = data.reshape( ( z_dim, - self.index_.tile_y_ + 2 * self.index_.tile_bdr_, - self.index_.tile_x_ + 2 * self.index_.tile_bdr_ + self._index.tile_y + 2 * self._index.tile_bdr, + self._index.tile_x + 2 * self._index.tile_bdr ) ) - if self.index_.row_order_ == "top_bottom": + if self._index.row_order == "top_bottom": data = np.flip( data, axis=1 ) - + return data def latlon_to_ij( self, lat, lon ): # Ignore staggering for now i = 0.0 j = 0.0 - if self.index_.projection_ == "regular_ll": - delta_lat = lat - self.index_.known_lat_ - delta_lon = lon - self.index_.known_lon_ + if self._index.projection == "regular_ll": + delta_lat = lat - self._index.known_lat + delta_lon = lon - self._index.known_lon - i = ( delta_lon / self.index_.dx_ + self.index_.known_x_ ) % int( 360.0 / self.index_.dx_ ) - j = ( delta_lat / self.index_.dy_ + self.index_.known_y_ ) + i = ( delta_lon / self._index.dx + self._index.known_x ) % int( 360.0 / self._index.dx ) + j = ( delta_lat / self._index.dy + self._index.known_y ) return i, j def ij_to_latlon( self, i, j ): lat = 0.0 lon = 0.0 - if self.index_.projection_ == "regular_ll": - lon = ( i - self.index_.known_x_ ) * self.index_.dx_ + self.index_.known_lon_ - lat = ( j - self.index_.known_y_ ) * self.index_.dy_ + self.index_.known_lat_ + if self._index.projection == "regular_ll": + lon = ( i - self._index.known_x ) * self._index.dx + self._index.known_lon + lat = ( j - self._index.known_y ) * self._index.dy + self._index.known_lat if lon > 180.0: lon -= 360.0 @@ -232,10 +229,10 @@ def ij_to_latlon( self, i, j ): return lat, lon def get_tile_extent( self, starti, startj ): - start_lat, start_lon = self.ij_to_latlon( starti - self.index_.tile_bdr_, startj - self.index_.tile_bdr_ ) - stop_lat, stop_lon = self.ij_to_latlon( starti + self.index_.tile_x_ - 1 + self.index_.tile_bdr_, startj + self.index_.tile_y_ - 1 + self.index_.tile_bdr_) - print( ( starti - self.index_.tile_bdr_, startj - self.index_.tile_bdr_ ) ) - print( ( starti + self.index_.tile_x_ - 1 + self.index_.tile_bdr_, startj + self.index_.tile_y_ - 1 + self.index_.tile_bdr_ ) ) + start_lat, start_lon = self.ij_to_latlon( starti - self._index.tile_bdr, startj - self._index.tile_bdr ) + stop_lat, stop_lon = self.ij_to_latlon( starti + self._index.tile_x - 1 + self._index.tile_bdr, startj + self._index.tile_y - 1 + self._index.tile_bdr ) + print( ( starti - self._index.tile_bdr, startj - self._index.tile_bdr ) ) + print( ( starti + self._index.tile_x - 1 + self._index.tile_bdr, startj + self._index.tile_y - 1 + self._index.tile_bdr ) ) return ( start_lon, stop_lon, start_lat, stop_lat ) def get_tile_name_ll( self, lat, lon ): @@ -243,43 +240,47 @@ def get_tile_name_ll( self, lat, lon ): return self.get_tile_name_ij( i, j ) def get_tile_name_ij( self, i, j ): - tile_i = self.index_.tile_x_ * int( int( i ) / self.index_.tile_x_ ) + 1 - tile_j = self.index_.tile_y_ * int( int( j ) / self.index_.tile_y_ ) + 1 + tile_i = self._index.tile_x * int( int( i ) / self._index.tile_x ) + 1 + tile_j = self._index.tile_y * int( int( j ) / self._index.tile_y ) + 1 - path = f"{self.source_path_}/{{xstart:0{self.index_.filename_digits_}d}}-{{xstop:0{self.index_.filename_digits_}d}}.{{ystart:0{self.index_.filename_digits_}d}}-{{ystop:0{self.index_.filename_digits_}d}}" + path = f"{self._source_path}/" + path += f"{{xstart:0{self._index.filename_digits}d}}-{{xstop:0{self._index.filename_digits}d}}." + path += f"{{ystart:0{self._index.filename_digits}d}}-{{ystop:0{self._index.filename_digits}d}}" return path.format( xstart=tile_i, - xstop =tile_i+self.index_.tile_x_-1, + xstop=tile_i + self._index.tile_x - 1, ystart=tile_j, - ystop =tile_j+self.index_.tile_y_-1 + ystop=tile_j + self._index.tile_y - 1 ), tile_i, tile_j + def get_box( self, lat, lon, size_x, size_y=None ): + if size_y is None: + size_y = size_x - def get_box( self, lat, lon, size ): # X span in points nx = 0 - if ( np.cos( np.deg2rad( lat ) ) > ( 2.0 * self.pts_per_deg_ * size * 180.0 ) / ( self.npts_x_ * np.pi * self.earth_radius_ ) ): - nx = int( np.ceil( ( 180.0 * size * self.pts_per_deg_ ) / ( np.pi * self.earth_radius_ * np.cos( np.deg2rad( lat ) ) ) ) ) + if ( + np.cos( np.deg2rad( lat ) ) + > ( 2.0 * self._pts_per_deg * size_x * 180.0 ) / ( self._npts_x * np.pi * self._earth_radius ) + ): + nx = int( + np.ceil( + ( 180.0 * size_x * self._pts_per_deg ) + / ( np.pi * self._earth_radius * np.cos( np.deg2rad( lat ) ) ) + ) + ) else: - nx = int( self.npts_x_ / 2 ) + nx = int( self._npts_x / 2 ) # Y span in points - ny = int( np.ceil( ( 180.0 * size * self.pts_per_deg_ ) / ( np.pi * self.earth_radius_ ) ) ) + ny = int( np.ceil( ( 180.0 * size_y * self._pts_per_deg ) / ( np.pi * self._earth_radius ) ) ) - - #################################################################################################################### - #################################################################################################################### - #################################################################################################################### - # ## Just load the tile in for now true_i, true_j = self.latlon_to_ij( lat, lon ) # Generate the indices for this box regardless of tile periodicity, let the tile data handle that indices = np.indices( ( ny, nx ) ) indices[0] += int( true_j - ny / 2 ) indices[1] += int( true_i - nx / 2 ) - box = self.tile_data_.get_box( indices ) - #################################################################################################################### - #################################################################################################################### - #################################################################################################################### - + box = self._tile_data.get_box( indices ) + # box is now a 2D box of data spanning size in meters centered on lat/lon return box diff --git a/util/TileData.py b/util/TileData.py index 06137415..098a573a 100644 --- a/util/TileData.py +++ b/util/TileData.py @@ -1,136 +1,143 @@ from collections import deque import numpy as np -class TileData: +class TileData: def __init__( self, nx, ny, tile_x, tile_y, load_func ): # https://docs.python.org/3/faq/programming.html#how-do-i-create-a-multidimensional-list # Beware multiplying lists - self.tiles_ = [ [None] * nx for y in range( ny ) ] - self.nx_ = nx - self.ny_ = ny - self.tile_x_ = tile_x - self.tile_y_ = tile_y - self.load_func_ = load_func - self.tile_cache_ = deque() - self.max_cache_ = 4 - self.debug_ = False - + self._tiles = [ [None] * nx for y in range( ny ) ] + self._nx = nx + self._ny = ny + self._tile_x = tile_x + self._tile_y = tile_y + self._load_func = load_func + self._tile_cache = deque() + self._max_cache = 4 + self._debug = False + def print_present_tile_grid( self ): - for l in self.tiles_: - r = [ 1 if c is not None else 0 for c in l ] + for tile_list in self._tiles: + r = [ 1 if c is not None else 0 for c in tile_list ] print( r ) def load_tile( self, tile_i, tile_j ): - if len( self.tile_cache_ ) == self.max_cache_: + if len( self._tile_cache ) == self._max_cache: # make space for incoming tile - tile = self.tile_cache_.popleft() - if self.debug_: print( f"Removing tile {tile} from loaded tiles cache" ) + tile = self._tile_cache.popleft() + if self._debug: + print( f"Removing tile {tile} from loaded tiles cache" ) - data = self.tiles_[tile[0]][tile[1]] - self.tiles_[tile[0]][tile[1]] = None + data = self._tiles[tile[0]][tile[1]] + self._tiles[tile[0]][tile[1]] = None del data - # Convert tile_i and tile_j to full ij index for the start of the tile # This assumes geogrid start of 1,1 - self.add_tile_by_index( tile_i, tile_j, self.load_func_( tile_i * self.tile_x_, tile_j * self.tile_y_ ) ) - - self.tile_cache_.append( ( tile_j, tile_i ) ) + self.add_tile_by_index( tile_i, tile_j, self._load_func( tile_i * self._tile_x, tile_j * self._tile_y ) ) + self._tile_cache.append( ( tile_j, tile_i ) ) def get_tile_index( self, i, j ): - tile_i = int( i / self.tile_x_ ) - tile_j = int( j / self.tile_y_ ) + tile_i = int( i / self._tile_x ) + tile_j = int( j / self._tile_y ) return tile_i, tile_j def add_tile_by_ij( self, i, j, data ): # i and j are raw indices tile_i, tile_j = self.get_tile_index( i, j ) - self.add_tile_by_tile_index( tile_i,tile_j, data ) - + self.add_tile_by_tile_index( tile_i, tile_j, data ) + def add_tile_by_index( self, tile_i, tile_j, data ): - self.tiles_[tile_j][tile_i] = data + self._tiles[tile_j][tile_i] = data - def get_tiles( self, i, j, tile_indices=np.array([[0,0]]) ): + def get_tiles( self, i, j, tile_indices=np.array([[0, 0]]) ): tiles = [] tile_i, tile_j = self.get_tile_index( i, j ) for ij in tiles: tiles.append( self.get_tile( tile_j + ij[1], tile_i + ij[0] ) ) - + return tiles - + def get_tile_by_index( self, tile_i, tile_j ): """ Go into the tiles and grab the corresponding tile, periodicity is resolved here """ # Lat-lon periodicity - true_tile_i = tile_i % self.nx_ + true_tile_i = tile_i % self._nx true_tile_j = tile_j rotate = False - if tile_j < 0 or tile_j > self.ny_: + if tile_j < 0 or tile_j > self._ny: # Flip at the poles so return a view of the data rotated 180 - true_tile_i = int( tile_i + self.nx_ / 2 ) % self.nx_ - true_tile_j = self.ny_ - ( tile_j % ( self.ny_ + 1 ) ) + true_tile_i = int( tile_i + self._nx / 2 ) % self._nx + true_tile_j = self._ny - ( tile_j % ( self._ny + 1 ) ) rotate = True - data = self.tiles_[ true_tile_j ][ true_tile_i ] + data = self._tiles[ true_tile_j ][ true_tile_i ] if data is None: - if self.debug_: print( f"Tile at {true_tile_i} {true_tile_j} is not loaded yet. Loading..." ) + if self._debug: + print( f"Tile at {true_tile_i} {true_tile_j} is not loaded yet. Loading..." ) self.load_tile( true_tile_i, true_tile_j ) # Grab the data again - data = self.tiles_[ true_tile_j ][ true_tile_i ] + data = self._tiles[ true_tile_j ][ true_tile_i ] else: - if self.debug_: print( f"Tile at {true_tile_i} {true_tile_j} is already loaded" ) + if self._debug: + print( f"Tile at {true_tile_i} {true_tile_j} is already loaded" ) - if rotate: - if self.debug_: print( f"Tile at {tile_i} {tile_j} must be rotated" ) + if self._debug: + print( f"Tile at {tile_i} {tile_j} must be rotated" ) data = np.rot90( data, 2 ) return data def get_box( self, indices ): - - # First transform indices to regions of distinct tiles by dividing the x and y + # First transform indices to regions of distinct tiles by dividing the x and y # indices by the respective tile sizes - as_tile_idx = np.floor( np.einsum( "ijk, i -> ijk", indices, [ 1/self.tile_y_, 1/self.tile_x_ ] ) ).astype( int ) + as_tile_idx = np.floor( + np.einsum( + "ijk, i -> ijk", + indices, + [ 1 / self._tile_y, 1 / self._tile_x ] + ) + ).astype( int ) # Create our unique set of tiles needed and load them into an analogous structure # Naive creation by checking for uniqueness # uniq_tiles = np.unique( np.unique( as_tile_idx, axis=1 ), axis=2 ).astype( int ) + # We know they will increase monotonically so just generate the range - uniq_tiles = np.indices( ( as_tile_idx[0,-1,0] - as_tile_idx[0,0,0] + 1, as_tile_idx[1,0,-1] - as_tile_idx[1,0,0] + 1 ) ) - uniq_tiles[0] += as_tile_idx[0,0,0] - uniq_tiles[1] += as_tile_idx[1,0,0] + uniq_tiles = np.indices( + ( + as_tile_idx[0, -1, 0] - as_tile_idx[0, 0, 0] + 1, + as_tile_idx[1, 0, -1] - as_tile_idx[1, 0, 0] + 1 + ) + ) + uniq_tiles[0] += as_tile_idx[0, 0, 0] + uniq_tiles[1] += as_tile_idx[1, 0, 0] uniq_tiles_list = list( map( tuple, uniq_tiles.transpose( 1, 2, 0 ).reshape( int( uniq_tiles.size / 2 ), 2 ) ) ) tiles = {} for tile in uniq_tiles_list: tiles[tile] = self.get_tile_by_index( tile[1], tile[0] ) - + box = np.zeros( indices[0].shape, dtype=next(iter(tiles.values())).dtype ) # Now we have our tile set easily accessible based on the 1-to-1 mapping of indices to tile idx - + for tile in uniq_tiles_list: - if self.debug_: print( f"Processing tile {tile}") + if self._debug: + print( f"Processing tile {tile}") # Get the box indices corresponding to this tile box_indices = np.logical_and( as_tile_idx[0] == tile[0], as_tile_idx[1] == tile[1] ) - # Bulk assign to box this region of uniq tile by getting the relative index + # Bulk assign to box this region of uniq tile by getting the relative index # into that tile from original indices - box[box_indices] = tiles[tile][indices[0,box_indices] % self.tile_y_, indices[1,box_indices] % self.tile_x_] - - # Naively assign one-by-one - # for j in range( indices.shape[1] ): - # for i in range( indices.shape[2] ): - # reli = indices[1,j,i] % self.tile_x_ - # relj = indices[0,j,i] % self.tile_y_ - # box[j,i] = tiles[tuple(as_tile_idx[:,j,i])][relj,reli] + box[box_indices] = tiles[tile][indices[0, box_indices] % self._tile_y, indices[1, box_indices] % self._tile_x] - if self.debug_: self.print_present_tile_grid() + if self._debug: + self.print_present_tile_grid() return box diff --git a/util/compute_gwdo.py b/util/compute_gwdo.py index 6616c117..aa229b30 100755 --- a/util/compute_gwdo.py +++ b/util/compute_gwdo.py @@ -1,65 +1,167 @@ #!/usr/bin/env python3 - import sys +import os +import re +import argparse +from collections import OrderedDict -import numpy as np from netCDF4 import Dataset from SourceData import SourceData from OrographyStats import OrographyStats -topo_data = "/home/aislas/wrf-model/DATA/WPS_GEOG/topo_gmted2010_30s" -topo_source = SourceData( "gmted2010_30s", topo_data ) - -# For now just get cmdline arg -geo_files = sys.argv[1:] -box_size = 30000 - -for geo in geo_files: - print( f"Processing {geo}" ) - geo_data = Dataset( geo, "r+" ) - - xlat = geo_data.variables[ "XLAT_M" ] - xlon = geo_data.variables[ "XLONG_M" ] - - # Fields to overwrite - con = geo_data.variables[ "CON" ] - var = geo_data.variables[ "VAR" ] - - oa1 = geo_data.variables[ "OA1" ] - oa2 = geo_data.variables[ "OA2" ] - oa3 = geo_data.variables[ "OA3" ] - oa4 = geo_data.variables[ "OA4" ] - - ol1 = geo_data.variables[ "OL1" ] - ol2 = geo_data.variables[ "OL2" ] - ol3 = geo_data.variables[ "OL3" ] - ol4 = geo_data.variables[ "OL4" ] - - ns_size = xlat.shape[1] - we_size = xlat.shape[2] - - for j in range( ns_size ): - for i in range( we_size ): - lat = xlat[0,j,i] - lon = xlon[0,j,i] - box = topo_source.get_box( lat, lon, box_size ) - - oro_stats = OrographyStats( box ) - con[0,j,i] = oro_stats.con_ - var[0,j,i] = oro_stats.std_ - - oa1[0,j,i] = oro_stats.oa_[0] - oa2[0,j,i] = oro_stats.oa_[1] - oa3[0,j,i] = oro_stats.oa_[2] - oa4[0,j,i] = oro_stats.oa_[3] - - ol1[0,j,i] = oro_stats.ol_[0] - ol2[0,j,i] = oro_stats.ol_[1] - ol3[0,j,i] = oro_stats.ol_[2] - ol4[0,j,i] = oro_stats.ol_[3] - - geo_data.close() - -print( "Done!" ) +# https://stackoverflow.com/a/34325723 +def progerss_bar( iteration, total, prefix="", suffix="", precision=1, length=100, fill="*", empty="-", end="\r"): + """Print out a progress bar based on iteration/total of fill + iteration : current iteration + total : total iterations + prefix : prefix string + suffix : suffix string + decimals : precision of float in printing + length : character length of bar + fill : bar fill character + end : end character (e.g. "\r", "\r\n") + """ + progress = iteration / float( total ) + percent = ( f"{{0:.{precision}f}}").format( 100 * progress ) + current_length = int( length * progress ) + bar = fill * current_length + empty * ( length - current_length ) + print(f'\r{prefix} [{bar}] {percent}% {suffix}', end=end) + if iteration == total: + print( "" ) + + +def read_nml( filename ): + """Read a strict subset of Fortran 90 namelist + This ONLY reads namelists in the format: + &group + key0 = value0, + key1 = value0, value1, + key2='value0', 'value1', + key3 = "value0", "value1", + / + ...and so on... + Index slices, multiple key-value pairs on the same line, + and other features are NOT supported. For full support + please utilize f90nml. + Spaces between key and value and commas does not matter. + Multivalue entries MUST be delimited by commas. + """ + contents = "" + with open( filename, "r" ) as f: + contents = f.read() + + nml = OrderedDict() + for match in re.finditer( r"&(?P\w+)(?P.*?)\n*^[ ]*/", contents, re.S | re.M ): + group = match.group( "group" ) + nml[group] = OrderedDict() + for kv_match in re.finditer( r"^[ ]+(?P\w+)[ ]*=[ ]*(?P.*?)$", match.group("kv_pairs"), re.M ): + key = kv_match.group( "key" ) + values = list( + filter( + None, + [ + s.strip("\"' " ) + for s in list( filter( None, kv_match.group( "value" ).split(",") ) ) + ] + ) + ) + nml[group][key] = values if len( values ) > 1 else values[0] + return nml + + +def main(): + parser = argparse.ArgumentParser( + "./compute_gwdo.py", + description="Compute GWDO ancillary fields directly and update geogrid files", + formatter_class=argparse.ArgumentDefaultsHelpFormatter + ) + parser.add_argument( + "-n", "--namelist", + help="WPS namelist", + type=str, + default="./namelist.wps" + ) + parser.add_argument( + "-d", "--dataset", + help="Topographic dataset to use, in geogrid format", + type=str, + default="topo_gmted2010_30s" + ) + parser.add_argument( + "-n", "--no_progress", + help="Hide the progress bar", + action="store_true" + ) + options = parser.parse_args() + nml = read_nml( options.namelist ) + datapath= nml["geogrid"]["geog_data_path"] + + topo_data = os.path.join( datapath, options.dataset ) + # "/home/aislas/ext_data/wrf/DATA/WPS_GEOG/topo_gmted2010_30s" + topo_source = SourceData( options.dataset, topo_data ) + + # For now just get cmdline arg + max_dom = 1 + if "max_dom" in nml["share"]: + max_dom = int(nml["share"]["max_dom"]) + geo_files = [ f"geo_em.d{dom_id + 1:02d}.nc" for dom_id in range( max_dom ) ] + box_size_x = float( nml["geogrid"]["dx"] ) + box_size_y = float( nml["geogrid"]["dy"] ) + + for geo in geo_files: + print( f"Processing {geo}" ) + geo_data = Dataset( geo, "r+" ) + + xlat = geo_data.variables[ "XLAT_M" ] + xlon = geo_data.variables[ "XLONG_M" ] + + # Fields to overwrite + con = geo_data.variables[ "CON" ] + var = geo_data.variables[ "VAR" ] + + oa1 = geo_data.variables[ "OA1" ] + oa2 = geo_data.variables[ "OA2" ] + oa3 = geo_data.variables[ "OA3" ] + oa4 = geo_data.variables[ "OA4" ] + + ol1 = geo_data.variables[ "OL1" ] + ol2 = geo_data.variables[ "OL2" ] + ol3 = geo_data.variables[ "OL3" ] + ol4 = geo_data.variables[ "OL4" ] + + ns_size = xlat.shape[1] + we_size = xlat.shape[2] + + if not options.no_progress: + progerss_bar( 0, ns_size * we_size, length=50 ) + for j in range( ns_size ): + for i in range( we_size ): + lat = xlat[0, j, i] + lon = xlon[0, j, i] + box = topo_source.get_box( lat, lon, box_size_x, box_size_y ) + + oro_stats = OrographyStats( box ) + con[0, j, i] = oro_stats.con + var[0, j, i] = oro_stats.std + + oa1[0, j, i] = oro_stats.oa[0] + oa2[0, j, i] = oro_stats.oa[1] + oa3[0, j, i] = oro_stats.oa[2] + oa4[0, j, i] = oro_stats.oa[3] + + ol1[0, j, i] = oro_stats.ol[0] + ol2[0, j, i] = oro_stats.ol[1] + ol3[0, j, i] = oro_stats.ol[2] + ol4[0, j, i] = oro_stats.ol[3] + if not options.no_progress: + progerss_bar( i + j * we_size + 1, ns_size * we_size, length=50 ) + + geo_data.close() + + print( "Done!" ) + + +if __name__ == "__main__": + main() From 8bcaae0dd5395eacc486b428e37c2e1ba9aef9e1 Mon Sep 17 00:00:00 2001 From: Anthony Islas Date: Mon, 13 Oct 2025 16:10:48 -0700 Subject: [PATCH 16/20] Add MAX_EL --- util/compute_gwdo.py | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/util/compute_gwdo.py b/util/compute_gwdo.py index aa229b30..e7e27b4e 100755 --- a/util/compute_gwdo.py +++ b/util/compute_gwdo.py @@ -1,11 +1,10 @@ #!/usr/bin/env python3 -import sys import os import re import argparse from collections import OrderedDict -from netCDF4 import Dataset +import netCDF4 as nc4 from SourceData import SourceData from OrographyStats import OrographyStats @@ -90,7 +89,7 @@ def main(): default="topo_gmted2010_30s" ) parser.add_argument( - "-n", "--no_progress", + "-hp", "--hide_progress", help="Hide the progress bar", action="store_true" ) @@ -112,7 +111,7 @@ def main(): for geo in geo_files: print( f"Processing {geo}" ) - geo_data = Dataset( geo, "r+" ) + geo_data = nc4.Dataset( geo, "r+" ) xlat = geo_data.variables[ "XLAT_M" ] xlon = geo_data.variables[ "XLONG_M" ] @@ -131,10 +130,15 @@ def main(): ol3 = geo_data.variables[ "OL3" ] ol4 = geo_data.variables[ "OL4" ] + if "MAX_EL" in geo_data.variables: + max_el = geo_data.variables["MAX_EL"] + else: + max_el = geo_data.createVariable( "MAX_EL", "f4", ( var.dimensions ) ) + ns_size = xlat.shape[1] we_size = xlat.shape[2] - if not options.no_progress: + if not options.hide_progress: progerss_bar( 0, ns_size * we_size, length=50 ) for j in range( ns_size ): for i in range( we_size ): @@ -155,7 +159,9 @@ def main(): ol2[0, j, i] = oro_stats.ol[1] ol3[0, j, i] = oro_stats.ol[2] ol4[0, j, i] = oro_stats.ol[3] - if not options.no_progress: + + max_el[0, j, i] = oro_stats.max + if not options.hide_progress: progerss_bar( i + j * we_size + 1, ns_size * we_size, length=50 ) geo_data.close() From 82e1dd43a2a1baa4897a343d00ba9368766eacaa Mon Sep 17 00:00:00 2001 From: Anthony Islas Date: Mon, 8 Dec 2025 12:20:50 -0700 Subject: [PATCH 17/20] Remove reuse of mean for older numpy versions --- util/OrographyStats.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/util/OrographyStats.py b/util/OrographyStats.py index 8402c615..221ef992 100644 --- a/util/OrographyStats.py +++ b/util/OrographyStats.py @@ -7,7 +7,7 @@ def __init__( self, box ): self.mean = np.mean( self.box ) # var (actulally stddev) - self.std = np.std( self.box, mean=self.mean ) + self.std = np.std( self.box ) self.max = np.max( box ) # Critical height used in calculation of orographic effective length From bb555e1bc67f484bdd2bed2c677a449ecdaad682 Mon Sep 17 00:00:00 2001 From: Anthony Islas Date: Fri, 23 Jan 2026 16:22:43 -0700 Subject: [PATCH 18/20] Update calcs to account for odd sized boxes and use mean for orographic length --- util/OrographyStats.py | 38 ++++++++++++++++++++++---------------- util/compute_gwdo.py | 4 ++-- 2 files changed, 24 insertions(+), 18 deletions(-) diff --git a/util/OrographyStats.py b/util/OrographyStats.py index 221ef992..1c1aa67c 100644 --- a/util/OrographyStats.py +++ b/util/OrographyStats.py @@ -10,9 +10,6 @@ def __init__( self, box ): self.std = np.std( self.box ) self.max = np.max( box ) - # Critical height used in calculation of orographic effective length - self.hc = 1116.2 - 0.878 * self.std - # This currently does not use the landuse to average over only land and zero # out on mostly water # con (convexity) @@ -38,39 +35,48 @@ def calc_oa( self ): # with [y,x] indices # oa1 is the orographic asymmetry in the West direction - nu = np.sum( self.box[:, :int(self.box.shape[1] / 2)] > self.mean ) - nd = np.sum( self.box[:, int(self.box.shape[1] / 2):] > self.mean ) + nu = np.sum( self.box[:, :int((self.box.shape[1] + self.box.shape[1] % 2) / 2)] > self.mean ) + nd = np.sum( self.box[:, int((self.box.shape[1] - self.box.shape[1] % 2) / 2):] > self.mean ) self.oa[0] = ( nu - nd ) / ( nu + nd ) if ( ( nu + nd ) > 0 ) else 0.0 # oa2 is the orographic asymmetry in the South direction - nu = np.sum( self.box[:int(self.box.shape[0] / 2), :] > self.mean ) - nd = np.sum( self.box[int(self.box.shape[0] / 2):, :] > self.mean ) + nu = np.sum( self.box[:int((self.box.shape[0] + self.box.shape[0] % 2) / 2), :] > self.mean ) + nd = np.sum( self.box[int((self.box.shape[0] - self.box.shape[0] % 2) / 2):, :] > self.mean ) self.oa[1] = ( nu - nd ) / ( nu + nd ) if ( ( nu + nd ) > 0 ) else 0.0 # Pre-compute the geometric diagonal of the box - slope = self.box.shape[1] / self.box.shape[0] + slope = self.box.shape[0] / self.box.shape[1] j, i = np.indices( self.box.shape ) - upstream = np.flip( i <= ( j * slope ), axis=0 ) + # Corrected - all indices that lie on diagonal are counted in both upstream and downstream + vals = ( i + 1 ) * slope - ( self.box.shape[0] - j ) + upstream = ( vals <= slope ) + downstream = ( vals >= -1.0 ) + + # MPAS calcs - slightly accounts for diagonal when i ~= j but increasingly off for larger discrepancies + # vals = np.rint( ( i + 1 ) * slope ) - ( self.box.shape[0] - j ) + # upstream = ( vals <= 0 ) + # downstream = ( vals >= 0 ) # oa3 is the orographic asymmetry in the South-West direction nu = np.sum( self.box[upstream] > self.mean ) - nd = np.sum( self.box[~upstream] > self.mean ) + nd = np.sum( self.box[downstream] > self.mean ) self.oa[2] = ( nu - nd ) / ( nu + nd ) if ( ( nu + nd ) > 0 ) else 0.0 # oa4 is the orographic asymmetry in the North-West direction - upstream = np.flip( upstream, axis=0 ) + upstream = np.flip( upstream, axis=1 ) + downstream = np.flip( downstream, axis=1 ) nu = np.sum( self.box[upstream] > self.mean ) - nd = np.sum( self.box[~upstream] > self.mean ) + nd = np.sum( self.box[downstream] > self.mean ) self.oa[3] = ( nu - nd ) / ( nu + nd ) if ( ( nu + nd ) > 0 ) else 0.0 def calc_ol( self ): # ol1 is the effective orographic length in the West direction interior = self.box[int(np.floor(self.box.shape[0] * .25)):int(np.ceil(self.box.shape[0] * .75)), :] - self.ol[0] = np.sum( interior > self.hc ) / interior.size + self.ol[0] = np.sum( interior > self.mean ) / interior.size # ol2 is the effective orographic length in the South direction interior = self.box[:, int(np.floor(self.box.shape[1] * .25)):int(np.ceil(self.box.shape[1] * .75))] - self.ol[1] = np.sum( interior > self.hc ) / interior.size + self.ol[1] = np.sum( interior > self.mean ) / interior.size # The prescribed methodology uses 4 quadrants to get the diagonals and effectively # test half of the box, however this does not actually test the interior half @@ -80,10 +86,10 @@ def calc_ol( self ): interiorA = self.box[:int(self.box.shape[0] / 2), :int(self.box.shape[1] / 2)] # first half of x first half of y interiorB = self.box[int(self.box.shape[0] / 2):, int(self.box.shape[1] / 2):] # second half of x second half of y - self.ol[2] = ( np.sum( interiorA > self.hc ) + np.sum( interiorB > self.hc ) ) / ( interiorA.size + interiorB.size ) + self.ol[2] = ( np.sum( interiorA > self.mean ) + np.sum( interiorB > self.mean ) ) / ( interiorA.size + interiorB.size ) # ol4 is the effective orographic length in the North-West direction interiorA = self.box[int(self.box.shape[0] / 2):, :int(self.box.shape[1] / 2)] # first half of x second half of y interiorB = self.box[:int(self.box.shape[0] / 2), int(self.box.shape[1] / 2):] # second half of x first half of y - self.ol[3] = ( np.sum( interiorA > self.hc ) + np.sum( interiorB > self.hc ) ) / ( interiorA.size + interiorB.size ) + self.ol[3] = ( np.sum( interiorA > self.mean ) + np.sum( interiorB > self.mean ) ) / ( interiorA.size + interiorB.size ) diff --git a/util/compute_gwdo.py b/util/compute_gwdo.py index e7e27b4e..4503b1dc 100755 --- a/util/compute_gwdo.py +++ b/util/compute_gwdo.py @@ -106,8 +106,8 @@ def main(): if "max_dom" in nml["share"]: max_dom = int(nml["share"]["max_dom"]) geo_files = [ f"geo_em.d{dom_id + 1:02d}.nc" for dom_id in range( max_dom ) ] - box_size_x = float( nml["geogrid"]["dx"] ) - box_size_y = float( nml["geogrid"]["dy"] ) + box_size_x = float( nml["geogrid"]["dx"] ) * 2 + box_size_y = float( nml["geogrid"]["dy"] ) * 2 for geo in geo_files: print( f"Processing {geo}" ) From 7a1bbd3b76d76f6bf8b2e57d443c07859a9aa5c7 Mon Sep 17 00:00:00 2001 From: Anthony Islas Date: Mon, 26 Jan 2026 14:27:05 -0700 Subject: [PATCH 19/20] Fix axis flip for oa4 --- util/OrographyStats.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/util/OrographyStats.py b/util/OrographyStats.py index 1c1aa67c..d55481cf 100644 --- a/util/OrographyStats.py +++ b/util/OrographyStats.py @@ -63,8 +63,8 @@ def calc_oa( self ): self.oa[2] = ( nu - nd ) / ( nu + nd ) if ( ( nu + nd ) > 0 ) else 0.0 # oa4 is the orographic asymmetry in the North-West direction - upstream = np.flip( upstream, axis=1 ) - downstream = np.flip( downstream, axis=1 ) + upstream = np.flip( upstream, axis=0 ) + downstream = np.flip( downstream, axis=0 ) nu = np.sum( self.box[upstream] > self.mean ) nd = np.sum( self.box[downstream] > self.mean ) self.oa[3] = ( nu - nd ) / ( nu + nd ) if ( ( nu + nd ) > 0 ) else 0.0 From 0eacdc042e7645be7ee0643d2a3047bf737be65f Mon Sep 17 00:00:00 2001 From: Anthony Islas Date: Wed, 4 Feb 2026 16:04:06 -0700 Subject: [PATCH 20/20] Code cleanup and revert unused changes --- util/CMakeLists.txt | 29 +--- util/SourceData.py | 29 ---- util/compute_gwdo.py | 1 - util/src/CMakeLists.txt | 30 ---- util/src/compute_gwdo.F | 374 ---------------------------------------- 5 files changed, 1 insertion(+), 462 deletions(-) delete mode 100644 util/src/compute_gwdo.F diff --git a/util/CMakeLists.txt b/util/CMakeLists.txt index ba2cbbae..684ecae3 100644 --- a/util/CMakeLists.txt +++ b/util/CMakeLists.txt @@ -11,7 +11,6 @@ add_executable( calc_ecmwf_p ) add_executable( mod_levs ) if ( DEFINED WRF_DIR ) add_executable( int2nc ) - add_executable( compute_gwdo ) endif() add_executable( height_ukmo ) @@ -57,27 +56,6 @@ if ( DEFINED WRF_DIR ) PRIVATE ${netCDF_INCLUDE_DIRS} ) - target_link_libraries( - compute_gwdo - PRIVATE - WRF::WRF_Core - netCDF::netcdff - ) - set_target_properties( - compute_gwdo - PROPERTIES - # Just dump everything in here - Fortran_MODULE_DIRECTORY ${CMAKE_INSTALL_PREFIX}/modules/util/compute_gwdo - Fortran_FORMAT FREE - ) - install( - TARGETS compute_gwdo - EXPORT ${EXPORT_NAME}Targets - RUNTIME DESTINATION bin/ - ARCHIVE DESTINATION lib/ - LIBRARY DESTINATION lib/ - ) - endif() # Every single target @@ -88,7 +66,7 @@ foreach( TARGET ${ALL_UTIL_TARGETS} util_common ) ${TARGET} PROPERTIES # Just dump everything in here - Fortran_MODULE_DIRECTORY ${CMAKE_INSTALL_PREFIX}/modules/util/${TARGET}/ + Fortran_MODULE_DIRECTORY ${CMAKE_INSTALL_PREFIX}/modules/util/ Fortran_FORMAT FREE ) @@ -99,11 +77,6 @@ foreach( TARGET ${ALL_UTIL_TARGETS} util_common ) PRIVATE util_common ) - target_include_directories( - ${TARGET} - PRIVATE - $ - ) endif() # Add these to the export targets diff --git a/util/SourceData.py b/util/SourceData.py index 8f4a38b2..81ffb3b5 100644 --- a/util/SourceData.py +++ b/util/SourceData.py @@ -114,34 +114,6 @@ def __init__( self, name, path ): self._pts_per_deg = int( 1.0 / self._index.dx ) self._subgrid_m_dx = 2.0 * np.pi * self._earth_radius / self._npts_x - # # These are all probably wrong - # if self._index.projection_ == "lambert": - # self._projection = cartopy.crs.LambertConformal( - # central_longitude =self._index.known_lon_, - # central_latitude =self._index.known_lat_, - # standard_parallels =( self._index.truelat1_, self._index.truelat2_ ) - # ) - # elif self._index.projection_ == "polar_wgs84" or self._index.projection_ == "polar": - # self._projection = cartopy.crs.Stereographic( - # central_longitude =self._index.known_lon_, - # central_latitude =self._index.known_lat_, - # true_scale_latitude =self._index.truelat1_ - # ) - # elif self._index.projection_ == "albers_nad83": - # self._projection = cartopy.crs.AlbersEqualArea( - # central_longitude =self._index.known_lon_, - # central_latitude =self._index.known_lat_, - # standard_parallels =( self._index.truelat1_, self._index.truelat2_ ) - # ) - # elif self._index.projection_ == "mercator": - # self._projection = cartopy.crs.Mercator( - # central_longitude =self._index.known_lon_, - # # central_latitude =self._index.known_lat_, - # latitude_true_scale =self._index.truelat1_ - # ) - # elif self._index.projection_ == "regular_ll": - # self._projection = cartopy.crs.LambertCylindrical( central_longitude =self._index.known_lon_ ) - self._tile_data = TileData( int( int( 360.0 / self._index.dx ) / self._index.tile_x ), int( int( 180.0 / self._index.dy ) / self._index.tile_y ), @@ -179,7 +151,6 @@ def read_geogrid( self, file ): (increasing latitude) NE """ - # print( "Reading " + file ) rawdata = np.fromfile( file, dtype=self._index.get_dtype() diff --git a/util/compute_gwdo.py b/util/compute_gwdo.py index 4503b1dc..05a8f239 100755 --- a/util/compute_gwdo.py +++ b/util/compute_gwdo.py @@ -98,7 +98,6 @@ def main(): datapath= nml["geogrid"]["geog_data_path"] topo_data = os.path.join( datapath, options.dataset ) - # "/home/aislas/ext_data/wrf/DATA/WPS_GEOG/topo_gmted2010_30s" topo_source = SourceData( options.dataset, topo_data ) # For now just get cmdline arg diff --git a/util/src/CMakeLists.txt b/util/src/CMakeLists.txt index 691aeff8..18e42d8c 100644 --- a/util/src/CMakeLists.txt +++ b/util/src/CMakeLists.txt @@ -70,33 +70,3 @@ target_sources( PRIVATE height_ukmo.F ) - -target_sources( - compute_gwdo - PRIVATE - compute_gwdo.F - # module_input_module.F90 - ${PROJECT_SOURCE_DIR}/geogrid/src/cio.c - ${PROJECT_SOURCE_DIR}/geogrid/src/wrf_debug.F - # ${PROJECT_SOURCE_DIR}/geogrid/src/bitarray_module.F - ${PROJECT_SOURCE_DIR}/geogrid/src/constants_module.F - # ${PROJECT_SOURCE_DIR}/geogrid/src/module_stringutil.F - # ${PROJECT_SOURCE_DIR}/geogrid/src/geogrid.F - ${PROJECT_SOURCE_DIR}/geogrid/src/gridinfo_module.F - # ${PROJECT_SOURCE_DIR}/geogrid/src/hash_module.F - # ${PROJECT_SOURCE_DIR}/geogrid/src/interp_module.F - ${PROJECT_SOURCE_DIR}/geogrid/src/list_module.F - ${PROJECT_SOURCE_DIR}/geogrid/src/llxy_module.F - ${PROJECT_SOURCE_DIR}/geogrid/src/misc_definitions_module.F - ${PROJECT_SOURCE_DIR}/geogrid/src/module_debug.F - ${PROJECT_SOURCE_DIR}/geogrid/src/module_map_utils.F - # ${PROJECT_SOURCE_DIR}/geogrid/src/module_map_utils.F - # ${PROJECT_SOURCE_DIR}/geogrid/src/output_module.F - # ${PROJECT_SOURCE_DIR}/geogrid/src/parallel_module.F - # ${PROJECT_SOURCE_DIR}/geogrid/src/process_tile_module.F - # ${PROJECT_SOURCE_DIR}/geogrid/src/proc_point_module.F - # ${PROJECT_SOURCE_DIR}/geogrid/src/queue_module.F - # ${PROJECT_SOURCE_DIR}/geogrid/src/read_geogrid.c - # ${PROJECT_SOURCE_DIR}/geogrid/src/smooth_module.F - # ${PROJECT_SOURCE_DIR}/geogrid/src/source_data_module.F - ) diff --git a/util/src/compute_gwdo.F b/util/src/compute_gwdo.F deleted file mode 100644 index 7de2bd24..00000000 --- a/util/src/compute_gwdo.F +++ /dev/null @@ -1,374 +0,0 @@ -module nc_helpers -contains - subroutine nf90_check_error( stat ) - use netcdf - implicit none - integer, intent( in ) :: stat - if ( stat /= nf90_noerr ) then - ! We had an error - write( 0, * ) __FILE__, ":", __LINE__, " : ", trim( nf90_strerror( stat ) ) - stop 1 - end if - end subroutine nf90_check_error -end module nc_helpers - -module source_data_module - type :: source_index_raw - character (len=256):: source_path = "required" - character (len=32) :: projection = "required" - character (len=32) :: source_type = "required" - logical :: signed = .false. - character (len=32) :: units = "required" - character (len=64) :: description = "required" - real :: dx = -1.0 !< required - real :: dy = -1.0 !< required - real :: known_x = 1.0 - real :: known_y = 1.0 - real :: known_lat = -1.0 !< required - real :: known_lon = -1.0 !< required - real :: stdlon = -2.0 !< not required - real :: truelat1 = -2.0 !< not required - real :: truelat2 = -2.0 !< not required - integer :: wordsize = -1 !< required - integer :: tile_x = -1 !< required - integer :: tile_y = -1 !< required - integer :: tile_z = -2 !< not required - integer :: tile_z_start = -2 !< not required - integer :: tile_z_end = -2 !< not required - integer :: category_min = -2 !< not required - integer :: category_max = -2 !< not required - integer :: tile_bdr = 0 - real :: missing_value = -2.0 !< not required - real :: scale_factor = 1.0 - character (len=32) :: row_order = "bottom_top" - character (len=32) :: endian = "big" - integer :: iswater = 16 - integer :: islake = -1 !< i.e. no separate inland water category (does not mean required) - integer :: isice = 24 - integer :: isurban = 1 - integer :: isoilwater = 14 - character (len=32) :: mminlu = "USGS" - integer :: filename_digits = 5 - end type source_index_raw - -contains - - subroutine read_index_file( source_path, index_data ) - use, intrinsic :: iso_fortran_env - - implicit none - - character (len=*), intent( in ) :: source_path - character (len=256) :: index_path, line_buffer, key, var - logical :: is_used, kv_exist - integer :: iostatus, funit - - type( source_index_raw ), intent( inout ) :: index_data - - write( index_data % source_path, "(a)" ) trim(source_path) - write( index_path, "(a)" ) trim(source_path) // 'index' - write( *, * ) "Reading 'index' file at ", index_data % source_path - - ! - ! Open the index file for the data source for this field, and read in metadata specs - ! - do funit=10,100 - inquire(unit=funit, opened=is_used) - if (.not. is_used) exit - end do - open(funit,file=trim(index_path),form='formatted',status='old',iostat=iostatus) - - ! TODO add optional to control when not found - if (iostatus /= 0 ) then - write( *, * ) 'Could not open ', trim(index_data % source_path) - end if - - do - read( funit, '(a)', iostat=iostatus ) line_buffer - if ( iostatus /= 0 ) then - if ( iostatus == iostat_end ) exit - end if - - ! Process line iff a key-value pair was found - kv_exist = get_kv_pair( line_buffer, key, var ) - if ( kv_exist ) then - ! consider multiline else if even if not neater - if ( key == "projection" ) index_data % projection = var - if ( key == "source_type" ) index_data % source_type = var - if ( key == "signed" .and. var == "yes" ) index_data % signed = .true. - if ( key == "units" ) index_data % units = var - if ( key == "description" ) index_data % description = var - if ( key == "dx" ) read( var, * ) index_data % dx - if ( key == "dy" ) read( var, * ) index_data % dy - if ( key == "known_x" ) read( var, * ) index_data % known_x - if ( key == "known_y" ) read( var, * ) index_data % known_y - if ( key == "known_lat" ) read( var, * ) index_data % known_lat - if ( key == "known_lon" ) read( var, * ) index_data % known_lon - if ( key == "stdlon" ) read( var, * ) index_data % stdlon - if ( key == "truelat1" ) read( var, * ) index_data % truelat1 - if ( key == "truelat2" ) read( var, * ) index_data % truelat2 - if ( key == "wordsize" ) read( var, * ) index_data % wordsize - if ( key == "tile_x" ) read( var, * ) index_data % tile_x - if ( key == "tile_y" ) read( var, * ) index_data % tile_y - if ( key == "tile_z" ) read( var, * ) index_data % tile_z - if ( key == "tile_z_start" ) read( var, * ) index_data % tile_z_start - if ( key == "tile_z_end" ) read( var, * ) index_data % tile_z_end - if ( key == "category_min" ) read( var, * ) index_data % category_min - if ( key == "category_max" ) read( var, * ) index_data % category_max - if ( key == "tile_bdr" ) read( var, * ) index_data % tile_bdr - if ( key == "missing_value" ) read( var, * ) index_data % missing_value - if ( key == "scale_factor" ) read( var, * ) index_data % scale_factor - if ( key == "row_order" ) index_data % row_order = var - if ( key == "endian" ) index_data % endian = var - if ( key == "iswater" ) read( var, * ) index_data % iswater - if ( key == "islake" ) read( var, * ) index_data % islake - if ( key == "isice" ) read( var, * ) index_data % isice - if ( key == "isurban" ) read( var, * ) index_data % isurban - if ( key == "isoilwater" ) read( var, * ) index_data % isoilwater - if ( key == "mminlu" ) index_data % mminlu = var - if ( key == "filename_digits" ) read( var, * ) index_data % filename_digits - end if - end do - end subroutine - - logical function get_kv_pair( line, key, var ) result( found ) - implicit none - character (len=*), intent( in ) :: line - character (len=*), intent( inout ) :: key, var - integer :: idx, idx_comment - ! logical, intent( out ) :: found - - found = .false. - - idx = index( line, "=" ) - idx_comment = index( line, "#" ) - - if ( idx == 0 ) return - - if ( idx_comment /= 0 .and. idx_comment < idx ) then - ! We have a comment before our = and that is not good - write( *, * ) "Key-value parse error, comment before = " - return - end if - - if ( idx_comment == 0 ) idx_comment = len( line ) - - ! We are fine to parse - key = trim(adjustl(line(1:idx-1))) - var = trim(adjustl(line(idx+1:idx_comment))) - found = .true. - - end function get_kv_pair - - integer function get_proj_code( index_data ) result( proj_code ) - use misc_definitions_module - implicit none - type( source_index_raw ), intent( in ) :: index_data - if ( index_data % projection == "lambert" ) then; proj_code = PROJ_LC - else if ( index_data % projection == "polar_wgs84" ) then; proj_code = PROJ_PS_WGS84 - else if ( index_data % projection == "albers_nad83" ) then; proj_code = PROJ_ALBERS_NAD83 - else if ( index_data % projection == "polar" ) then; proj_code = PROJ_PS - else if ( index_data % projection == "mercator" ) then; proj_code = PROJ_MERC - else if ( index_data % projection == "regular_ll" ) then; proj_code = PROJ_LATLON - else; write( *, * ) "ERROR: Unknown projection '", index_data % projection, "'" - end if - end function get_proj_code - - - ! Taken from the souce_data_module of geogrid, but simplified to take stateful index data - subroutine get_tile_fname(index_data, lat, lon, tile_fname, istatus) - - use llxy_module - use gridinfo_module - - implicit none - - type( source_index_raw ), intent( in ) :: index_data - real, intent(in) :: lat, lon - character (len=*), intent(out) :: tile_fname - integer, intent(out) :: istatus - - ! Local variables - integer :: current_domain, idx - real :: rx, ry - - istatus = 0 - write(tile_fname, '(a)') 'TILE.OUTSIDE.DOMAIN' - - ! do idx=1,num_entries - ! if ((index(source_fieldname(idx),trim(fieldname)) /= 0) .and. & - ! (len_trim(source_fieldname(idx)) == len_trim(fieldname))) then - ! if (ilevel == source_priority(idx)) then - ! istatus = 0 - ! exit - ! end if - ! end if - ! end do - - ! if (istatus /= 0) return - - current_domain = iget_selected_domain() - call select_domain(SOURCE_PROJ) - call lltoxy(lat, lon, rx, ry, M) - call select_domain(current_domain) - - ! rx = real(source_tile_x(idx)) * real(floor((rx-0.5*source_dx(idx))/ real(source_tile_x(idx)))) + 1.0 - ! ry = real(source_tile_y(idx)) * real(floor((ry-0.5*source_dy(idx))/ real(source_tile_y(idx)))) + 1.0 - rx = real( index_data % tile_x ) * real(floor((rx-0.5) / real( index_data % tile_x ))) + 1.0 - ry = real( index_data % tile_y ) * real(floor((ry-0.5) / real( index_data % tile_y ))) + 1.0 - - if (rx > 0. .and. ry > 0.) then - select case(index_data % filename_digits) - case (5) - write(tile_fname, '(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim( index_data % source_path ), & - nint(rx),'-',nint(rx)+index_data % tile_x -1,'.',nint(ry),'-',nint(ry)+index_data % tile_y -1 - case (6) - write(tile_fname, '(a,i6.6,a1,i6.6,a1,i6.6,a1,i6.6)') trim( index_data % source_path ), & - nint(rx),'-',nint(rx)+index_data % tile_x -1,'.',nint(ry),'-',nint(ry)+index_data % tile_y -1 - case default - call mprintf(.true., ERROR, 'In GEOGRID.TBL, for index file of data at '// & - 'entry %i, filename_digits must be either 5 or 6.', i1=idx) - istatus = 1 - return - end select - end if - - end subroutine get_tile_fname -end module source_data_module - -program compute_gwdo - use netcdf - use nc_helpers - use source_data_module - use gridinfo_module - use llxy_module - - implicit none - character (len=128), parameter :: geo_em_prefix = "geo_em.d" - character (len=128) :: geo_em_file, tmp - integer :: i, j, k, l - integer :: ierr, ncid, ndims, dim1, dim2, dim3 - integer, dimension(nf90_max_var_dims) :: dimids - integer :: start_i, end_i, start_j, end_j, start_k, end_k - integer :: mapfacxid, mapfacyid - integer :: latvarid, lonvarid - integer :: convarid, varvarid - integer, dimension(4) :: olvarid, oavarid - - ! I will probably need more - real, dimension(:,:,:), allocatable :: con_data - real, dimension(:,:,:), allocatable :: lats, lons - - real :: cur_lat, cur_lon - type( source_index_raw ) :: index_data - - - - - call get_grid_params() - - write( *, * ) "Found geogrid : ", n_domains - write( *, * ) "Using data from : ", geog_data_path - - write( tmp, "(a)" ) trim(geog_data_path) // "topo_gmted2010_30s/" - call read_index_file( tmp, index_data ) - - ! For now set the source data projection here - call push_source_projection( get_proj_code( index_data ), index_data % stdlon, & - index_data % truelat1, index_data % truelat2, & - index_data % dx, index_data % dy, & - index_data % dy, index_data % dx, & - index_data % known_x, index_data % known_y, & - index_data % known_lat, index_data % known_lon ) - - - - write( *, * ) "Source data dx / dy : ", index_data % dx, ", ", index_data % dy - - do i = 1, n_domains - write( geo_em_file, "(A8,I0.2,A3)" ) trim( geo_em_prefix ), i, ".nc" - write( *, * ) "Processing ", geo_em_file - ierr = nf90_open( geo_em_file, NF90_WRITE, ncid ) - call nf90_check_error( ierr ) - - ierr = nf90_inq_varid( ncid, "XLAT_M", latvarid ) - call nf90_check_error( ierr ) - - ierr = nf90_inq_varid( ncid, "XLONG_M", lonvarid ) - call nf90_check_error( ierr ) - - ierr = nf90_inq_varid( ncid, "MAPFAC_MX", mapfacxid ) - call nf90_check_error( ierr ) - - ierr = nf90_inq_varid( ncid, "MAPFAC_MY", mapfacyid ) - call nf90_check_error( ierr ) - - ierr = nf90_inquire_variable( ncid, lonvarid, ndims = ndims, dimids = dimids ) - call nf90_check_error( ierr ) - - write( *, * ) "Variable CLONG has ", ndims, " dims" - - ierr = nf90_inq_varid( ncid, "CON", convarid ) - call nf90_check_error( ierr ) - - ierr = nf90_inq_varid( ncid, "VAR", varvarid ) - call nf90_check_error( ierr ) - - do j = 1, 4 - write( tmp, "(A2,I1)" ) "OL", j - ierr = nf90_inq_varid( ncid, tmp, olvarid(j) ) - write( tmp, "(A2,I1)" ) "OA", j - ierr = nf90_inq_varid( ncid, tmp, oavarid(j) ) - end do - - ierr = nf90_inquire_dimension( ncid, dimids(1), len = dim1 ) - ierr = nf90_inquire_dimension( ncid, dimids(2), len = dim2 ) - ierr = nf90_inquire_dimension( ncid, dimids(3), len = dim3 ) - - write( *, * ) "Dimensions sizes are ", dim1, ", ", dim2, ", ", dim3 - - allocate( con_data( dim1, dim2, dim3 ) ) - allocate( lats( dim1, dim2, dim3 ) ) - allocate( lons( dim1, dim2, dim3 ) ) - - ierr = nf90_get_var( ncid, latvarid, lats ) - ierr = nf90_get_var( ncid, lonvarid, lons ) - - ! Find correct file and location within file to grab - - - ! Just stuff for now - do j = 1, dim1 - do k = 1, dim2 - do l = 1, dim3 - con_data(j,k,l) = 3.1415926535 - end do - end do - end do - - ! Just stuff for now - do j = 1, dim1 - do k = 1, dim2 - do l = 1, dim3 - cur_lat = lats(j,k,l) - cur_lon = lons(j,k,l) - call get_tile_fname( index_data, cur_lat, cur_lon, tmp, ierr ) - if ( j == 10 .and. k == 10 ) then - write( *, * ) "Filename for [", i, ", ", j, ", ", k, "] at LL [", cur_lat, ", ", cur_lon, "] :", tmp - end if - end do - end do - end do - - ierr = nf90_put_var( ncid, convarid, con_data ) - call nf90_check_error( ierr ) - - - ierr = nf90_close( ncid ) - call nf90_check_error( ierr ) - - end do - - return -end program compute_gwdo -