Develop fix vegetation geodata (#943)

* Update gridlabd-geodata

* Update gridlabd-geodata

* Update gridlabd-geodata

* Update gridlabd-geodata

* Update gridlabd-geodata

* Update gridlabd-geodata

* Update geodata_example.glm

* Update geodata_example

* Work update

* Update geodata_address.py

* Update gridlabd-geodata

* Fix traceback output when DEBUG is True

* Add general geodata autotest

* Update gridlabd-geodata

* Update test_geodata_address.glm

* Update gridlabd-geodata

* Add geodata dataset autotest

* Update geodata_utility.py

* Update autotest.sh

* Rename unittest.sh

* Update geodata_elevation.py

* Update geodata_address.py

* Update geodata_address.py

* Update geodata_elevation.py

* Update geodata_utility.py

* Update test_geodata_elevation.glm

* Update test_geodata_utility.glm

* Update unittest.sh

* Update gridlabd-geodata

* Update test_geodata.txt

* Update unittest.sh

* Update geodata_powerline.py

* Fix index resolution conflicts

* Update geodata_powerline.py

* Add haversine function support

* Update requirements.txt

* Update gridlabd-geodata

* Update gridlabd-geodata

* Update geodata_powerline.py

* Update gridlabd-geodata

* Update test_geodata.glm

* Update gridlabd-geodata

* Update gridlabd-geodata

* Update main.cpp

* Update test_geodata.glm

* Update geodata package call method

* Update geodata_powerline.py

* Update geodata_powerline.py

* Update gridlabd-geodata

* Update geodata_census.py

* Add cabletypes for powerline geodata package

* Update gridlabd-geodata

* Update geodata_powerline.py

* Update geodata_powerline.py

* Update geodata_powerline.py

* Update geodata_powerline_cabletypes.csv

* Update geodata_powerline.py

* Update geodata_powerline.py

* Add csv-table to glm-library converter

* Update geodata_powerline.py

* Refactor gridlabd convert implementation

* Update Convert.md

* Update Convert.md

* Update gridlabd-convert

* Update csv-table2glm-library.py

* Update gridlabd-convert

* Update csv-table2glm-library.py

* Update geodata_powerline.py

* Update geodata_powerline.py

* Update path_example.csv

* Update geodata_powerline.py

* Update test_geodata.txt

* Add verbose output for get_args()

* Add census state data processing

* Update geodata_census.py

* Add more support for plot format

* Add support for uuid key

* Update gridlabd-geodata

* Refactor python requirements handling

* Update gridlabd-geodata

* Update geodata_census.py

* Fix redirect issue

* Update test_geodata.glm

* Fix permissions on scripts

* Update Makefile.am

* Update gridlabd-geodata

* Update geodata_census.py

* Update geodata_census.py

* Update gridlabd-geodata

* Added zipcode-level census data

* Update zipcodes.py

* Update gridlabd-geodata

* Update geodata_census.py

* Create path_example.csv

* Create GridLAB-D GeoData Tutorial.ipynb

* Update geodata_distance.py

* Update build_number

* Update gridlabd-geodata

* Update GridLAB-D GeoData Tutorial.ipynb

* Update geodata_elevation.py

* Update gridlabd-geodata

* Update geodata_utility.py

* Update GridLAB-D GeoData Tutorial.ipynb

* Update geodata_utility.py

* Update gridlabd-geodata

* Update GridLAB-D GeoData Tutorial.ipynb

* Update gridlabd-geodata

* Update geodata_powerline.py

* Update GridLAB-D GeoData Tutorial.ipynb

* Update geodata_powerline.py

* Update GridLAB-D GeoData Tutorial.ipynb

* Update GridLAB-D GeoData Tutorial.ipynb

* Update gridlabd-geodata

* Update GridLAB-D GeoData Tutorial.ipynb

* Update gridlabd-geodata

* Update geodata_distance.py

* Update GridLAB-D GeoData Tutorial.ipynb

* Update gridlabd-geodata

* Update GridLAB-D GeoData Tutorial.ipynb

* Update GridLAB-D GeoData Tutorial.ipynb

* Update GridLAB-D GeoData Tutorial.ipynb

* Update gridlabd-geodata

* Update GridLAB-D GeoData Tutorial.ipynb

* Update gridlabd-geodata

* Update gridlabd-geodata

* Update geodata_distance.py

* Update Makefile.am

* Update gridlabd-geodata

* Update GridLAB-D GeoData Tutorial.ipynb

* Update gridlabd-geodata

* Update geodata_example.sh

* Update GridLAB-D GeoData Tutorial.ipynb

* Update version check so it works from the command line too

* Add more support for JSON and Excel file I/O

* Add relative distance option to geodata distance package

* Update gridlabd-geodata

* Update GridLAB-D GeoData Tutorial.ipynb

* Update Geodata.md

* Update Geodata.md

* Update test_geodata_elevation.glm

* Update requirements.txt

* Update requirements.txt

* Update Geodata.md

* Update gridlabd-geodata

* Update gridlabd-geodata

* Update GridLAB-D GeoData Tutorial.ipynb

* Fix input delimiter problem with address (again)

* Fix autotest errors

* Fix geodata autotest problems on github actions

* Create Geodata.md

* Update Geodata.md

* Add tutorial docs generation

* Tutorial update

* Update test_geodata.glm

* Making pole test optional until it's fixed

* Update the geodata tutorial

* Update geodata tutorial

* Update docs and add --resolution_id option

* Update Geodata.md

* Fix waypoints id calcs

* Update geodata_elevation.py

* Create geodata_vegetation.py

* Add geodata_vegetation.py

* Update geodata_vegetation.py

* Add vegetation geodata package

* Update geodata_vegetation.py

* Update geodata_vegetation.py

* Update vegetation geodata package

* Update geodata_vegetation.py

* Update geodata_vegetation.py

* Add link to CFO repo for user registration

* Update vegetation package and docs

* Update geodata tutorial

* Tutorial update

* Update tutorial docs

* Update distance package to support relative distances

* Update with a second example

* Fix distance calcs when resolution is too large

* Merge branch 'develop' into develop-add-vegetation-geodata

* Fix geodata autotest

* fix distance calulation issue

* fix the bug for calculation linesag at arbitrary points

* Update geodata_vegetation.py

* Create update_vegetation.py

* Update update_vegetation.py

* Update update_vegetation.py

* Fix image tiling

* Code for downloading 0.1 degree vegetation tiles and uploading them to AWS

* Update geodata_vegetation.py

* Update Geodata.ipynb

* Update Geodata.md

* Update Geodata.ipynb

* Update Geodata.md

* Update geodata_powerline.py

* Update Geodata.ipynb

* Update Geodata.md

* Update docs

* change pole height

* change pole geight

* Added line contact probability test

* Update docs

* Update Geodata.ipynb

* Update CFO authenticator

* Added interpolation of elevation

* Update geodata_elevation.py

* update elevation intepolation for geodata_elevation

* Update the initial tension calculation

* Fixed python notebook for geodata

* Update Geodata.md

* Update develop.yml

* merge updates

* Update path_example.csv

* Update gridlabd-library

* Update geodata_vegetation.py

* Update geodata_vegetation.py

* Create geodata.conf

* Update update_vegetation.py

* Update Geodata.ipynb

* Update CFO authentication

* Update Geodata.ipynb

* Update geodata_vegetation.py

* Fix autotest validation

* Update test_geodata.txt

* Update geodata_vegetation.py

* Update geodata_vegetation.py

* Update geodata_powerline.py

* Update geodata_powerline.py

* Update geodata_powerline.py

Co-authored-by: Fuhong Xie <fxie2@stanford.edu>
Co-authored-by: Fuhong Xie <fxie2@slac.stanford.edu>
Co-authored-by: Alyona Teyber <Ivanova.alyona5@gmail.com>

gldcore/geodata/geodata_powerline.py

@@ -163,21 +163,12 @@ def get_air_properties(temp):
     """
     # calculate the constant parameters
     # reference: http://home.eng.iastate.edu/~jdm/wind/TransmissionLineLoadingDesignCriteriaAndHTS.pdf
-    temp_x = np.arange(0.0, 110.0, 10) # unit: DegC
-    specific_mass = [1.29, 1.25, 1.20, 1.17, 1.13, 1.09, 1.06, 1.03, 1.0, 0.97, 0.95] # unit: kg/m3
-    thermal_conductivity = [0.0243, 0.025, 0.0257, 0.0265, 0.0272, 0.028, 0.0287,0.0294, 0.0301, 0.0309, 0.0316] # unit: W/(m*DegC)
-    dynamic_viscosity = [0.175e-4, 0.18e-4, 0.184e-4, 0.189e-4, 0.194e-4, 0.199e-4, 0.203e-4, 0.208e-4, 0.213e-4, 0.217e-4, 0.222e-4] # unit: N*s/m2
 
-    f_specific_mass = interp1d(temp_x, specific_mass, fill_value = "extrapolate")
-    specific_mass_temp = f_specific_mass(temp)
+    specific_mass = -0.00342727*temp + 1.275 
+    thermal_conductivity = 7.32727273e-05*temp + 0.02428182
+    dynamic_viscosity = 4.7e-08*temp + 1.75045455e-05
 
-    f_thermal_conductivity = interp1d(temp_x, thermal_conductivity, fill_value = "extrapolate")
-    thermal_conductivity_temp = f_thermal_conductivity(temp)
-
-    f_dynamic_viscosity = interp1d(temp_x, dynamic_viscosity, fill_value = "extrapolate")
-    dynamic_viscosity_temp = f_dynamic_viscosity(temp)
-
-    return specific_mass_temp, thermal_conductivity_temp, dynamic_viscosity_temp
+    return specific_mass, thermal_conductivity, dynamic_viscosity
 
 hold_values = {}
 def hold0(name,value=None):
@@ -221,8 +212,8 @@ def get_line_tension_coefficient(d_hori):
     """
     # reference: IEC 60826:2017
     # fit of data is used
-    # 0.2106718346 + 0.0003126614987*d_hori
-    return 0.50723325 - 0.29656325/(0.05*d_hori+1)
+    # 0.50723325 - 0.29656325/(0.05*d_hori+1)
+    return 0.2106718346 + 0.0003126614987*d_hori
 
 def get_distance(pos1, pos2):
     """Compute haversine distance between two locations
@@ -366,28 +357,29 @@ def linesag(data):
                 p0 = [line['latitude'],line['longitude']]
                 d0 = 0.0
                 z0 = elevation + line['pole_height']
-                sag = z0 - elevation
+                sag = line['pole_height']
                 result[id] = sag
 
             else: # end of a line segment
 
                 p1 = [line['latitude'],line['longitude']]
                 d1 = get_distance(p0,p1)
                 z1 = elevation + line['pole_height']
-                sag = z1 - elevation
+                sag = line['pole_height']
                 result[id] = sag
 
                 # compute linesag at waypoints
                 for n,l in ld.items():
                     p = [l['latitude'],l['longitude']]
-                    d_hori = get_distance(p0,p)
+                    d_hori = get_distance(p0,p1)
                     try:
                         elevation = l['elevation']
                     except:
                         pass
                     sag = get_sag_value(d_hori,line,cable,p0,p,z0,z1,
                         power_flow,global_horizontal_irradiance,ground_reflectance,
                         ice_thickness,wind_direction,air_temperature,wind_speed,ice_density)
+                    # the line sag is defined as the distance between the interested point on the line and the ground
                     sag = sag - elevation
                     result[n] = round(sag,OPTIONS["precision"]["linesag"])
 
@@ -423,52 +415,54 @@ def get_sag_value(d_hori,line,cable,p0,p1,z0,z1,
     span = sqrt(d_hori*d_hori + d_vert*d_vert)
     rts = cable['rated_tensile_strength']
     unit_weight = cable['unit_weight']
-    # if d_vert > 30:
-    #     k_init = min(k_init, 3.0*unit_weight*span*span/(2*d_vert*rts))
+    diameter = cable['diameter']
+    air_mass, k_f, air_viscosity = get_air_properties(air_temperature)
+    # calculate the line heading and the angle between line heading and wind direction
+    try:
+        line_angle = line['heading']
+    except:
+        line_angle = 180*atan2(p1[0]-p0[0],p1[1]-p0[1])/np.pi
+    phi = (wind_direction - line_angle)*np.pi/180
+    # calculate the new line sag at loaded condition
+    ice_unit_weight = ice_density*np.pi*ice_thickness*(diameter+ice_thickness)*g
+    wind_unit_weight = 0.5*air_mass*(wind_speed*sin(phi))**2 *(diameter+2*ice_thickness)
+    total_unit_weight = sqrt(wind_unit_weight**2+(unit_weight+ice_unit_weight)**2)
+    if d_vert/d_hori > 0.1:
+        k_init = min(k_init, total_unit_weight*d_hori*d_hori/(2*d_vert*rts))
     H_init = rts*k_init
     sag_init = unit_weight*span*span/(8*H_init)
+    # for Q_I
     P_rated = power_flow
     Vll_rated = cable['voltage_rating']
     Irms =  P_rated/(sqrt(3)*Vll_rated)
     R_20C = cable['nominal_resistance']
     coeff_Al = cable['resistivity']
     Q_I_coeff_first = Irms*Irms*R_20C*coeff_Al
     Q_I_coeff_constant = Irms*Irms*R_20C*(1-coeff_Al*(20.0+273.0))
+    # for Q_S
     k_a = 1.0 - cable['reflectivity'] # solar radiation absorption coefficient
     GHI = global_horizontal_irradiance # unit: W/m2
     k_g = ground_reflectance # ground reflect
-    diameter = cable['diameter']
     Q_S_constant = k_a*(diameter+2*ice_thickness)*(1+k_g)*GHI
-
-    try:
-        line_angle = line['heading']
-    except:
-        line_angle = 180*atan2(p1[0]-p0[0],p1[1]-p0[1])/np.pi
-
-    phi = (wind_direction - line_angle)*np.pi/180
+    # for Q_C
     k_angle = 1.194 - cos(phi) + 0.194*cos(2*phi) + 0.368*sin(2*phi)
-    air_mass, k_f, air_viscosity = get_air_properties(air_temperature)
     Nre = wind_speed * (diameter + 2*ice_thickness) * air_mass / air_viscosity
-    if wind_speed > 1.1:
+    if wind_speed < 0.82:
         Q_C_coeff_first = k_angle*(1.01+1.35*Nre**0.52)*k_f
         Q_C_constant = -k_angle*(1.01+1.35*Nre**0.52)*k_f*(air_temperature+273.0)
     else:
         Q_C_coeff_first = k_angle*0.754*(Nre**0.6)*k_f
         Q_C_constant = -k_angle*0.754*(Nre**0.6)*k_f*(air_temperature+273.0)
     # for Q_R
     k_e = cable['emissivity']
-    Q_R_constant = -5.6704e-8*k_e*(diameter+2.0*ice_thickness)*(air_temperature+273.0)**4
-    Q_R_coeff_fourth = 5.6704e-8*k_e*(diameter+2.0*ice_thickness)
+    Q_R_constant = -5.6704e-8*k_e*(diameter+2.0*ice_thickness)*np.pi*(air_temperature+273.0)**4
+    Q_R_coeff_fourth = 5.6704e-8*k_e*(diameter+2.0*ice_thickness)*np.pi
     # for new conductor temp under loading
     coef_sag = [-Q_R_coeff_fourth,0.0,0.0,Q_I_coeff_first-Q_C_coeff_first,Q_I_coeff_constant+Q_S_constant-Q_C_constant-Q_R_constant]
     r = np.roots(coef_sag)
     r = r[~np.iscomplex(r)]
     temp_load = np.absolute(r[r > 0.0]) # unit: K
     temp_load = temp_load - 273.0 # unit: DegC
-    # calculate the new line sag at loaded condition
-    ice_unit_weight = ice_density*np.pi*ice_thickness*(diameter+ice_thickness)*g
-    wind_unit_weight = 0.5*air_mass*(wind_speed*sin(phi))**2 *(diameter+2*ice_thickness)
-    total_unit_weight = sqrt(wind_unit_weight**2+(unit_weight+ice_unit_weight)**2)
 
     area = cable['conductor_crosssection_area']
     elasticity = cable['elasticity']
@@ -480,18 +474,20 @@ def get_sag_value(d_hori,line,cable,p0,p1,z0,z1,
     r = np.roots(coef_H)
     r = r[~np.iscomplex(r)]
     H_load = np.absolute(r[r > 0.0])
-    #
+
     sag_load = total_unit_weight*span*span/(8*H_load)
     sag_angle = atan(wind_unit_weight/(ice_unit_weight+unit_weight))
     C_catenary = H_load / (ice_unit_weight+unit_weight)
     if z0 > z1:
-        d0_hori = d_hori*(1+d_vert/(4*sag_load))/2
+        d0_hori = d_hori/2 + H_load*d_vert/(total_unit_weight*d_hori)
+        d1_hori = d_hori - d0_hori
         sag0 = total_unit_weight*d0_hori**2 /(2*H_load)
         dt = get_distance(p0,p1)
         sag0_cosh = sag0 - C_catenary*(np.cosh((dt-d0_hori)/C_catenary)-1)
         sag_elevation = z0 - sag0_cosh*cos(sag_angle)
     else:
-        d0_hori = d_hori*(1-d_vert/(4*sag_load))/2
+        d0_hori = d_hori/2 - H_load*d_vert/(total_unit_weight*d_hori)
+        d1_hori = d_hori - d0_hori
         sag0 = total_unit_weight*d0_hori**2 /(2*H_load)
         dt = get_distance(p0,p1)
         sag0_cosh = sag0 - C_catenary*(np.cosh((dt-d0_hori)/C_catenary)-1)