diff --git a/Jenkinsfile b/Jenkinsfile new file mode 100644 index 0000000..ed58203 --- /dev/null +++ b/Jenkinsfile @@ -0,0 +1,90 @@ +pipeline { + agent any + stages { + stage('develop build and upload') { + when { + branch "develop" + } + steps { + script { + try { + sh ''' + docker login http://172.16.102.3:30648 -u ${reg_username} -p ${reg_passwd} + docker build -t 172.16.102.3:30648/maimai/${service_name}:latest . + docker push 172.16.102.3:30648/maimai/${service_name}:latest + kubectl rollout restart deployment ${service_name} -n dev + ''' + } catch(err) { + sh 'exit 1' + } + } + } + } + stage('test build only') { + when { + expression {ref ==~ 'refs/tags/test-.*' } + } + steps { + script { + try { + sh ''' + docker login http://172.16.102.3:30648 -u ${reg_username} -p ${reg_passwd} + docker build -t 172.16.102.3:30648/maimai/${service_name}:${result} . + docker push 172.16.102.3:30648/maimai/${service_name}:${result} + ''' + } catch(err) { + exit 1 + } + } + } + } + stage('master build only') { + when { + expression {ref ==~ 'refs/tags/prod-.*' } + } + steps { + script { + try { + sh ''' + docker login http://172.16.102.3:30648 -u ${reg_username} -p ${reg_passwd} + docker build -t 172.16.102.3:30648/maimai/${service_name}:${result} . + docker push 172.16.102.3:30648/maimai/${service_name}:${result} + ''' + } catch(err) { + exit 1 + } + } + } + } + } + environment { + reg_username = 'maimai' + reg_passwd = 'M9hUQk4Ti0l0lHZi' + service_name = 'xj-irrigation-model' + result = sh(script: """echo $ref | awk -F"/" '{print \$NF}'""", returnStdout: true).trim() + git_commit_msg = sh (script: 'git log -1 --pretty=%B ${GIT_COMMIT}', returnStdout: true).trim() + git_commit_user = sh (script: 'git show -s --pretty=%an', returnStdout: true).trim() + } + post { + success { + sh """curl --location --request POST 'https://qyapi.weixin.qq.com/cgi-bin/webhook/send?key=d31913db-bb53-4081-9e2d-f707111dbee1' \ + --data-raw '{\"msgtype\":\"text\",\"text\":{\"content\":\"成功: [${env.JOB_NAME} [${env.BUILD_NUMBER}]](${git_commit_msg}) @${git_commit_user}\"}}'""" + } + failure { + sh """curl --location --request POST 'https://qyapi.weixin.qq.com/cgi-bin/webhook/send?key=d31913db-bb53-4081-9e2d-f707111dbee1' \ + --data-raw '{\"msgtype\":\"text\",\"text\":{\"content\":\"失败: [${env.JOB_NAME} [${env.BUILD_NUMBER}]](${git_commit_msg}) @${git_commit_user}\"}}'""" + } + } + triggers { + GenericTrigger( + genericVariables: [[key: 'ref', value: '$.ref']], + causeString: 'Triggered on $ref', + token: 'xj-irrigation-model', + printContributedVariables: true, + printPostContent: true, + silentResponse: false, + regexpFilterText: '$ref/'+BRANCH_NAME, + regexpFilterExpression: '(^refs/heads/develop/develop$)|(^refs/tags/test-.*/test$)|(^refs/tags/prod-.*/master$)' + ) + } +} \ No newline at end of file diff --git a/README.en.md b/README.en.md new file mode 100644 index 0000000..80ab61d --- /dev/null +++ b/README.en.md @@ -0,0 +1,36 @@ +# xj-irrigation-model + +#### Description +智能水肥模型 + +#### Software Architecture +Software architecture description + +#### Installation + +1. xxxx +2. xxxx +3. xxxx + +#### Instructions + +1. xxxx +2. xxxx +3. xxxx + +#### Contribution + +1. Fork the repository +2. Create Feat_xxx branch +3. Commit your code +4. Create Pull Request + + +#### Gitee Feature + +1. You can use Readme\_XXX.md to support different languages, such as Readme\_en.md, Readme\_zh.md +2. Gitee blog [blog.gitee.com](https://blog.gitee.com) +3. Explore open source project [https://gitee.com/explore](https://gitee.com/explore) +4. The most valuable open source project [GVP](https://gitee.com/gvp) +5. The manual of Gitee [https://gitee.com/help](https://gitee.com/help) +6. The most popular members [https://gitee.com/gitee-stars/](https://gitee.com/gitee-stars/) diff --git a/README.md b/README.md index 5e0d62d..302c50a 100644 --- a/README.md +++ b/README.md @@ -1,3 +1,37 @@ -# irrigation-model +# xj-irrigation-model -判断是否需要灌溉 \ No newline at end of file +#### 介绍 +智能水肥模型 + +#### 软件架构 +软件架构说明 + + +#### 安装教程 + +1. xxxx +2. xxxx +3. xxxx + +#### 使用说明 + +1. xxxx +2. xxxx +3. xxxx + +#### 参与贡献 + +1. Fork 本仓库 +2. 新建 Feat_xxx 分支 +3. 提交代码 +4. 新建 Pull Request + + +#### 特技 + +1. 使用 Readme\_XXX.md 来支持不同的语言，例如 Readme\_en.md, Readme\_zh.md +2. Gitee 官方博客 [blog.gitee.com](https://blog.gitee.com) +3. 你可以 [https://gitee.com/explore](https://gitee.com/explore) 这个地址来了解 Gitee 上的优秀开源项目 +4. [GVP](https://gitee.com/gvp) 全称是 Gitee 最有价值开源项目，是综合评定出的优秀开源项目 +5. Gitee 官方提供的使用手册 [https://gitee.com/help](https://gitee.com/help) +6. Gitee 封面人物是一档用来展示 Gitee 会员风采的栏目 [https://gitee.com/gitee-stars/](https://gitee.com/gitee-stars/) diff --git a/config.nifi b/config.nifi new file mode 100644 index 0000000..fd218ec --- /dev/null +++ b/config.nifi @@ -0,0 +1,11 @@ +[postgresql] +host = localhost +port = 5432 +database = datastore +user = postgres +password = postgres +schema = public +tablename=irrigation_data + +[logging] +level = INFO \ No newline at end of file diff --git a/dockerfile b/dockerfile new file mode 100644 index 0000000..7f9796c --- /dev/null +++ b/dockerfile @@ -0,0 +1,9 @@ +FROM python:3.10.5 + +WORKDIR /app +RUN pip config set global.index-url https://mirrors.aliyun.com/pypi/simple + +COPY . . +RUN pip install -r requirements.txt +ENTRYPOINT [ "python", "main.py" ] +EXPOSE 8001 \ No newline at end of file diff --git a/irrgiation/__init__.py b/irrgiation/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/irrgiation/__pycache__/__init__.cpython-310.pyc b/irrgiation/__pycache__/__init__.cpython-310.pyc new file mode 100644 index 0000000..957bca8 Binary files /dev/null and b/irrgiation/__pycache__/__init__.cpython-310.pyc differ diff --git a/irrgiation/__pycache__/crop_growth_stages.cpython-310.pyc b/irrgiation/__pycache__/crop_growth_stages.cpython-310.pyc new file mode 100644 index 0000000..c28890a Binary files /dev/null and b/irrgiation/__pycache__/crop_growth_stages.cpython-310.pyc differ diff --git a/irrgiation/__pycache__/dailyEvaporation.cpython-310.pyc b/irrgiation/__pycache__/dailyEvaporation.cpython-310.pyc new file mode 100644 index 0000000..74b779e Binary files /dev/null and b/irrgiation/__pycache__/dailyEvaporation.cpython-310.pyc differ diff --git a/irrgiation/__pycache__/db_connect.cpython-310.pyc b/irrgiation/__pycache__/db_connect.cpython-310.pyc new file mode 100644 index 0000000..434ea44 Binary files /dev/null and b/irrgiation/__pycache__/db_connect.cpython-310.pyc differ diff --git a/irrgiation/__pycache__/irrigateDecisionDemo.cpython-310.pyc b/irrgiation/__pycache__/irrigateDecisionDemo.cpython-310.pyc new file mode 100644 index 0000000..d4c56e8 Binary files /dev/null and b/irrgiation/__pycache__/irrigateDecisionDemo.cpython-310.pyc differ diff --git a/irrgiation/__pycache__/mathuntils.cpython-310.pyc b/irrgiation/__pycache__/mathuntils.cpython-310.pyc new file mode 100644 index 0000000..84f6833 Binary files /dev/null and b/irrgiation/__pycache__/mathuntils.cpython-310.pyc differ diff --git a/irrgiation/__pycache__/soil_Ke.cpython-310.pyc b/irrgiation/__pycache__/soil_Ke.cpython-310.pyc new file mode 100644 index 0000000..597fc32 Binary files /dev/null and b/irrgiation/__pycache__/soil_Ke.cpython-310.pyc differ diff --git a/irrgiation/__pycache__/weatherAndSoilDataRequest.cpython-310.pyc b/irrgiation/__pycache__/weatherAndSoilDataRequest.cpython-310.pyc new file mode 100644 index 0000000..bf3190b Binary files /dev/null and b/irrgiation/__pycache__/weatherAndSoilDataRequest.cpython-310.pyc differ diff --git a/irrgiation/crop_growth_stages.py b/irrgiation/crop_growth_stages.py new file mode 100644 index 0000000..2cd8f10 --- /dev/null +++ b/irrgiation/crop_growth_stages.py @@ -0,0 +1,68 @@ +crop_growth_stages={ + '甘草': { + 'initial': { + 'root_length': {5: 0.8, 10: 0.2}, # 主根长：80%概率在0-15cm，平均约8cm + 'duration_days': 20, # 6月播种后1-2周 + 'kc': 0.40, # 作物系数 + 'time_range': '2025-07-08至2025-07-28', + 'hight':0.03,#m + 'k_rain':0.25 + }, + 'development': { + 'root_length': {20: 0.6, 45: 0.4}, # 主根长：60%概率在20cm，40%概率在45cm，平均约30cm + 'duration_days': 42, # 发芽后-8月20日左右（约2个月） + 'kc': 0.9, # 作物系数 + 'time_range': '2025-07-29至2025-09-10', + 'hight':0.15, + 'k_rain':0.25 + }, + 'mid': { + 'root_length': {30: 0.5, 80: 0.5}, # 主根长：50%概率在30cm，50%概率在80cm，平均约55cm + 'duration_days': 30, # 8月20日-9月下旬 + 'kc': 1.1, # 作物系数 + 'time_range': '2025-09-11至2025-10-10', + 'hight':0.45, + 'soil_evaporation_depth': 0.65, + 'k_rain':0.4 + }, + 'growth_slowly': { + 'root_length':{40: 0.4,60:0.4, 80: 0.2}, # 根据数据估算平均主根长 + 'duration_days': 31, # 9月下旬-10月下旬 + 'kc': 0.6, + 'time_range': '2025-10-11至2025-11-10', + 'hight':0.6 + }, + 'dormant_period': { + 'root_length': {40: 0.3,60:0.4, 80: 0.2}, # 冬季休眠，主根基本停止生长 + 'duration_days': 140, # 11月-次年3月 + 'kc': 0.4, + 'time_range': '2025-11-11至2026-03-31', + 'hight':0.8, + 'k_rain':0.5 + }, + 'fanqing_year2': { + 'root_length': {40: 0.2,60:0.6, 80: 0.2}, # 返青期主根长变化不大 + 'duration_days': 30, # 3月下旬-4月下旬 + 'kc': 0.35, + 'time_range': '2026-04-01至2026-04-30', + 'hight':0.5, + 'k_rain':0.5 + }, + 'development_year2': { + 'root_length': {40: 0.1,60:0.7, 80: 0.2}, # 根据数据估算平均主根长 + 'duration_days': 90, # 4月下旬-7月下旬 + 'kc': 0.9, + 'time_range': '2026-05-01至2026-07-30', + 'hight':0.45, + 'k_rain':0.65 + }, + 'maturity_year2': { + 'root_length':{40: 0.1,60:0.6, 80: 0.3}, # 根据数据估算平均主根长 + 'duration_days': 40, # 7月下旬-8月下旬 + 'kc': 0.6, + 'time_range': '2026-08-01至2026-09-10', + 'hight':0.65, + 'k_rain':0.5 + } + } +} \ No newline at end of file diff --git a/irrgiation/dailyEvaporation.py b/irrgiation/dailyEvaporation.py new file mode 100644 index 0000000..f68c8ed --- /dev/null +++ b/irrgiation/dailyEvaporation.py @@ -0,0 +1,365 @@ +from datetime import datetime, timedelta +from pyeto import fao,convert +from pyeto.fao import svp_from_t, avp_from_rhmean, delta_svp, psy_const,atm_pressure +import pandas as pd +#1、计算日参考蒸散量 +def calculateET0(tmean, tmax, tmin, rh_mean, wind_speed, solar_rad, elevation, lat, doy): + # 计算饱和水汽压和实际水汽压 + svp_max = svp_from_t(tmax) + svp_min = svp_from_t(tmin) + svp_mean =(svp_max+svp_min)/2 + avp = avp_from_rhmean(svp_max, svp_min, rh_mean) + + # 计算净辐射(简化版) + # 将纬度从度转换为弧度 + lat_rad = convert.deg2rad(lat) + sol_dec = fao.sol_dec(doy) + sha = fao.sunset_hour_angle(lat_rad, sol_dec) + ird = fao.inv_rel_dist_earth_sun(doy) + et_rad = fao.et_rad(lat_rad, sol_dec, sha, ird) + cs_rad = fao.cs_rad(elevation, et_rad) + nisr = fao.net_in_sol_rad(solar_rad) + nolr = fao.net_out_lw_rad(tmin, tmax, solar_rad, cs_rad, avp) + rn = fao.net_rad(nisr, nolr) + deltasvp= delta_svp(tmean) + psy = psy_const(atm_pressure(elevation)) + # 计算ET₀ + et0 = fao.fao56_penman_monteith( + net_rad=rn, + t=tmean, + ws=wind_speed, + svp=svp_mean, + avp=avp, + delta_svp=deltasvp, + psy=psy + ) + return et0 +def cal_et0_List(planting_date,monitor_date,date_params): + et0_list = [] + current_date = planting_date + while current_date < monitor_date: + result_params=date_params[current_date.strftime("%Y-%m-%d")] + et0=calculateET0(result_params['tmean'], result_params['tmax'], result_params['tmin'], result_params['rh_mean'], result_params['wind_speed'], + result_params['solar_rad'], result_params['elevation'],result_params['lat'], result_params['doy']) + et0_list.append(et0) + current_date += timedelta(days=1) # 增加一天 + return et0_list + + +#计算实际每天作物蒸散量 +# 1. 定义甘草生育阶段与作物系数（参考苜蓿） +# growth_stages = { +# "initial": {"duration_days": 20, "kc": 0.40}, # 生长初期 +# "development": {"duration_days": 40, "kc": 0.90}, #生长快速期 +# "mid": {"duration_days": 20, "kc": 1.1}, # 生长缓慢期 +# "growth_slowly": {"duration_days": 30, "kc": 0.6}, # 生长缓慢期 +# "dormant_period": {"duration_days": 150, "kc": 0.40}, # 营养储备期 +# "fanqing_year2":{"duration_days": 30, "kc": 0.35},#返青期 +# "development_year2":{"duration_days": 90, "kc": 0.9}, #第二年快速生长期 +# "maturity _year2": {"duration_days": 40, "kc": 0.6}, # 第二年成熟期 +# } +growth_stages_roots1 = { + 'initial': { + 'root_length': {5: 0.8, 10: 0.2}, # 主根长：80%概率在0-15cm，平均约8cm + 'duration_days': 23, # 6月播种后1-2周 + 'kc': 0.40, # 作物系数 + 'time_range': '2025-06-17至2025-07-10', + 'hight':0.03,#m + 'k_rain':0.25 + }, + 'development': { + 'root_length': {20: 0.6, 45: 0.4}, # 主根长：60%概率在20cm，40%概率在45cm，平均约30cm + 'duration_days': 40, # 发芽后-8月20日左右（约2个月） + 'kc': 0.9, # 作物系数 + 'time_range': '2025-07-10至2025-08-20', + 'hight':0.15, + 'k_rain':0.25 + }, + 'mid': { + 'root_length': {30: 0.5, 80: 0.5}, # 主根长：50%概率在30cm，50%概率在80cm，平均约55cm + 'duration_days': 30, # 8月20日-9月下旬 + 'kc': 1.1, # 作物系数 + 'time_range': '2025-08-21至2025-09-20', + 'hight':0.45, + 'soil_evaporation_depth': 0.65, + 'k_rain':0.4 + }, + 'growth_slowly': { + 'root_length':{40: 0.4,60:0.4, 80: 0.2}, # 根据数据估算平均主根长 + 'duration_days': 30, # 9月下旬-10月下旬 + 'kc': 0.6, + 'time_range': '2025-09-21至2025-10-20', + 'hight':0.6 + }, + 'dormant_period': { + 'root_length': {40: 0.3,60:0.4, 80: 0.2}, # 冬季休眠，主根基本停止生长 + 'duration_days': 160, # 11月-次年3月 + 'kc': 0.4, + 'time_range': '2025-10-21至2026-03-30', + 'hight':0.8, + 'k_rain':0.5 + }, + 'fanqing_year2': { + 'root_length': {40: 0.2,60:0.6, 80: 0.2}, # 返青期主根长变化不大 + 'duration_days': 30, # 3月下旬-4月下旬 + 'kc': 0.35, + 'time_range': '2026-04-01至2026-04-30', + 'hight':0.5, + 'k_rain':0.5 + }, + 'development_year2': { + 'root_length': {40: 0.1,60:0.7, 80: 0.2}, # 根据数据估算平均主根长 + 'duration_days': 90, # 4月下旬-7月下旬 + 'kc': 0.9, + 'time_range': '2026-05-01至2026-07-30', + 'hight':0.45, + 'k_rain':0.65 + }, + 'maturity_year2': { + 'root_length':{40: 0.1,60:0.6, 80: 0.3}, # 根据数据估算平均主根长 + 'duration_days': 40, # 7月下旬-8月下旬 + 'kc': 0.6, + 'time_range': '2026-08-01至2026-09-10', + 'hight':0.65, + 'k_rain':0.5 + } +} +# growth_stages_roots = { +# 'initial': { +# 'root_length': {5: 0.8, 10: 0.2}, # 主根长：80%概率在0-15cm，平均约8cm +# 'duration_days': 20, # 6月播种后1-2周 +# 'kc': 0.40, # 作物系数 +# 'time_range': '2025-07-08至2025-07-28', +# 'hight':0.03,#m +# 'k_rain':0.25 +# }, +# 'development': { +# 'root_length': {20: 0.6, 45: 0.4}, # 主根长：60%概率在20cm，40%概率在45cm，平均约30cm +# 'duration_days': 42, # 发芽后-8月20日左右（约2个月） +# 'kc': 0.9, # 作物系数 +# 'time_range': '2025-07-29至2025-09-10', +# 'hight':0.15, +# 'k_rain':0.25 +# }, +# 'mid': { +# 'root_length': {30: 0.5, 80: 0.5}, # 主根长：50%概率在30cm，50%概率在80cm，平均约55cm +# 'duration_days': 30, # 8月20日-9月下旬 +# 'kc': 1.1, # 作物系数 +# 'time_range': '2025-09-11至2025-10-10', +# 'hight':0.45, +# 'soil_evaporation_depth': 0.65, +# 'k_rain':0.4 +# }, +# 'growth_slowly': { +# 'root_length':{40: 0.4,60:0.4, 80: 0.2}, # 根据数据估算平均主根长 +# 'duration_days': 31, # 9月下旬-10月下旬 +# 'kc': 0.6, +# 'time_range': '2025-10-11至2025-11-10', +# 'hight':0.6 +# }, +# 'dormant_period': { +# 'root_length': {40: 0.3,60:0.4, 80: 0.2}, # 冬季休眠，主根基本停止生长 +# 'duration_days': 140, # 11月-次年3月 +# 'kc': 0.4, +# 'time_range': '2025-11-11至2026-03-31', +# 'hight':0.8, +# 'k_rain':0.5 +# }, +# 'fanqing_year2': { +# 'root_length': {40: 0.2,60:0.6, 80: 0.2}, # 返青期主根长变化不大 +# 'duration_days': 30, # 3月下旬-4月下旬 +# 'kc': 0.35, +# 'time_range': '2026-04-01至2026-04-30', +# 'hight':0.5, +# 'k_rain':0.5 +# }, +# 'development_year2': { +# 'root_length': {40: 0.1,60:0.7, 80: 0.2}, # 根据数据估算平均主根长 +# 'duration_days': 90, # 4月下旬-7月下旬 +# 'kc': 0.9, +# 'time_range': '2026-05-01至2026-07-30', +# 'hight':0.45, +# 'k_rain':0.65 +# }, +# 'maturity_year2': { +# 'root_length':{40: 0.1,60:0.6, 80: 0.3}, # 根据数据估算平均主根长 +# 'duration_days': 40, # 7月下旬-8月下旬 +# 'kc': 0.6, +# 'time_range': '2026-08-01至2026-09-10', +# 'hight':0.65, +# 'k_rain':0.5 +# } +# } +# --------------------------- +# 2. 计算每日Kc值（根据生育阶段动态分配） +# --------------------------- +def get_daily_kc(growth_stages, total_days=423): + daily_kc = [] + for stage in growth_stages.values(): + daily_kc.extend([stage["kc"]] * stage["duration_days"]) + # 如果总天数超出定义，用末期Kc填充剩余天数 + if len(daily_kc) < total_days: + daily_kc.extend([growth_stages["maturity_year2"]["kc"]] * (total_days - len(daily_kc))) + return daily_kc[:total_days] # 截断至总天数 + +# --------------------------- +# 3. 计算每日实际蒸散量ETc(得到历史列表) +# --------------------------- +def calculate_etc_List(et0_daily, planting_date,monitor_date, kc_daily): + """ + 计算每日ETc（实际蒸散量） + :param eto_daily: 每日参考蒸散量列表（单位：cm/d） + :param planting_date: 播种日期，格式"YYYY-MM-DD" + :param kc_daily: 每日Kc值列表 + :return: 日期列表、ETc列表 + """ + start_date = planting_date + delta=monitor_date-start_date + days = delta.days + dates = [start_date + timedelta(days=i) for i in range(days)] + date_strings = [date.strftime("%Y-%m-%d") for date in dates] + kc_daily = kc_daily[:days] + etc = [eto * kc for eto, kc in zip(et0_daily, kc_daily)] + return etc,date_strings,kc_daily +#基于作物蒸散和土壤蒸发计算实际蒸散量 +def calculate_etc_List_new(et0_daily, planting_date,monitor_date, kc_daily,ke_daily): + """ + 计算每日ETc（实际蒸散量） + :param eto_daily: 每日参考蒸散量列表（单位：cm/d） + :param planting_date: 播种日期，格式"YYYY-MM-DD" + :param kc_daily: 每日Kc值列表 + :return: 日期列表、ETc列表 + """ + start_date = planting_date + delta=monitor_date-start_date + days = delta.days + dates = [start_date + timedelta(days=i) for i in range(days)] + date_strings = [date.strftime("%Y-%m-%d") for date in dates] + kc_daily = kc_daily[:days] + etc = [eto * (kc+ke) for eto, kc,ke in zip(et0_daily, kc_daily,ke_daily)] + return etc,date_strings,kc_daily +def calculate_etc_dayily(et0_daily, planting_date,monitor_date, kc_daily): + """ + 计算每日ETc（实际蒸散量） + :param eto_daily: 每日参考蒸散量列表（单位：cm/d） + :param planting_date: 播种日期，格式"YYYY-MM-DD" + :param kc_daily: 每日Kc值列表 + :return: 日期列表、ETc列表 + """ + start_date = planting_date + delta=monitor_date-start_date + days = delta.days + # dates = [start_date + timedelta(days=i) for i in range(len(kc_daily))] + kc_daily = kc_daily[days] + etc = et0_daily*kc_daily + return etc +def export_etc_et0_date(et0_list,etc_list,dates,kc_list,ke_daily,soil_data,rain_daily,irrigate_daily,dk_name): + df = pd.DataFrame({ + "Date": dates, + "ET0": [float(f"{x:.3f}") for x in et0_list], + "KC": kc_list, + "KE":ke_daily, + "ETc": [float(f"{x:.3f}") for x in etc_list], + "rain_effective":rain_daily, + "I_effective":irrigate_daily, + "VWC_fc":0.2884, + "VWC_wp":0.1086, + "dkbm":dk_name + }) + df['VWC_SUR'] = 0.0 + df['VWC_MID'] = 0.0 + df['VWC_BOT'] = 0.0 + for i in range(len(df)): + date=df.at[i, 'Date'] + df.at[i, 'VWC_SUR']=round(float(soil_data["SOIL_MOISTURE_SURFACE"][date]),3) + df.at[i, 'VWC_MID'] = round(float(soil_data["SOIL_MOISTURE_MIDDLE"][date]),3) + df.at[i, 'VWC_BOT'] = round(float(soil_data["SOIL_MOISTURE_BOTTOM"][date]),3) + return df + + +def estimate_mad_from_et(max_et): + + # 根据最大日蒸散量查表确定MAD值 + if max_et < 2: + mad = 0.68 + elif 2 <= max_et < 3: + mad = 0.58 + elif 3 <= max_et < 4: + mad = 0.48 + elif 4 <= max_et < 5: + mad = 0.4 + elif 5 <= max_et < 6: + mad = 0.35 + elif 6 <= max_et < 7: + mad = 0.33 + elif 7 <= max_et < 8: + mad = 0.28 + elif 8 <= max_et < 9: + mad = 0.25 + else: # max_et >= 8 + mad = 0.23 + return mad +def calculate_weighted_soil_parameters(root_distribution, soil_layers): + # 初始化累计变量 + sum_fc = 0.0 + sum_wp = 0.0 + # 处理每个根系分布层 + for root_depth, root_ratio in root_distribution['root_length'].items(): + # 找到包含该深度的土壤层 + for layer in soil_layers: + min_depth, max_depth = layer['depth_range'] + if min_depth <= root_depth < max_depth: + # 计算该层贡献 (简单比例分配) + layer_thickness = max_depth - min_depth + depth_in_layer = min(root_depth, max_depth) - min_depth + contribution_ratio = depth_in_layer / layer_thickness * root_ratio + sum_fc += layer['theta_33'] * contribution_ratio + sum_wp += layer['theta_1500'] * contribution_ratio + break + + return sum_fc, sum_wp +def cal_stagebyDate(monitor_date,zwlx_growth_stages): + for stage, details in zwlx_growth_stages.items(): + time_range_str = details['time_range'] + if type(monitor_date) is str: + monitor_date= datetime.strptime(monitor_date, '%Y-%m-%d') + if '-' in time_range_str: + start_str, end_str = time_range_str.split('至') + start_date = datetime.strptime(start_str, '%Y-%m-%d') + end_date = datetime.strptime(end_str, '%Y-%m-%d') + if start_date <= monitor_date <= end_date: + return stage,details['root_length'],details['hight'],details['k_rain'] + raise ValueError(f"监测日期 {monitor_date.strftime('%Y-%m-%d')} 不在任何生育阶段范围内") +#根据日期计算主根长度 +def calculate_average_root_length(root_length_dist): + return sum(length * prob for length, prob in root_length_dist.items()) +def cal_soil_fc_wp(): + # 土壤分层数据 (33kPa和1500kPa对应的含水率) + soil_data = [ + {'depth_range': (0, 5), 'theta_33': 0.296, 'theta_1500': 0.103}, + {'depth_range': (5, 15), 'theta_33': 0.284, 'theta_1500': 0.104}, + {'depth_range': (15, 30), 'theta_33': 0.286, 'theta_1500': 0.105}, + {'depth_range': (30, 60), 'theta_33': 0.290, 'theta_1500': 0.114}, + {'depth_range': (60, 100), 'theta_33': 0.286, 'theta_1500': 0.117}, + {'depth_range': (100, 200), 'theta_33': 0.284, 'theta_1500': 0.114} + ] + # theta_fc, theta_wp = calculate_weighted_soil_parameters( + # growth_stages_roots[stage], + # soil_data + # ) + theta_fc=0.2884 + theta_wp=0.1086 + return theta_fc,theta_wp + +if __name__ == "__main__": + # #cal_et0_List(datetime(2026, 6, 18),datetime(2026, 6, 22)) + et0_daily=[23,45,67,89] + # etc,dates,kc_daily=calculate_etc_List(et0_daily, datetime(2026, 6, 18),datetime(2026, 6, 22), get_daily_kc(growth_stages_roots)) + # calculate_etc_dayily(45, datetime(2026, 6, 18), datetime(2026, 6, 22), get_daily_kc(growth_stages)) + # df=export_etc_et0_date(et0_daily,etc,dates,kc_daily) + # print(df) + # cal_soil_fc_wp('initial') + # stage,details = cal_stagebyDate(datetime(2025, 7, 18)) + # print(details) + # rd=calculate_average_root_length(details) + # print(rd) \ No newline at end of file diff --git a/irrgiation/db_connect.py b/irrgiation/db_connect.py new file mode 100644 index 0000000..c580557 --- /dev/null +++ b/irrgiation/db_connect.py @@ -0,0 +1,199 @@ +import configparser +import os +from pathlib import Path +from typing import Optional, Union, List + +import pandas as pd +import psycopg2 +from dotenv import load_dotenv +from psycopg2 import sql +from psycopg2.extras import execute_batch +from sqlalchemy import create_engine, text + + +def get_config(): + config = configparser.ConfigParser() + try: + project_nifi = Path(__file__).parent.parent / 'config.nifi' + config.read(project_nifi) + # 确保所有必要键都存在 + required_keys = ['host', 'port', 'database', 'user', 'password','schema'] + for key in required_keys: + if key not in config['postgresql']: + raise ValueError(f"Missing required config key: {key}") + + return { + "host": config['postgresql']['host'], + "port": config['postgresql']['port'], + "database": config['postgresql']['database'], + "user": config['postgresql']['user'], + "password": config['postgresql']['password'], + "schema": config['postgresql']['schema'], + "tablename": config['postgresql']['tablename'] + } + except Exception as e: + print(f"配置读取错误: {e}") + # 返回默认配置或退出 + return None +def connect(): + """建立数据库连接""" + try: + config_params=get_config() + conn = psycopg2.connect(host=config_params['host'],port=config_params['port'],dbname=config_params['database'], + user=config_params['user'],password=config_params['password'] + ) + cursor = conn.cursor() + except Exception as e: + print(f"❌ 连接失败: {e}") + raise + +# 根据条件查询 +def query_postgresql_to_dataframe( + condition: Optional[dict] = None, + condition_operator: str = "AND", + columns: Union[str, List[str]] = "*",): + + # 连接数据库 + config_params = get_config() + conn = psycopg2.connect(host=config_params['host'], port=config_params['port'], dbname=config_params['database'], + user=config_params['user'], password=config_params['password'] + ) + try: + # 构建列选择部分 + if isinstance(columns, list): + columns_sql = sql.SQL(", ").join(sql.Identifier(col) for col in columns) + else: + columns_sql = sql.SQL(columns) + + # 构建基础查询 + query = sql.SQL("SELECT {} FROM {}").format( + columns_sql, + sql.Identifier(config_params['tablename']) + ) + + # 构建条件部分（参数化） + params = {} + if condition: + conditions = [] + for i, (key, value) in enumerate(condition.items()): + param_name = f"param_{i}" + conditions.append(sql.SQL("{} = %({})s").format( + sql.Identifier(key), + sql.SQL(param_name) + )) + params[param_name] = value + + where_clause = sql.SQL(" WHERE {}").format( + sql.SQL(f" {condition_operator} ").join(conditions) + ) + query = query + where_clause + print(query) + # 执行查询 + with conn.cursor() as cursor: + cursor.execute(query, params) + if cursor.description: + columns = [desc[0] for desc in cursor.description] + data = cursor.fetchall() + return pd.DataFrame(data, columns=columns) + return pd.DataFrame() + + finally: + conn.close() + +# 插入库表 +def dataframe_to_postgresql_batch(df,batch_size=1000): + """ + 使用execute_batch批量插入DataFrame数据 + + 参数: + df: 要插入的DataFrame + table_name: 目标表名 + config_params: 数据库连接配置 + batch_size: 每批插入的行数 + """ + if df.empty: + print("DataFrame为空，无需插入") + return + + # 获取列名 + columns = df.columns.tolist() + data = [tuple(x) for x in df.to_numpy()] + + try: + config_params = get_config() + conn = psycopg2.connect(host=config_params['host'], port=config_params['port'], + dbname=config_params['database'], + user=config_params['user'], password=config_params['password'] + ) + + cursor = conn.cursor() + + # 构建INSERT语句 + insert_query = sql.SQL("INSERT INTO {} ({}) VALUES ({})").format( + sql.Identifier(config_params['tablename']), + sql.SQL(', ').join(map(sql.Identifier, columns)), + sql.SQL(', ').join([sql.Placeholder()] * len(columns)) + ) + + # 批量执行 + execute_batch(cursor, insert_query, data, batch_size) + + conn.commit() + + except Exception as e: + conn.rollback() + print(f"批量插入时出错: {e}") + finally: + if conn is not None: + conn.close() +# 根据条件更新某个字段 +def update_irrigation_data(date_value, id_value, field_to_update, new_value): + """ + 更新表中单条记录的单个字段 + + 参数: + db_url: 数据库连接字符串 + table_name: 表名 + date_field: 日期字段名 + id_field: ID/地块字段名 + date_value: 日期值 + id_value: ID/地块值 + field_to_update: 要更新的字段名 + new_value: 新值 + """ + try: + # 数据库配置 + db_config = { + "db_url": "postgresql://postgres:postgres@localhost:5432/datastore", + "table_name": "irrigation_data", + "date_field": "Date", + "id_field": "dkbm" + } + engine = create_engine(db_config["db_url"]) + table_name=db_config["table_name"] + date_field=db_config["date_field"] + id_field=db_config["id_field"] + with engine.begin() as conn: # 自动提交事务 + # 使用参数化查询防止SQL注入 + update_sql = text(f""" + UPDATE {table_name} + SET "{field_to_update}" = :new_value + WHERE "{date_field}" = :date_value AND "{id_field}" = :id_value""") + conn.execute(update_sql, { + 'new_value': new_value, + 'date_value': date_value, + 'id_value': id_value + }) + + return True + except Exception as e: + print(f"更新失败: {str(e)}") + return False + finally: + engine.dispose() # 确保连接关闭 + + +if __name__ == '__main__': + config = configparser.ConfigParser() + project_nifi = Path(__file__).parent.parent / 'config.nifi' + config.read(project_nifi) \ No newline at end of file diff --git a/irrgiation/irrigateDecisionDemo.py b/irrgiation/irrigateDecisionDemo.py new file mode 100644 index 0000000..4906775 --- /dev/null +++ b/irrgiation/irrigateDecisionDemo.py @@ -0,0 +1,620 @@ +import statistics +import uuid + +from collections import defaultdict +from datetime import datetime, timedelta + +import pandas as pd + +from irrgiation.db_connect import query_postgresql_to_dataframe, dataframe_to_postgresql_batch,update_irrigation_data + +from irrgiation.weatherAndSoilDataRequest import getWeatherAndSoilData, loginAuth, dateToTimestamp, get_soil_data_dk1, \ + getFutureWeather, weather_params, get_soil_data_dk2 +from irrgiation.dailyEvaporation import get_daily_kc, cal_et0_List, calculate_etc_List, export_etc_et0_date, \ + cal_soil_fc_wp, estimate_mad_from_et, cal_stagebyDate, calculate_average_root_length,calculate_etc_List_new +from irrgiation.mathuntils import day_of_year +from irrgiation.soil_Ke import calculate_ke,calculate_D_ei +from irrgiation.crop_growth_stages import crop_growth_stages + +#解析数据接口，计算均值、最大、最小值 +def analy_soil_data(data): + soil_result = {} + for key in data: + daily_data = defaultdict(lambda: defaultdict(list)) + for metric, entries in data[key].items(): + for entry in entries: + dt = datetime.strptime(entry["datetime"], "%Y-%m-%d %H:%M:%S") + date_key = dt.date() # 按日期分组 + value = float(entry["value"]) + daily_data[date_key][metric].append(value) + result = {} + for date, metrics in daily_data.items(): + date_stats = {} + for metric, values in metrics.items(): + date_stats[metric] = { + "max": max(values), + "min": min(values), + "mean": statistics.mean(values), + 'count': len(values), + 'sum': sum(values) + } + result[str(date)] = date_stats # 将日期转为字符串作为键 + soil_result[key] = result + return soil_result + + +#获取地块6个传感器的土壤湿度均值 +def analy_soil_data_mean(soil_data): + result = defaultdict(lambda: { + "SOIL_MOISTURE_SURFACE": [], + "SOIL_MOISTURE_MIDDLE": [], + "SOIL_MOISTURE_BOTTOM": [] + + }) + # 填充数据 + for sensor, dates in soil_data.items(): + for date, layers in dates.items(): + result[date]["SOIL_MOISTURE_SURFACE"].append(layers["SOIL_MOISTURE_SURFACE"]["mean"]) + result[date]["SOIL_MOISTURE_MIDDLE"].append(layers["SOIL_MOISTURE_MIDDLE"]["mean"]) + result[date]["SOIL_MOISTURE_BOTTOM"].append(layers["SOIL_MOISTURE_BOTTOM"]["mean"]) + # 转换为普通字典 + final_result = dict(result) + result_soil = {} + for date, layers in final_result.items(): + result_soil[date] = { + "SOIL_MOISTURE_SURFACE_MEAN": sum(layers["SOIL_MOISTURE_SURFACE"]) / len(layers["SOIL_MOISTURE_SURFACE"]), + "SOIL_MOISTURE_MIDDLE_MEAN": sum(layers["SOIL_MOISTURE_MIDDLE"]) / len(layers["SOIL_MOISTURE_MIDDLE"]), + "SOIL_MOISTURE_BOTTOM_MEAN": sum(layers["SOIL_MOISTURE_BOTTOM"]) / len(layers["SOIL_MOISTURE_BOTTOM"]) + } + + return final_result, result_soil + + +def analy_soil_dataNight_mean(soil_data): + daily_means = defaultdict(lambda: defaultdict(list)) + # 计算平均值 + for sensor, measurements in soil_data.items(): + for measure_type, dates in measurements.items(): + for date, value in dates.items(): + daily_means[date][measure_type].append(float(value)) + + # 计算并格式化结果 + result = {} + for date, measures in daily_means.items(): + result[date] = { + f"{measure_type}_MEAN": round(sum(values) / len(values), 2) + for measure_type, values in measures.items() + } + + return result + + +#获取地块6个传感器中土壤湿度最大值 +def analy_soil_data_max(soil_data): + result = defaultdict(lambda: { + "SOIL_MOISTURE_SURFACE": [], + "SOIL_MOISTURE_MIDDLE": [], + "SOIL_MOISTURE_BOTTOM": [] + + }) + # 填充数据 + for sensor, dates in soil_data.items(): + for date, layers in dates.items(): + result[date]["SOIL_MOISTURE_SURFACE"].append(layers["SOIL_MOISTURE_SURFACE"]["mean"]) + result[date]["SOIL_MOISTURE_MIDDLE"].append(layers["SOIL_MOISTURE_MIDDLE"]["mean"]) + result[date]["SOIL_MOISTURE_BOTTOM"].append(layers["SOIL_MOISTURE_BOTTOM"]["mean"]) + # 转换为普通字典 + final_result = dict(result) + result_soil = {} + for date, layers in final_result.items(): + result_soil[date] = { + "SOIL_MOISTURE_SURFACE_MEAN": max(layers["SOIL_MOISTURE_SURFACE"]), + "SOIL_MOISTURE_MIDDLE_MEAN": max(layers["SOIL_MOISTURE_MIDDLE"]), + "SOIL_MOISTURE_BOTTOM_MEAN": max(layers["SOIL_MOISTURE_BOTTOM"]) + } + + return final_result, result_soil + + +def get_soil_data_night_eyveryDay(data): + daily_data = {} + for sersor in data: + data_depth = {} + for key in data[sersor]: + target_time = (21, 59, 59) + result = {} # 时、分、秒元组 + daily_records = {} # 临时存储每天的所有记录，用于找不到22点数据时获取最新记录 + for record in data[sersor][key]: + try: + dt = datetime.strptime(record["datetime"], "%Y-%m-%d %H:%M:%S") + date_str = record["datetime"].split(" ")[0] # 提取日期部分 + # 严格匹配 02:00:00 + # 提取时间部分（时、分、秒） + current_time = (dt.hour, dt.minute, dt.second) + # 存储每天的所有记录（按时间排序） + if date_str not in daily_records: + daily_records[date_str] = [] + daily_records[date_str].append((dt, record['value'])) + # 比较时间部分 + if current_time >= target_time: + if date_str not in result: + result[date_str] = record['value'] + except: + continue + # 对于没有22点数据的日期，取当天最新的一条数据 + for date_str in daily_records: + if date_str not in result and daily_records[date_str]: + daily_records[date_str].sort(key=lambda x: x[0]) # 按datetime排序 + result[date_str] = daily_records[date_str][-1][1] # 取最后一个值 + data_depth[key] = result + daily_data[sersor] = data_depth + return daily_data + + +def analyze_weather_data1(data): + # 1. 按日期分组数据 + daily_data = defaultdict(lambda: defaultdict(list)) + for metric, entries in data.items(): + for entry in entries: + dt = datetime.strptime(entry["datetime"], "%Y-%m-%d %H:%M:%S") + date_key = dt.date() # 按日期分组 + if metric == "AIR_LUX": + value = float(entry["value_solar"]) + else: + value = float(entry["value"]) + daily_data[date_key][metric].append(value) + # 2. 计算每日统计值 + result = {} + for date, metrics in daily_data.items(): + date_stats = {} + for metric, values in metrics.items(): + date_stats[metric] = { + "max": max(values), + "min": min(values), + "mean": statistics.mean(values), + 'count': len(values), + 'sum': sum(values) + } + result[str(date)] = date_stats # 将日期转为字符串作为键 + + return result + + +# 获取彭曼公式的入参 +def get_params(start_date, monitor_date, data): + delta = monitor_date - start_date + days = delta.days + dates = [start_date + timedelta(days=i) for i in range(days)] + date_strings = [date.strftime("%Y-%m-%d") for date in dates] + date_params = {} + for date in date_strings: + result_params = get_weather_params(data, date) + date_params[date] = result_params + return date_params + + +def get_weather_params(data, date): + result = {} + elevation = 865 # 海拔(m) + lat = 41.34252110189814 # 纬度(度) + #计算当前日期位于当年的多少天 + doy = day_of_year(date) + + for item in data['AIR_LUX']: + item['value_solar'] = float(item['value']) * 0.0081 * 3600 + data_tj = analyze_weather_data1(data) + data_tj_date = data_tj[date]; + tmax = data_tj_date['AIR_TEMPERATURE']['max'] + tmin = data_tj_date['AIR_TEMPERATURE']['min'] + tmean = data_tj_date['AIR_TEMPERATURE']['mean'] + rh_mean = data_tj_date['AIR_HUMIDITY']['mean'] + wind_speed = data_tj_date['WIND_SPEED']['mean'] + solar = data_tj_date['AIR_LUX']['sum'] / 1000000 + result['tmax'] = round(tmax, 3) + result['tmin'] = round(tmin, 3) + result['tmean'] = round(tmean, 3) + result['rh_mean'] = round(rh_mean, 3) + result['wind_speed'] = round(wind_speed, 3)*0.94 + result['solar_rad'] = round(solar, 3) + result['elevation'] = elevation + result['lat'] = round(lat, 3) + result['doy'] = doy + return result + + +# 判断是否需要灌溉 +def isNeedirrigate(df, qx_coff, soil_coff,Di_1,zwlx_growth_stages): + df['RD'] = 0.0 + df['Di'] = 0.0 + df['Di-1'] = Di_1 + df['mad'] = 0.0 + df['irrigation_depth'] = 0.0 + df['irrigation_acre'] = 0.0 + df['isNeedirrigate'] = False + df['Reset_Flag'] = False + # df['irrigate_sd_predict'] = 0.0 + df['f_t'] = 0.0 + df['f_js'] = 0.0 + df['f_soil'] = 0.0 + for i in range(len(df)): + # 设置Di-1（前一天累计值） + if i > 0: + df.at[i, 'Di-1'] = df.at[i - 1, 'Di'] + # 计算当前累计值（Di-1 + 当天ETc） + raw, rd, mad, k_rain = cal_raw(df.at[i, 'ETc'], df.at[i, 'Date'],zwlx_growth_stages) + if df.at[i, 'rain_effective']<5: + df.at[i, 'rain_effective'] = 0 + else: + df.at[i, 'rain_effective'] = df.at[i, 'rain_effective'] * k_rain + current_cum = df.at[i, 'Di-1'] + df.at[i, 'ETc']- df.at[i, 'rain_effective'] + df['RD'] = rd + df.at[i,'mad'] = mad + df.at[i, 'raw'] = round(raw, 3) + f_t = qx_coff[df.at[i, 'Date']]['qw'] + f_js = qx_coff[df.at[i, 'Date']]['js'] + f_soil = soil_coff[df.at[i, 'Date']] + raw_adj = cal_raw_adjust(raw, f_t, f_js, f_soil) + df.at[i, 'f_t'] = f_t + df.at[i, 'f_js'] = f_js + df.at[i, 'f_soil'] = f_soil + df.at[i, 'raw_adjust'] = round(raw_adj, 3) + # 判断是否需要重置 + if current_cum > raw_adj: + df.at[i, 'Reset_Flag'] = True + df.at[i, 'irrigation_depth'] = round(current_cum, 3) + df.at[i, 'isNeedirrigate'] = True + df.at[i, 'Di'] = 0 # 保留当天ETc值 + df.at[i, 'irrigation_acre'] = round((current_cum*667/1000),3) + else: + df.at[i, 'Di'] = round(current_cum, 3) # 正常累加 + return df + + +def cal_raw(etc, monitor_date,zwlx_growth_stages): + theta_fc, theta_wp = cal_soil_fc_wp() + stage, details,hight,k_rain= cal_stagebyDate(monitor_date,zwlx_growth_stages) + rd = calculate_average_root_length(details) + mad = estimate_mad_from_et(etc) + raw = mad * (theta_fc - theta_wp) * rd * 10 + return raw, rd,mad,k_rain + + +def cal_raw_adjust(raw, f_t, f_rain, f_soil): + raw_adjust = raw * f_soil * f_t * f_rain + return raw_adjust + + +def cal_qxcoefficient_adjust(start_date, monitor_date, qx_data): + + qx_coefficient = {} + current_date = start_date + qxData_filter = {} + for key in qx_data: + if key == "AIR_TEMPERATURE" or key == "RAIN_FALL_REALTIME": + qxData_filter[key] = qx_data[key] + qxData_filter_tj = analyze_weather_data1(qxData_filter) + now_date = datetime.now().date() + while current_date < monitor_date: + qx_futureData = [] + result = {} + if (current_date + timedelta(days=3)).date() < now_date: + start = datetime.strptime((current_date + timedelta(days=1)).strftime("%Y-%m-%d"), "%Y-%m-%d") + end = datetime.strptime((current_date + timedelta(days=3)).strftime("%Y-%m-%d"), "%Y-%m-%d") + # 筛选日期范围内的数据 + filter = { + date: values + for date, values in qxData_filter_tj.items() + if start <= datetime.strptime(date, "%Y-%m-%d") <= end + } + qx_futureData = transform_data(filter) + else: + try: + # print(f"进入else分支: {current_date}") + qx_futureData = getFutureWeather(current_date) + except Exception as e: + # print(f"else分支异常: {e}") + raise # 保留异常以便调试 + qw_mean = (qx_futureData['qw_day1']['mean'] + qx_futureData['qw_day2']['mean'] + qx_futureData['qw_day3']['mean']) / 3 + js_sum = qx_futureData['js_day1']['max'] + qx_futureData['js_day2']['max'] + qx_futureData['js_day3']['max'] + qw_coff = 1.0 + js_coff = 1.0 + if qw_mean < 10: + qw_coff = 1.1 + elif 10 < qw_mean <= 25: + qw_coff = 1 + elif 25 < qw_mean <= 30: + qw_coff = 0.85 + else: + qw_coff = 0.75 + if js_sum <= 5: + js_coff = 1.0 # 少雨 + elif 5 < js_sum <= 10: + js_coff = 1.2 # 中雨 + elif 10 < js_sum <= 20: + js_coff = 1.4 # 大雨 + else: + js_coff = 1.6 # 暴雨 + result['qw'] = qw_coff + result['js'] = js_coff + qx_coefficient[current_date.strftime("%Y-%m-%d")] = result + current_date = current_date + timedelta(days=1) + return qx_coefficient + + +def transform_data(original_data): + transformed = {} + # 处理温度数据 (qw_day) + for i, (date, metrics) in enumerate(original_data.items(), 1): + temp_data = metrics["AIR_TEMPERATURE"] + transformed[f"qw_day{i}"] = { + "max": round(temp_data["max"], 2), + "min": round(temp_data["min"], 2), + "mean": round(temp_data["mean"], 2), + "sum": round(temp_data["sum"], 2) + } + + # 处理降雨数据 (js_day) + for i, (date, metrics) in enumerate(original_data.items(), 1): + rain_data = metrics["RAIN_FALL_REALTIME"] + transformed[f"js_day{i}"] = { + "max": round(rain_data["max"], 2), + "min": round(rain_data["min"], 2), + "mean": round(rain_data["mean"], 2), + "sum": round(rain_data["sum"], 2) + } + + return transformed + + +def cal_soil_coff(soil_data): + soil_data_coff = {} + for sensor in soil_data: + for date in soil_data[sensor]: + soil_cof = 1.0 + sd_surf = float(soil_data['SOIL_MOISTURE_SURFACE'][date]) + sd_mid = float(soil_data['SOIL_MOISTURE_MIDDLE'][date]) + if sd_surf < 10: + soil_cof = 0.5 + if 10 0 and rain < 4: + fw = 0.35 + elif I > 0 and rain < 4: + fw = 1 + elif rain > 4 and I == 0: + fw = 1 + else: + fw = ke_data[prev_date_key]['fw'] + et0 = 0.0 + kc = 0.0 + if i < len(et0_list): + et0 = et0_list[i] + kc = kc_list[i] + ke,few=calculate_ke(hight, u2, rh_min, kc, DE_i_1, REW, fw, theta_fc,theta_wp) + et_soil=ke*et0 + DE_i = calculate_D_ei(DE_i_1, rain, I,et_soil, fw,few) + #创建当天的数据字典 + daily_data = { + 'ke': ke, + 'DE_i': DE_i, + 'et0': et0, + 'et_soil': et_soil, + 'fw':fw, + "rain":rain, + "irrigation":I + } + ke_data[date_key]=daily_data + current_date = current_date + timedelta(days=1) + i=i+1 + return ke_data +# 单天监测 +def cal_irrigationByDay(zwlx_name,dkname,start_date,end_date,soil_key,weather_key,irrigation_really): + yesterday_date=(start_date - timedelta(days=1)).strftime("%Y-%m-%d") + condition = {"Date": yesterday_date,"dkbm":dkname} + df=query_postgresql_to_dataframe(condition) + Di_1 = 0 + if irrigation_really>0: + Di_1 = 0 + update_irrigation_data(yesterday_date,dkname,"irrigation_really",irrigation_really) + else: + # 将Date列转换为日期时间类型 + df['Date'] = pd.to_datetime(df['Date']) + if len(df[df.Date == (start_date - timedelta(days=1))])!=0: + isNeedjg=df[df.Date == (start_date- timedelta(days=1))]['isNeedirrigate'].values[0] + if isNeedjg: + Di_1 = df[df.Date == (start_date - timedelta(days=1))]['Di-1'].values[0]+df[df.Date == (start_date - timedelta(days=1))]['ETc'].values[0] + else: + Di_1 = df[df.Date == (start_date - timedelta(days=1))]['Di'].values[0] + dataframe=cal_dk_isNeegirrigate_Day(zwlx_name,dkname,soil_key,start_date,end_date,Di_1,weather_key) + condition_new = {"Date": start_date.strftime("%Y-%m-%d"), "dkbm": dkname} + df_query=query_postgresql_to_dataframe(condition_new) + if len(df_query)==0: + dataframe_to_postgresql_batch(dataframe) + else: + print("库表中已存在该数据") + return dataframe +if __name__ == "__main__": + #调用气象数据接口--获取具体日期的气象数据(气象) + # auth_token = loginAuth() + start_date = datetime(2025, 7, 16) + monitor_date = datetime(2025, 7, 17) + dk_name = "5" + zwlx_name="甘草" + soil_key = "8637b970-561b-11f0-a556-4f10f26fc07f" + weather_key="18d121f0-561b-11f0-a556-4f10f26fc07f" + irrigation_really=5 + cal_irrigationByDay(zwlx_name,dk_name,start_date,monitor_date,soil_key,weather_key,irrigation_really) + # cal_dk_isNeegirrigate(dk_name, planting_date, monitor_date,0) diff --git a/irrgiation/mathuntils.py b/irrgiation/mathuntils.py new file mode 100644 index 0000000..d9c981b --- /dev/null +++ b/irrgiation/mathuntils.py @@ -0,0 +1,371 @@ +#根据日期计算当年天数 +import numpy as np +import pandas as pd +from matplotlib import pyplot as plt + + +def is_leap_year(year): + return (year % 4 == 0 and year % 100 != 0) or (year % 400 == 0) + +def day_of_year(date_str): + date_part = date_str.split()[0] + year = int(date_str[:4]) # 取前4位 → 2025 + month = int(date_str[5:7]) # 取第5-6位 → 06 → 6 + day = int(date_str[8:10]) # 取第8-9位 → 30 + month_days = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] + if is_leap_year(year): + month_days[1] = 29 # 闰年2月29天 + total = sum(month_days[:month-1]) + day + return total +#计算每天有效日照时长 +def count_above_threshold(data, threshold,date): + filtered = filter(lambda x: float(x['value_cal']) >= threshold and x['datetime'].split(" ")[0] == date, data) + return len(list(filtered)) +# 示例计算函数：Lux转W/m² +def lux_to_wm2(lux_value): + """将Lux值转换为W/m²（使用标准系数0.0081）""" + return lux_value * 0.0081 +def convert_lux_to_Radiant_Flux(data_dict): + """ + 递归地将字典中的所有'timestamp'字段转换为可读的日期时间格式 + 保留原始时间戳的同时添加新的'datetime'字段 + """ + if 'value' in data_dict: + radiant_Flux=data_dict["value"]*0.001496 + data_dict["value_radiant"] =radiant_Flux + return data_dict +def lux_to_solar_rad(lux,h): + + return lux * 0.001496*h * 3.6 / 1000 +def plot(soil_data,name): + import matplotlib.dates as mdates + # 创建画布 + plt.figure(figsize=(14, 8)) + + # 处理并绘制地表湿度数据 + surface_df = pd.DataFrame(soil_data["SOIL_MOISTURE_SURFACE"]) + surface_df['datetime'] = pd.to_datetime(surface_df['datetime']) + surface_df['value'] = surface_df['value'].astype(float) + plt.plot(surface_df['datetime'], surface_df['value'], + label='Surface Moisture (0-5cm)', color='#FF6B6B', linewidth=2, marker='o', markersize=5) + + # 处理并绘制中层湿度数据 + middle_df = pd.DataFrame(soil_data["SOIL_MOISTURE_MIDDLE"]) + middle_df['datetime'] = pd.to_datetime(middle_df['datetime']) + middle_df['value'] = middle_df['value'].astype(float) + plt.plot(middle_df['datetime'], middle_df['value'], + label='Middle Layer Moisture', color='#4ECDC4', linewidth=2, marker='s', markersize=5) + # 处理并绘制中层湿度数据 + middle_df = pd.DataFrame(soil_data["SOIL_MOISTURE_BOTTOM"]) + middle_df['datetime'] = pd.to_datetime(middle_df['datetime']) + middle_df['value'] = middle_df['value'].astype(float) + plt.plot(middle_df['datetime'], middle_df['value'], + label='BOTTOM Layer Moisture', color='blue', linewidth=2, marker='s', markersize=5) + + # 设置图表格式 + plt.title('Soil Moisture Variation (Jun 18 - Jul 3, 2025)', fontsize=16, pad=20) + plt.xlabel('Date', fontsize=12) + plt.ylabel('Soil Moisture (%)', fontsize=12) + plt.grid(True, linestyle='--', alpha=0.7) + + # 设置x轴日期格式 + ax = plt.gca() + ax.xaxis.set_major_locator(mdates.DayLocator(interval=2)) + ax.xaxis.set_major_formatter(mdates.DateFormatter('%m/%d')) + plt.xticks(rotation=45) + + # 添加图例和注释 + plt.legend(loc='upper right', fontsize=10) + plt.tight_layout() + # 保存图片到本地 + output_path = 'D:\悦柠\遥感\甘草\ceshi/'+name+".png" + plt.savefig(output_path, dpi=300, bbox_inches='tight') + # 显示图表 + # plt.show() + + +def rh_to_theta_vg(rh, soil_params, temperature=25): + """使用Van Genuchten模型将相对湿度转换为体积含水量""" + # 计算土壤水势(kPa) + psi = -145 * np.log(rh / 100) + + # 提取土壤参数 + theta_s = soil_params["theta_s"] # 饱和含水量 + theta_r = soil_params["theta_r"] # 残余含水量 + alpha = soil_params["alpha"] # Van Genuchten参数 + n = soil_params["n"] # Van Genuchten参数 + m = 1 - 1 / n + + # 计算体积含水量 + theta = theta_r + (theta_s - theta_r) / ((1 + (alpha * abs(psi)) ** n) ** m) + print(theta) + return theta + + +def calculate_soil_moisture_after_irrigation( + initial_moisture: float, + field_capacity: float, + irrigation_volume: float, + soil_volume: float, + area: float = None, + verbose: bool = False +) -> float: + """ + 不考虑降雨量和蒸散量时，计算灌溉后的土壤相对湿度 + + 参数: + initial_moisture: 初始土壤相对湿度 (%) + field_capacity: 田间持水量 (%) + irrigation_volume: 灌溉量 (m³) + soil_volume: 土壤体积 (m³) + area: 灌溉面积 (㎡，可选，仅用于打印提示) + verbose: 是否打印详细计算过程 + + 返回: + float: 灌溉后土壤相对湿度 (%) + """ + # 参数验证 + if not (0 <= initial_moisture <= 100): + raise ValueError("初始湿度必须在0-100%之间") + if not (0 < field_capacity <= 100): + raise ValueError("田间持水量必须在0-100%之间且大于0") + if irrigation_volume < 0: + raise ValueError("灌溉量不能为负数") + if soil_volume <= 0: + raise ValueError("土壤体积必须大于0") + + # 计算初始土壤含水量（m³） + initial_water = (initial_moisture / 100) * (field_capacity / 100) * soil_volume + + # 灌溉后总含水量（m³） + total_water = initial_water + irrigation_volume + + # 计算田间持水总量（m³） + max_water = (field_capacity / 100) * soil_volume + + # 计算灌溉后湿度（超过田间持水量时取100%） + final_moisture = min((total_water / max_water) * 100, 100.0) + + # if verbose: + # area_info = f"（面积: {area}㎡）" if area else "" + # print(f"===== 不考虑降雨与蒸散的灌溉计算 {area_info} =====") + # print(f"初始湿度: {initial_moisture}%") + # print(f"田间持水量: {field_capacity}%") + # print(f"土壤体积: {soil_volume}m³") + # print(f"灌溉量: {irrigation_volume}m³") + # print(f"初始含水量: {initial_water:.2f}m³") + # print(f"灌溉后含水量: {total_water:.2f}m³") + # print(f"田间最大持水量: {max_water:.2f}m³") + # print(f"灌溉后土壤湿度: {final_moisture:.2f}%") + + return final_moisture +def calculate_irrigation_volume( + initial_moisture: float, + target_moisture: float, + field_capacity: float, + soil_volume: float +) -> float: + """ + 计算为达到目标土壤相对湿度所需的灌溉量 + + 参数: + initial_moisture: 初始土壤相对湿度 (%) + target_moisture: 目标土壤相对湿度 (%) + field_capacity: 田间持水量 (%) + soil_volume: 土壤体积 (m³) + area: 灌溉面积 (㎡，可选，仅用于打印提示) + verbose: 是否打印详细计算过程 + + 返回: + float: 需要的灌溉量 (m³) + """ + # 参数验证 + if not (0 <= initial_moisture <= 100): + raise ValueError("初始湿度必须在0-100%之间") + if not (0 <= target_moisture <= 100): + raise ValueError("目标湿度必须在0-100%之间") + if not (0 < field_capacity <= 100): + raise ValueError("田间持水量必须在0-100%之间且大于0") + if soil_volume <= 0: + raise ValueError("土壤体积必须大于0") + if target_moisture <= initial_moisture: + raise ValueError("目标湿度必须大于初始湿度") + + # 计算初始土壤含水量（m³） + initial_water = (initial_moisture / 100) * (field_capacity / 100) * soil_volume + + # 计算目标含水量（m³） + target_water = (target_moisture / 100) * (field_capacity / 100) * soil_volume + + # 需要的灌溉量 + irrigation_volume = target_water - initial_water + + return irrigation_volume +def irriDepth_irrvolume(area,irrDepth): + volume=area*667*(irrDepth/1000) + return volume + + +def calculate_irrigation(dry_soil_weight, current_humidity, target_humidity, + irrigation_efficiency=0.8, water_density=1000): + """ + 计算达到目标湿度所需的灌溉量 + + Args: + dry_soil_weight (float): 土壤干重（kg） + current_humidity (float): 当前相对湿度（%） + target_humidity (float): 目标相对湿度（%） + irrigation_efficiency (float): 灌溉效率（0~1） + water_density (float): 水密度（kg/m³） + + Returns: + float: 灌溉量（立方米） + """ + if current_humidity >= target_humidity: + return 0.0 + + delta_humidity = (target_humidity - current_humidity) / 100 + water_weight = dry_soil_weight * delta_humidity + water_volume = water_weight / water_density + return water_volume / irrigation_efficiency +def calculate_rh_target( + initial_moisture: float, + root: float, + field_capacity: float, + irrigation:float +) -> float: + initial_moisture=(initial_moisture / 100) * field_capacity + iir_volume=irrigation/(10*root) + target_moisture = initial_moisture + (iir_volume/field_capacity) + return target_moisture + + +def calculate_rh_target( + initial_moisture: float, # 初始土壤相对湿度(%) + root: float, # 根系深度(cm) + field_capacity: float, # 田间持水量(体积含水量 cm³/cm³) + irrigation: float # 灌溉量(mm) +) -> float: + """ + 计算灌溉后的目标土壤相对湿度(%) + + 参数: + initial_moisture: 初始土壤相对湿度(%) + root: 根系深度(cm) + field_capacity: 田间持水量(体积含水量 cm³/cm³) + irrigation: 灌溉量(mm) + + 返回: + float: 灌溉后的目标土壤相对湿度(%) + """ + # 参数验证 + if initial_moisture < 0 or initial_moisture > 100: + raise ValueError("初始湿度必须在0-100%之间") + if root <= 0: + raise ValueError("根系深度必须大于0") + if field_capacity <= 0 or field_capacity > 1: + raise ValueError("田间持水量必须在0-1 cm³/cm³之间") + if irrigation < 0: + raise ValueError("灌溉量不能为负值") + + # 将相对湿度转为体积含水量 + initial_volume_moisture = (initial_moisture / 100) * field_capacity + + # 计算灌溉增加的体积含水量 + # 将mm转为cm，并除以根系深度得到体积含水量增量 + moisture_increase = (irrigation / 10) / root + + # 计算灌溉后的体积含水量 + target_volume_moisture = initial_volume_moisture + moisture_increase + # 转换回相对湿度(%) + target_rh = (target_volume_moisture / field_capacity) * 100 + + # 确保不超过100% + target_rh = min(target_rh, 100.0) + + # # 3. 计算目标体积含水量 (不超过田间持水量) + # target_volume = min(initial_volume + moisture_increase, field_capacity) + # + # # 4. 转换回相对湿度 (%) + # target_rh = (target_volume / field_capacity) * 100 + + return target_rh + + +def calculate_rh_target_rz( + initial_moisture: float, # 初始土壤相对湿度(%) + soil_root: float, # 根系深度(cm) + area: float, # 灌溉面积(亩) + irrigation: float, # 灌溉量(mm) + field_capacity: float = 28.84, # 田间持水量(%)，默认35% + soil_bulk_density: float = 1.5 # 土壤容重(g/cm³)，默认1.5 +) -> float: + """ + 计算灌溉后的目标根区土壤相对湿度 + + 参数: + initial_moisture: 初始土壤相对湿度(%) + soil_root: 根系深度(cm) + area: 灌溉面积(亩) + irrigation: 灌溉量(mm) + field_capacity: 田间持水量(%)，默认35% + soil_bulk_density: 土壤容重(g/cm³)，默认1.5 + + 返回: + float: 灌溉后的目标土壤相对湿度(%) + """ + # 参数验证 + if initial_moisture < 0 or initial_moisture > 100: + raise ValueError("初始湿度必须在0-100%之间") + if soil_root <= 0: + raise ValueError("根系深度必须大于0") + if area <= 0: + raise ValueError("面积必须大于0") + if irrigation < 0: + raise ValueError("灌溉量不能为负值") + + # 计算根区体积 + root_volume = area * 667 * (soil_root / 100) # m³，将cm转为m + + # 计算灌溉水体积 + water_volume = (irrigation / 1000) * area * 667 # m³，将mm转为m + + # 计算灌溉增加的土壤湿度(%) + moisture_increase = (water_volume / (root_volume * soil_bulk_density)) * 100 + + # 计算目标湿度 + target_moisture = initial_moisture + moisture_increase + + # 确保不超过田间持水量 + target_moisture = min(target_moisture, field_capacity) + + return target_moisture +if __name__ == '__main__': + # volume=irriDepth_irrvolume(10.95,11.123) + # soil_volume_dk1 = 10.95 * 667 * 0.15 + # soil_volume_dk2 = 8.92 * 667 * 0.15 + # result1 = calculate_soil_moisture_after_irrigation( + # initial_moisture=15.2, + # field_capacity=28.84, + # irrigation_volume=168, + # soil_volume=soil_volume_dk1, + # verbose=True + # ) + # print(f"\n示例1结果: 灌溉后湿度 = {result1:.2f}%\n") + # print(soil_volume_dk1) + # result2=calculate_irrigation_volume(initial_moisture=15.44, + # field_capacity=28.84, + # target_moisture=70, + # soil_volume=soil_volume_dk1) + # print(f"\n达到目标湿度所需灌溉量 = {result2:.2f}\n") + # # 调用示例 + # dry_soil_weight=10.95*667*0.1*1200 + # volume = calculate_irrigation( + # dry_soil_weight=dry_soil_weight, + # current_humidity=15.4, + # target_humidity=70 + # ) + # print(f"灌溉量: {volume:.2f} m³") # 输出: 0.09 m³ + rh=calculate_rh_target(67.96,6,0.2884,33.14) + rh1 = rh*0.7 + print(rh,rh1) diff --git a/irrgiation/soil_Ke.py b/irrgiation/soil_Ke.py new file mode 100644 index 0000000..c320f93 --- /dev/null +++ b/irrgiation/soil_Ke.py @@ -0,0 +1,133 @@ +# 计算土壤蒸发系数 +def calculate_ke(h, u2, rh_min, k_cb, DE_i_1,REW,fw,theta_fc,theta_wp): + """ + 计算土壤蒸发系数 ke + 参数: + h: 株高(m) + u2: 2米高处风速(m/s) + rh_min: 最小相对湿度(%) + k_cb: 基础作物系数 + DE_i_1: 前一天从土壤表层蒸发的累积深度(mm) + TEW: 总可蒸发水量(mm) + REW: 易蒸发水量(mm) + few: 裸露湿润的土壤的比值 + 返回: + ke: 土壤蒸发系数 + """ + K_cmin = 0.175 + Z_e=0.12 + kc_max = calculate_kc_max(h, u2, rh_min, k_cb) + fc = calculate_fc(k_cb, K_cmin, kc_max, h) + few = calculate_fev(fc, fw) + TEW = calculate_tew(theta_fc, theta_wp, Z_e) + kr = calculate_Kr(DE_i_1, TEW, REW) + term1 = kr * (kc_max - k_cb) + term2 = few * k_cb + return min(term1, term2),few + +# 每日土壤蒸发量累积量 +#土壤蒸发有效部分 +def calculate_fev(fc, fw): + fev = min(1 - fc, fw) + return fev +#计算作物覆盖率fc +def calculate_fc(Kcb, Kcmin, Kcmax, h): + #Kcb：作物基础系数，Kcmin：无地表覆盖的干燥土壤最小 值（约为 0.15－0.20），Kcmax：紧随湿润过程的土壤最大 值 + + fc = ((Kcb- Kcmin) / (Kcmax - Kcmin)) ** (1 + 0.5 * h) + return fc +# 计算蒸发减少系数 +def calculate_Kr(DE_i_1, TEW, REW): + if DE_i_1 <= REW: + Kr = 1 + else: + print("dd") + Kr = (TEW - DE_i_1) / (TEW - REW) + return Kr +#计算Kcmax +def calculate_kc_max(h: float, u2: float, rh_min: float, k_cb: float) -> float: + """ + 计算作物系数上限 kc_max + 参数: + h: 株高(m) + u2: 2米高处风速(m/s) + rh_min: 最小相对湿度(%) + k_cb: 基础作物系数 + 返回: + kc_max: 作物系数上限 + """ + term1 = 1.2 + (0.04 * (u2 - 2) - 0.004 * (rh_min - 45)) * (h / 3) ** 0.3 + term2 = k_cb + 0.5 + return max(term1, term2) +#计算TEW +def calculate_tew(theta_fc: float, theta_wp: float, Z_e: float) -> float: + """ + 计算总可蒸发水量 TEW + 参数: + theta_fc: 田间持水量(m³/m³) + theta_wp: 萎蔫系数(m³/m³) + Z_e: 土壤蒸发层深度(cm) + 返回: + TEW: 总可蒸发水量(mm) + """ + return 1000 * (theta_fc - theta_wp) * Z_e +# 计算深层渗漏损失 +def calculate_DP_ei(P_i, RO_i, I_i, D_e_prev, f_w): + if (P_i - RO_i) + I_i / f_w - D_e_prev >= 0: + return (P_i - RO_i) + I_i / f_w - D_e_prev + else: + return 0 + +# 计算每日土壤蒸发量累积量 +def calculate_D_ei(D_e_prev, P_i, I_i, E_i, f_w, f_ew): + RO_i=0 + T_ew_i=0 + DP_ei = calculate_DP_ei(P_i, RO_i, I_i, D_e_prev, f_w) + D_ei = D_e_prev - (P_i - RO_i) - I_i / f_w + E_i / f_ew + T_ew_i + DP_ei + return D_ei +# 示例计算 +if __name__ == "__main__": + # # 示例参数 + # h = 0.07 # 株高(m) + # u2 = 3.0 # 2米高处风速(m/s) + # rh_min = 20.0 # 最小相对湿度(%) + # k_cb = 0.4 # 基础作物系数 + # K_s = 0.6 # 土壤导水率 + # K_cmin = 0.175 + # + # fw = 0.35 # 土壤表层被灌溉浸湿的部分 + # DE_i = 6.24 # 前一天从土壤表层蒸发的累积深度(mm) + # DE_i_1 = 2.965 # 前一天从土壤表层蒸发的累积深度(mm) + # theta_fc = 0.2884 # 田间持水量(m³/m³) + # theta_wp = 0.1086 # 萎蔫系数(m³/m³) + # Z_e = 0.1 # 土壤蒸发层深度(m) + # + # # # # 计算 TEW + # # TEW = calculate_tew(theta_fc, theta_wp, Z_e) + # REW = 5.178 # 假设REW为TEW的40% + # # # + # # # # 计算各参数 + # # kc_max = calculate_kc_max(h, u2, rh_min, k_cb) + # # fc = calculate_fc(k_cb, K_cmin, kc_max, h) + # # few= calculate_fev(fc, fw) + # # kr = calculate_Kr(DE_i_1,DE_i, TEW, REW) + # ke=calculate_ke(h, u2, rh_min, k_cb, DE_i_1,DE_i,REW) + # # 输出结果 + # + # print(f"土壤蒸发系数 Ke: {ke:.4f}") + # 假设初始参数 + TEW = 1000 # 总可蒸发水量，单位：mm + D_e_start = TEW # 模拟开始时的初始D_e + D_e_prev = 0 # 模拟第一天的D_e，假设充分灌溉后为0 + # 假设一些每日数据，实际应用中需替换为真实数据 + P_i = 10 # 第i天的降雨量，单位：mm + RO_i = 0 # 第i天的降雨形成地表径流量，这里假定为0，单位：mm + I_i = 19 # 第i天渗入土壤的灌溉量，单位：mm + k_e = 0.14 # 土壤蒸发系数 + ET_o = 8.188 # 参考作物蒸散量，单位：mm + E_i = k_e * ET_o # 第i天的蒸发量 + T_ew_i = 0 # 第i天湿润裸露表层土壤的蒸腾量，这里假定为0，单位：mm + f_w = 0.8 # 经灌溉湿润的土壤面积比值 + f_ew = 0.6 # 裸露湿润土壤面积比值 + D_ei = calculate_D_ei(D_e_prev, P_i, RO_i, I_i, f_w, f_ew) + print(f"第i天末土壤完全湿润后的累积蒸发（消耗）量为: {D_ei} mm") \ No newline at end of file diff --git a/irrgiation/utils.py b/irrgiation/utils.py new file mode 100644 index 0000000..54d3eb4 --- /dev/null +++ b/irrgiation/utils.py @@ -0,0 +1,10 @@ +dk_irriationInfo={ + "1号地":{ + "dkbm":1, + "irriationFile":"D:/悦柠/遥感/甘草/model_result/1号地_灌溉数据.xlsx" + }, + "5号地":{ + "dkbm":5, + "irriationFile":"D:/悦柠/遥感/甘草/model_result/5号地_灌溉数据.xlsx" + } +} \ No newline at end of file diff --git a/irrgiation/weatherAndSoilDataRequest.py b/irrgiation/weatherAndSoilDataRequest.py new file mode 100644 index 0000000..297ddbd --- /dev/null +++ b/irrgiation/weatherAndSoilDataRequest.py @@ -0,0 +1,336 @@ +from datetime import datetime, timezone, timedelta, date + +import pandas as pd +import requests +import json + +from matplotlib import pyplot as plt + + +#登录接口认证 + +def loginAuth(): + # 接口地址 + global auth_token + url = "https://xj-tb.maimaiag.com/api/auth/login" # 请替换为实际地址 + + # 请求头 + headers = { + "Content-Type": "application/json" + } + + # 请求体数据 + payload = { + "username": "read@xj.com", + "password": "maimai" + } + + try: + # 发送POST请求 + response = requests.post( + url, + headers=headers, + data=json.dumps(payload) # 或者直接使用 json=payload + ) + + # 检查响应状态码 + if response.status_code == 200: + print("登录成功！") + print("响应数据:", response.json()) + + # 通常登录接口会返回token，可以这样获取 + response_data = response.json() + auth_token = response_data.get("token") # 根据实际返回字段调整 + if auth_token: + print("获取到的认证Token:", auth_token) + else: + print(f"登录失败，状态码: {response.status_code}") + print("错误信息:", response.text) + + except requests.exceptions.RequestException as e: + print(f"请求发生异常: {str(e)}") + except json.JSONDecodeError as e: + print(f"JSON解析错误: {str(e)}") + return auth_token +# 调用数据接口--历史数据 +def getWeatherAndSoilData(auth_token,device_id,startTs,endTs,agg,interval,limit,orderBy,keys): + # 配置信息 + BASE_URL = "https://xj-tb.maimaiag.com/" # 请替换为实际API地址 + TOKEN = auth_token # 替换为实际的Bearer token + DEVICE_ID = device_id # 替换为实际的设备ID + # 构造请求URL + endpoint = f"/api/plugins/telemetry/DEVICE/{DEVICE_ID}/values/timeseries" + url = BASE_URL + endpoint + # 请求头 + headers = { + "Content-Type": "application/json", + "X-Authorization": f"Bearer+{TOKEN}" + } + # 查询参数 + params = { + "keys": keys, # 要查询的参数，多个用逗号分隔 + "startTs": startTs, # 开始时间戳（毫秒）- 可选 + "endTs": endTs, # 结束时间戳（毫秒）- 可选 + "agg": agg, # 聚合函数 - 可选 + "interval": interval, # 聚合间隔（毫秒）- 可选 + "limit": limit, # 返回条数 - 可选 + "orderBy": orderBy # 排序方式 - 可选 + } + try: + # 发送GET请求 + response = requests.get( + url, + headers=headers, + params=params + ) + # 检查响应状态码 + if response.status_code == 200: + # print("请求成功！") + data = response.json() + # print("响应数据:", json.dumps(data, indent=2)) + converted_data = {k: convert_timestamp_to_datetime(v) for k, v in data.items()} + return converted_data + elif response.status_code == 401: + print("认证失败，请检查Token是否正确") + elif response.status_code == 404: + print("设备不存在或路径错误") + else: + print(f"请求失败，状态码: {response.status_code}") + print("错误信息:", response.text) + + except requests.exceptions.RequestException as e: + print(f"请求发生异常: {str(e)}") + except json.JSONDecodeError as e: + print(f"JSON解析错误: {str(e)}") + print("原始响应:", response.text) +def convert_timestamp_to_datetime(data_dict): + """ + 递归地将字典中的所有'timestamp'字段转换为可读的日期时间格式 + 保留原始时间戳的同时添加新的'datetime'字段 + """ + if isinstance(data_dict, dict): + if "ts" in data_dict: + # 转换毫秒时间戳为datetime对象 + timestamp_sec = data_dict["ts"] / 1000 + dt = datetime.fromtimestamp(timestamp_sec) + # 添加新字段，保留原始时间戳 + data_dict["datetime"] = dt.strftime("%Y-%m-%d %H:%M:%S") + return data_dict + elif isinstance(data_dict, list): + return [convert_timestamp_to_datetime(item) for item in data_dict] + else: + return data_dict + +def dateToTimestamp(date): + # dt = datetime.strptime(date, "%Y-%m-%d %H:%M:%S") + if isinstance(date, datetime): + # 如果已经是datetime对象，直接使用 + dt = date + elif isinstance(date, str): + # 如果是字符串，解析为datetime对象 + dt = datetime.strptime(date, "%Y-%m-%d %H:%M:%S.%f") + else: + raise TypeError("date参数必须是字符串或datetime对象") + timestamp_ms = int(dt.timestamp() * 1000) + return timestamp_ms + +#调用接口获取未来气象数据 +def getFutureWeather(monitor_date): + base_url_future="https://api.open-meteo.com/" + end_point="v1/forecast" + current_date = datetime.now().date() + params = { + "latitude": "41.34252110189814", # 纬度 + "longitude": "86.28810755462379", # 经度 + "hourly": "temperature_2m,precipitation,precipitation_probability" + } + # if monitor_date.date()= 22] + plt.scatter(night_data['datetime'], night_data['value'], + color='red', + s=30, + label='22:00-06:00 Data') + + # 设置标题和坐标轴标签 + plt.title('Soil Moisture Time Series (2025-06-22 to 2025-06-30)', + fontsize=14, + pad=15) + plt.xlabel('Datetime', fontsize=12, labelpad=10) + plt.ylabel('Soil Moisture (%)', fontsize=12, labelpad=10) + + # 设置x轴刻度旋转，避免重叠 + plt.xticks(rotation=45, ha='right', fontsize=10) + + # 添加网格线 + plt.grid(axis='y', linestyle='--', alpha=0.7) + + # 添加图例 + plt.legend(loc='upper right', fontsize=10) + + # 调整布局 + plt.tight_layout() + + # 显示图表 + plt.show() +if __name__ == "__main__": + auth_token = loginAuth() + device_id="18d121f0-561b-11f0-a556-4f10f26fc07f" + agg="AVG" + interval=21600000 + limit="" + orderBy="" + keys="RAIN_FALL_REALTIME" + date_start = "2025-06-25 00:00:00.123" + date_end = "2025-07-01 23:59:59.123" + startTs=dateToTimestamp(date_start) + endTs=dateToTimestamp(date_end) + # data=get_soil_data_dk1(auth_token,startTs,endTs) + # print(data) + data=getWeatherAndSoilData(auth_token,device_id,startTs,endTs,"AVG",3600000,limit,orderBy,keys) + print(data['RAIN_FALL']) + # get_plot(data['1号传感器']['SOIL_MOISTURE_SURFACE']) + get_plot(data['RAIN_FALL_REALTIME']) \ No newline at end of file diff --git a/main.py b/main.py new file mode 100644 index 0000000..44eee24 --- /dev/null +++ b/main.py @@ -0,0 +1,42 @@ +from datetime import datetime, timedelta + +import uvicorn +from fastapi import FastAPI, Query +from pydantic import BaseModel + +from irrgiation.irrigateDecisionDemo import cal_irrigationByDay +# 创建 FastAPI 应用 +app = FastAPI() +# 定义请求模型 +class IrrigationRequest(BaseModel): + zwlx_name:str + dkbm: str + monitor_date: str + soil_key: str + weather_key: str + irrigation:float +@app.post("/calculate_irrigation/") +async def create_item(data: IrrigationRequest): + # 调用气象数据接口--获取具体日期的气象数据(气象) + try: + start_date = datetime.strptime(data.monitor_date, '%Y-%m-%d') + end_date =start_date+timedelta(days=1) + result = cal_irrigationByDay( + data.zwlx_name, + data.dkbm, + start_date, + end_date, + data.soil_key, + data.weather_key, + data.irrigation + ) + return {"status": "success", "data": result.to_dict(orient="records")} + except Exception as e: + return {"status": "error", "message": str(e)} + +@app.get("/items/lib") +async def get_lib(): + return {"message": "Library data"} + return result +if __name__ == "__main__": + uvicorn.run(app, host="127.0.0.1", port=8001) \ No newline at end of file diff --git a/pyeto/__init__.py b/pyeto/__init__.py new file mode 100644 index 0000000..d74a91e --- /dev/null +++ b/pyeto/__init__.py @@ -0,0 +1,99 @@ + +import pyeto + +from pyeto.convert import ( + celsius2kelvin, + kelvin2celsius, + deg2rad, + rad2deg, +) + +from pyeto.fao import ( + atm_pressure, + avp_from_tmin, + avp_from_rhmin_rhmax, + avp_from_rhmax, + avp_from_rhmean, + avp_from_tdew, + avp_from_twet_tdry, + cs_rad, + daily_mean_t, + daylight_hours, + delta_svp, + energy2evap, + et_rad, + fao56_penman_monteith, + hargreaves, + inv_rel_dist_earth_sun, + mean_svp, + monthly_soil_heat_flux, + monthly_soil_heat_flux2, + net_in_sol_rad, + net_out_lw_rad, + net_rad, + psy_const, + psy_const_of_psychrometer, + rh_from_avp_svp, + SOLAR_CONSTANT, + sol_dec, + sol_rad_from_sun_hours, + sol_rad_from_t, + sol_rad_island, + STEFAN_BOLTZMANN_CONSTANT, + sunset_hour_angle, + svp_from_t, + wind_speed_2m, +) + +from pyeto.thornthwaite import ( + thornthwaite, + monthly_mean_daylight_hours, +) + +__all__ = [ + # Unit conversions + 'celsius2kelvin', + 'deg2rad', + 'kelvin2celsius', + 'rad2deg', + + # FAO equations + 'atm_pressure', + 'avp_from_tmin', + 'avp_from_rhmin_rhmax', + 'avp_from_rhmax', + 'avp_from_rhmean', + 'avp_from_tdew', + 'avp_from_twet_tdry', + 'cs_rad', + 'daily_mean_t', + 'daylight_hours', + 'delta_svp', + 'energy2evap', + 'et_rad', + 'fao56_penman_monteith', + 'hargreaves', + 'inv_rel_dist_earth_sun', + 'mean_svp', + 'monthly_soil_heat_flux', + 'monthly_soil_heat_flux2', + 'net_in_sol_rad', + 'net_out_lw_rad', + 'net_rad', + 'psy_const', + 'psy_const_of_psychrometer', + 'rh_from_avp_svp', + 'SOLAR_CONSTANT', + 'sol_dec', + 'sol_rad_from_sun_hours', + 'sol_rad_from_t', + 'sol_rad_island', + 'STEFAN_BOLTZMANN_CONSTANT', + 'sunset_hour_angle', + 'svp_from_t', + 'wind_speed_2m', + + # Thornthwaite method + 'thornthwaite', + 'monthly_mean_daylight_hours', +] diff --git a/pyeto/__pycache__/__init__.cpython-310.pyc b/pyeto/__pycache__/__init__.cpython-310.pyc new file mode 100644 index 0000000..ee47222 Binary files /dev/null and b/pyeto/__pycache__/__init__.cpython-310.pyc differ diff --git a/pyeto/__pycache__/_check.cpython-310.pyc b/pyeto/__pycache__/_check.cpython-310.pyc new file mode 100644 index 0000000..592489e Binary files /dev/null and b/pyeto/__pycache__/_check.cpython-310.pyc differ diff --git a/pyeto/__pycache__/convert.cpython-310.pyc b/pyeto/__pycache__/convert.cpython-310.pyc new file mode 100644 index 0000000..36cfff5 Binary files /dev/null and b/pyeto/__pycache__/convert.cpython-310.pyc differ diff --git a/pyeto/__pycache__/fao.cpython-310.pyc b/pyeto/__pycache__/fao.cpython-310.pyc new file mode 100644 index 0000000..8a64cc0 Binary files /dev/null and b/pyeto/__pycache__/fao.cpython-310.pyc differ diff --git a/pyeto/__pycache__/thornthwaite.cpython-310.pyc b/pyeto/__pycache__/thornthwaite.cpython-310.pyc new file mode 100644 index 0000000..ab84565 Binary files /dev/null and b/pyeto/__pycache__/thornthwaite.cpython-310.pyc differ diff --git a/pyeto/_check.py b/pyeto/_check.py new file mode 100644 index 0000000..1c39fea --- /dev/null +++ b/pyeto/_check.py @@ -0,0 +1,69 @@ +""" +Internal validation functions. + +:copyright: (c) 2015 by Mark Richards. +:license: BSD 3-Clause, see LICENSE.txt for more details. +""" +from pyeto.convert import deg2rad + +# Internal constants +# Latitude +_MINLAT_RADIANS = deg2rad(-90.0) +_MAXLAT_RADIANS = deg2rad(90.0) + +# Solar declination +_MINSOLDEC_RADIANS = deg2rad(-23.5) +_MAXSOLDEC_RADIANS = deg2rad(23.5) + +# Sunset hour angle +_MINSHA_RADIANS = 0.0 +_MAXSHA_RADIANS = deg2rad(180) + + +def check_day_hours(hours, arg_name): + """ + Check that *hours* is in the range 1 to 24. + """ + if not 0 <= hours <= 24: + raise ValueError( + '{0} should be in range 0-24: {1!r}'.format(arg_name, hours)) + + +def check_doy(doy): + """ + Check day of the year is valid. + """ + if not 1 <= doy <= 366: + raise ValueError( + 'Day of the year (doy) must be in range 1-366: {0!r}'.format(doy)) + + +def check_latitude_rad(latitude): + if not _MINLAT_RADIANS <= latitude <= _MAXLAT_RADIANS: + raise ValueError( + 'latitude outside valid range {0!r} to {1!r} rad: {2!r}' + .format(_MINLAT_RADIANS, _MAXLAT_RADIANS, latitude)) + + +def check_sol_dec_rad(sd): + """ + Solar declination can vary between -23.5 and +23.5 degrees. + + See http://mypages.iit.edu/~maslanka/SolarGeo.pdf + """ + if not _MINSOLDEC_RADIANS <= sd <= _MAXSOLDEC_RADIANS: + raise ValueError( + 'solar declination outside valid range {0!r} to {1!r} rad: {2!r}' + .format(_MINSOLDEC_RADIANS, _MAXSOLDEC_RADIANS, sd)) + + +def check_sunset_hour_angle_rad(sha): + """ + Sunset hour angle has the range 0 to 180 degrees. + + See http://mypages.iit.edu/~maslanka/SolarGeo.pdf + """ + if not _MINSHA_RADIANS <= sha <= _MAXSHA_RADIANS: + raise ValueError( + 'sunset hour angle outside valid range {0!r} to {1!r} rad: {2!r}' + .format(_MINSHA_RADIANS, _MAXSHA_RADIANS, sha)) diff --git a/pyeto/convert.py b/pyeto/convert.py new file mode 100644 index 0000000..0ded485 --- /dev/null +++ b/pyeto/convert.py @@ -0,0 +1,52 @@ +""" +Unit conversion functions. + +:copyright: (c) 2015 by Mark Richards. +:license: BSD 3-Clause, see LICENSE.txt for more details. +""" + +import math + + +def celsius2kelvin(celsius): + """ + Convert temperature in degrees Celsius to degrees Kelvin. + + :param celsius: Degrees Celsius + :return: Degrees Kelvin + :rtype: float + """ + return celsius + 273.15 + + +def kelvin2celsius(kelvin): + """ + Convert temperature in degrees Kelvin to degrees Celsius. + + :param kelvin: Degrees Kelvin + :return: Degrees Celsius + :rtype: float + """ + return kelvin - 273.15 + + +def deg2rad(degrees): + """ + Convert angular degrees to radians + + :param degrees: Value in degrees to be converted. + :return: Value in radians + :rtype: float + """ + return degrees * (math.pi / 180.0) + + +def rad2deg(radians): + """ + Convert radians to angular degrees + + :param radians: Value in radians to be converted. + :return: Value in angular degrees + :rtype: float + """ + return radians * (180.0 / math.pi) diff --git a/pyeto/fao.py b/pyeto/fao.py new file mode 100644 index 0000000..48d0a51 --- /dev/null +++ b/pyeto/fao.py @@ -0,0 +1,737 @@ +""" +Library of functions for estimating reference evapotransporation (ETo) for +a grass reference crop using the FAO-56 Penman-Monteith and Hargreaves +equations. The library includes numerous functions for estimating missing +meteorological data. + +:copyright: (c) 2015 by Mark Richards. +:license: BSD 3-Clause, see LICENSE.txt for more details. +""" + +import math + +from ._check import ( + check_day_hours as _check_day_hours, + check_doy as _check_doy, + check_latitude_rad as _check_latitude_rad, + check_sol_dec_rad as _check_sol_dec_rad, + check_sunset_hour_angle_rad as _check_sunset_hour_angle_rad, +) + +#: Solar constant [ MJ m-2 min-1] +SOLAR_CONSTANT = 0.0820 + +# Stefan Boltzmann constant [MJ K-4 m-2 day-1] +STEFAN_BOLTZMANN_CONSTANT = 0.000000004903 # +"""Stefan Boltzmann constant [MJ K-4 m-2 day-1]""" + + +def atm_pressure(altitude): + """ + Estimate atmospheric pressure from altitude. + + Calculated using a simplification of the ideal gas law, assuming 20 degrees + Celsius for a standard atmosphere. Based on equation 7, page 62 in Allen + et al (1998). + + :param altitude: Elevation/altitude above sea level [m] + :return: atmospheric pressure [kPa] + :rtype: float + """ + tmp = (293.0 - (0.0065 * altitude)) / 293.0 + return math.pow(tmp, 5.26) * 101.3 + + +def avp_from_tmin(tmin): + """ + Estimate actual vapour pressure (*ea*) from minimum temperature. + + This method is to be used where humidity data are lacking or are of + questionable quality. The method assumes that the dewpoint temperature + is approximately equal to the minimum temperature (*tmin*), i.e. the + air is saturated with water vapour at *tmin*. + + **Note**: This assumption may not hold in arid/semi-arid areas. + In these areas it may be better to subtract 2 deg C from the + minimum temperature (see Annex 6 in FAO paper). + + Based on equation 48 in Allen et al (1998). + + :param tmin: Daily minimum temperature [deg C] + :return: Actual vapour pressure [kPa] + :rtype: float + """ + return 0.611 * math.exp((17.27 * tmin) / (tmin + 237.3)) + + +def avp_from_rhmin_rhmax(svp_tmin, svp_tmax, rh_min, rh_max): + """ + Estimate actual vapour pressure (*ea*) from saturation vapour pressure and + relative humidity. + + Based on FAO equation 17 in Allen et al (1998). + + :param svp_tmin: Saturation vapour pressure at daily minimum temperature + [kPa]. Can be estimated using ``svp_from_t()``. + :param svp_tmax: Saturation vapour pressure at daily maximum temperature + [kPa]. Can be estimated using ``svp_from_t()``. + :param rh_min: Minimum relative humidity [%] + :param rh_max: Maximum relative humidity [%] + :return: Actual vapour pressure [kPa] + :rtype: float + """ + tmp1 = svp_tmin * (rh_max / 100.0) + tmp2 = svp_tmax * (rh_min / 100.0) + return (tmp1 + tmp2) / 2.0 + + +def avp_from_rhmax(svp_tmin, rh_max): + """ + Estimate actual vapour pressure (*e*a) from saturation vapour pressure at + daily minimum temperature and maximum relative humidity + + Based on FAO equation 18 in Allen et al (1998). + + :param svp_tmin: Saturation vapour pressure at daily minimum temperature + [kPa]. Can be estimated using ``svp_from_t()``. + :param rh_max: Maximum relative humidity [%] + :return: Actual vapour pressure [kPa] + :rtype: float + """ + return svp_tmin * (rh_max / 100.0) + + +def avp_from_rhmean(svp_tmin, svp_tmax, rh_mean): + """ + Estimate actual vapour pressure (*ea*) from saturation vapour pressure at + daily minimum and maximum temperature, and mean relative humidity. + + Based on FAO equation 19 in Allen et al (1998). + + :param svp_tmin: Saturation vapour pressure at daily minimum temperature + [kPa]. Can be estimated using ``svp_from_t()``. + :param svp_tmax: Saturation vapour pressure at daily maximum temperature + [kPa]. Can be estimated using ``svp_from_t()``. + :param rh_mean: Mean relative humidity [%] (average of RH min and RH max). + :return: Actual vapour pressure [kPa] + :rtype: float + """ + return (rh_mean / 100.0) * ((svp_tmax + svp_tmin) / 2.0) + + +def avp_from_tdew(tdew): + """ + Estimate actual vapour pressure (*ea*) from dewpoint temperature. + + Based on equation 14 in Allen et al (1998). As the dewpoint temperature is + the temperature to which air needs to be cooled to make it saturated, the + actual vapour pressure is the saturation vapour pressure at the dewpoint + temperature. + + This method is preferable to calculating vapour pressure from + minimum temperature. + + :param tdew: Dewpoint temperature [deg C] + :return: Actual vapour pressure [kPa] + :rtype: float + """ + return 0.6108 * math.exp((17.27 * tdew) / (tdew + 237.3)) + + +def avp_from_twet_tdry(twet, tdry, svp_twet, psy_const): + """ + Estimate actual vapour pressure (*ea*) from wet and dry bulb temperature. + + Based on equation 15 in Allen et al (1998). As the dewpoint temperature + is the temperature to which air needs to be cooled to make it saturated, the + actual vapour pressure is the saturation vapour pressure at the dewpoint + temperature. + + This method is preferable to calculating vapour pressure from + minimum temperature. + + Values for the psychrometric constant of the psychrometer (*psy_const*) + can be calculated using ``psyc_const_of_psychrometer()``. + + :param twet: Wet bulb temperature [deg C] + :param tdry: Dry bulb temperature [deg C] + :param svp_twet: Saturated vapour pressure at the wet bulb temperature + [kPa]. Can be estimated using ``svp_from_t()``. + :param psy_const: Psychrometric constant of the pyschrometer [kPa deg C-1]. + Can be estimated using ``psy_const()`` or + ``psy_const_of_psychrometer()``. + :return: Actual vapour pressure [kPa] + :rtype: float + """ + return svp_twet - (psy_const * (tdry - twet)) + + +def cs_rad(altitude, et_rad): + """ + Estimate clear sky radiation from altitude and extraterrestrial radiation. + + Based on equation 37 in Allen et al (1998) which is recommended when + calibrated Angstrom values are not available. + + :param altitude: Elevation above sea level [m] + :param et_rad: Extraterrestrial radiation [MJ m-2 day-1]. Can be + estimated using ``et_rad()``. + :return: Clear sky radiation [MJ m-2 day-1] + :rtype: float + """ + return (0.00002 * altitude + 0.75) * et_rad + + +def daily_mean_t(tmin, tmax): + """ + Estimate mean daily temperature from the daily minimum and maximum + temperatures. + + :param tmin: Minimum daily temperature [deg C] + :param tmax: Maximum daily temperature [deg C] + :return: Mean daily temperature [deg C] + :rtype: float + """ + return (tmax + tmin) / 2.0 + + +def daylight_hours(sha): + """ + Calculate daylight hours from sunset hour angle. + + Based on FAO equation 34 in Allen et al (1998). + + :param sha: Sunset hour angle [rad]. Can be calculated using + ``sunset_hour_angle()``. + :return: Daylight hours. + :rtype: float + """ + _check_sunset_hour_angle_rad(sha) + return (24.0 / math.pi) * sha + + +def delta_svp(t): + """ + Estimate the slope of the saturation vapour pressure curve at a given + temperature. + + Based on equation 13 in Allen et al (1998). If using in the Penman-Monteith + *t* should be the mean air temperature. + + :param t: Air temperature [deg C]. Use mean air temperature for use in + Penman-Monteith. + :return: Saturation vapour pressure [kPa degC-1] + :rtype: float + """ + tmp = 4098 * (0.6108 * math.exp((17.27 * t) / (t + 237.3))) + return tmp / math.pow((t + 237.3), 2) + + +def energy2evap(energy): + """ + Convert energy (e.g. radiation energy) in MJ m-2 day-1 to the equivalent + evaporation, assuming a grass reference crop. + + Energy is converted to equivalent evaporation using a conversion + factor equal to the inverse of the latent heat of vapourisation + (1 / lambda = 0.408). + + Based on FAO equation 20 in Allen et al (1998). + + :param energy: Energy e.g. radiation or heat flux [MJ m-2 day-1]. + :return: Equivalent evaporation [mm day-1]. + :rtype: float + """ + return 0.408 * energy + + +def et_rad(latitude, sol_dec, sha, ird): + """ + Estimate daily extraterrestrial radiation (*Ra*, 'top of the atmosphere + radiation'). + + Based on equation 21 in Allen et al (1998). If monthly mean radiation is + required make sure *sol_dec*. *sha* and *irl* have been calculated using + the day of the year that corresponds to the middle of the month. + + **Note**: From Allen et al (1998): "For the winter months in latitudes + greater than 55 degrees (N or S), the equations have limited validity. + Reference should be made to the Smithsonian Tables to assess possible + deviations." + + :param latitude: Latitude [radians] + :param sol_dec: Solar declination [radians]. Can be calculated using + ``sol_dec()``. + :param sha: Sunset hour angle [radians]. Can be calculated using + ``sunset_hour_angle()``. + :param ird: Inverse relative distance earth-sun [dimensionless]. Can be + calculated using ``inv_rel_dist_earth_sun()``. + :return: Daily extraterrestrial radiation [MJ m-2 day-1] + :rtype: float + """ + _check_latitude_rad(latitude) + _check_sol_dec_rad(sol_dec) + _check_sunset_hour_angle_rad(sha) + + tmp1 = (24.0 * 60.0) / math.pi + tmp2 = sha * math.sin(latitude) * math.sin(sol_dec) + tmp3 = math.cos(latitude) * math.cos(sol_dec) * math.sin(sha) + return tmp1 * SOLAR_CONSTANT * ird * (tmp2 + tmp3) + + +def fao56_penman_monteith(net_rad, t, ws, svp, avp, delta_svp, psy, shf=0.0): + """ + Estimate reference evapotranspiration (ETo) from a hypothetical + short grass reference surface using the FAO-56 Penman-Monteith equation. + + Based on equation 6 in Allen et al (1998). + + :param net_rad: Net radiation at crop surface [MJ m-2 day-1]. If + necessary this can be estimated using ``net_rad()``. + :param t: Air temperature at 2 m height [deg Kelvin]. + :param ws: Wind speed at 2 m height [m s-1]. If not measured at 2m, + convert using ``wind_speed_at_2m()``. + :param svp: Saturation vapour pressure [kPa]. Can be estimated using + ``svp_from_t()''. + :param avp: Actual vapour pressure [kPa]. Can be estimated using a range + of functions with names beginning with 'avp_from'. + :param delta_svp: Slope of saturation vapour pressure curve [kPa degC-1]. + Can be estimated using ``delta_svp()``. + :param psy: Psychrometric constant [kPa deg C]. Can be estimatred using + ``psy_const_of_psychrometer()`` or ``psy_const()``. + :param shf: Soil heat flux (G) [MJ m-2 day-1] (default is 0.0, which is + reasonable for a daily or 10-day time steps). For monthly time steps + *shf* can be estimated using ``monthly_soil_heat_flux()`` or + ``monthly_soil_heat_flux2()``. + :return: Reference evapotranspiration (ETo) from a hypothetical + grass reference surface [mm day-1]. + :rtype: float + """ + a1 = (0.408 * (net_rad - shf) * delta_svp / + (delta_svp + (psy * (1 + 0.34 * ws)))) + a2 = (900 * ws / (t+273) * (svp - avp) * psy / + (delta_svp + (psy * (1 + 0.34 * ws)))) + return a1 + a2 + + +def hargreaves(tmin, tmax, tmean, et_rad): + """ + Estimate reference evapotranspiration over grass (ETo) using the Hargreaves + equation. + + Generally, when solar radiation data, relative humidity data + and/or wind speed data are missing, it is better to estimate them using + the functions available in this module, and then calculate ETo + the FAO Penman-Monteith equation. However, as an alternative, ETo can be + estimated using the Hargreaves ETo equation. + + Based on equation 52 in Allen et al (1998). + + :param tmin: Minimum daily temperature [deg C] + :param tmax: Maximum daily temperature [deg C] + :param tmean: Mean daily temperature [deg C]. If emasurements not + available it can be estimated as (*tmin* + *tmax*) / 2. + :param et_rad: Extraterrestrial radiation (Ra) [MJ m-2 day-1]. Can be + estimated using ``et_rad()``. + :return: Reference evapotranspiration over grass (ETo) [mm day-1] + :rtype: float + """ + # Note, multiplied by 0.408 to convert extraterrestrial radiation could + # be given in MJ m-2 day-1 rather than as equivalent evaporation in + # mm day-1 + return 0.0023 * (tmean + 17.8) * (tmax - tmin) ** 0.5 * 0.408 * et_rad + + +def inv_rel_dist_earth_sun(day_of_year): + """ + Calculate the inverse relative distance between earth and sun from + day of the year. + + Based on FAO equation 23 in Allen et al (1998). + + :param day_of_year: Day of the year [1 to 366] + :return: Inverse relative distance between earth and the sun + :rtype: float + """ + _check_doy(day_of_year) + return 1 + (0.033 * math.cos((2.0 * math.pi / 365.0) * day_of_year)) + + +def mean_svp(tmin, tmax): + """ + Estimate mean saturation vapour pressure, *es* [kPa] from minimum and + maximum temperature. + + Based on equations 11 and 12 in Allen et al (1998). + + Mean saturation vapour pressure is calculated as the mean of the + saturation vapour pressure at tmax (maximum temperature) and tmin + (minimum temperature). + + :param tmin: Minimum temperature [deg C] + :param tmax: Maximum temperature [deg C] + :return: Mean saturation vapour pressure (*es*) [kPa] + :rtype: float + """ + return (svp_from_t(tmin) + svp_from_t(tmax)) / 2.0 + + +def monthly_soil_heat_flux(t_month_prev, t_month_next): + """ + Estimate monthly soil heat flux (Gmonth) from the mean air temperature of + the previous and next month, assuming a grass crop. + + Based on equation 43 in Allen et al (1998). If the air temperature of the + next month is not known use ``monthly_soil_heat_flux2()`` instead. The + resulting heat flux can be converted to equivalent evaporation [mm day-1] + using ``energy2evap()``. + + :param t_month_prev: Mean air temperature of the previous month + [deg Celsius] + :param t_month2_next: Mean air temperature of the next month [deg Celsius] + :return: Monthly soil heat flux (Gmonth) [MJ m-2 day-1] + :rtype: float + """ + return 0.07 * (t_month_next - t_month_prev) + + +def monthly_soil_heat_flux2(t_month_prev, t_month_cur): + """ + Estimate monthly soil heat flux (Gmonth) [MJ m-2 day-1] from the mean + air temperature of the previous and current month, assuming a grass crop. + + Based on equation 44 in Allen et al (1998). If the air temperature of the + next month is available, use ``monthly_soil_heat_flux()`` instead. The + resulting heat flux can be converted to equivalent evaporation [mm day-1] + using ``energy2evap()``. + + Arguments: + :param t_month_prev: Mean air temperature of the previous month + [deg Celsius] + :param t_month_cur: Mean air temperature of the current month [deg Celsius] + :return: Monthly soil heat flux (Gmonth) [MJ m-2 day-1] + :rtype: float + """ + return 0.14 * (t_month_cur - t_month_prev) + + +def net_in_sol_rad(sol_rad, albedo=0.23): + """ + Calculate net incoming solar (or shortwave) radiation from gross + incoming solar radiation, assuming a grass reference crop. + + Net incoming solar radiation is the net shortwave radiation resulting + from the balance between incoming and reflected solar radiation. The + output can be converted to equivalent evaporation [mm day-1] using + ``energy2evap()``. + + Based on FAO equation 38 in Allen et al (1998). + + :param sol_rad: Gross incoming solar radiation [MJ m-2 day-1]. If + necessary this can be estimated using functions whose name + begins with 'sol_rad_from'. + :param albedo: Albedo of the crop as the proportion of gross incoming solar + radiation that is reflected by the surface. Default value is 0.23, + which is the value used by the FAO for a short grass reference crop. + Albedo can be as high as 0.95 for freshly fallen snow and as low as + 0.05 for wet bare soil. A green vegetation over has an albedo of + about 0.20-0.25 (Allen et al, 1998). + :return: Net incoming solar (or shortwave) radiation [MJ m-2 day-1]. + :rtype: float + """ + return (1 - albedo) * sol_rad + + +def net_out_lw_rad(tmin, tmax, sol_rad, cs_rad, avp): + """ + Estimate net outgoing longwave radiation. + + This is the net longwave energy (net energy flux) leaving the + earth's surface. It is proportional to the absolute temperature of + the surface raised to the fourth power according to the Stefan-Boltzmann + law. However, water vapour, clouds, carbon dioxide and dust are absorbers + and emitters of longwave radiation. This function corrects the Stefan- + Boltzmann law for humidity (using actual vapor pressure) and cloudiness + (using solar radiation and clear sky radiation). The concentrations of all + other absorbers are assumed to be constant. + + The output can be converted to equivalent evaporation [mm day-1] using + ``energy2evap()``. + + Based on FAO equation 39 in Allen et al (1998). + + :param tmin: Absolute daily minimum temperature [degrees Kelvin] + :param tmax: Absolute daily maximum temperature [degrees Kelvin] + :param sol_rad: Solar radiation [MJ m-2 day-1]. If necessary this can be + estimated using ``sol+rad()``. + :param cs_rad: Clear sky radiation [MJ m-2 day-1]. Can be estimated using + ``cs_rad()``. + :param avp: Actual vapour pressure [kPa]. Can be estimated using functions + with names beginning with 'avp_from'. + :return: Net outgoing longwave radiation [MJ m-2 day-1] + :rtype: float + """ + tmax_k = tmax + 273.15 + tmin_k = tmin + 273.15 + tmp1 = (STEFAN_BOLTZMANN_CONSTANT * + ((math.pow(tmax_k, 4) + math.pow(tmin_k, 4)) / 2)) + tmp2 = (0.34 - (0.14 * math.sqrt(avp))) + tmp3 = 1.35 * (sol_rad / cs_rad) - 0.35 + return tmp1 * tmp2 * tmp3 + + +def net_rad(ni_sw_rad, no_lw_rad): + """ + Calculate daily net radiation at the crop surface, assuming a grass + reference crop. + + Net radiation is the difference between the incoming net shortwave (or + solar) radiation and the outgoing net longwave radiation. Output can be + converted to equivalent evaporation [mm day-1] using ``energy2evap()``. + + Based on equation 40 in Allen et al (1998). + + :param ni_sw_rad: Net incoming shortwave radiation [MJ m-2 day-1]. Can be + estimated using ``net_in_sol_rad()``. + :param no_lw_rad: Net outgoing longwave radiation [MJ m-2 day-1]. Can be + estimated using ``net_out_lw_rad()``. + :return: Daily net radiation [MJ m-2 day-1]. + :rtype: float + """ + return ni_sw_rad - no_lw_rad + + +def psy_const(atmos_pres): + """ + Calculate the psychrometric constant. + + This method assumes that the air is saturated with water vapour at the + minimum daily temperature. This assumption may not hold in arid areas. + + Based on equation 8, page 95 in Allen et al (1998). + + :param atmos_pres: Atmospheric pressure [kPa]. Can be estimated using + ``atm_pressure()``. + :return: Psychrometric constant [kPa degC-1]. + :rtype: float + """ + return 0.000665 * atmos_pres + + +def psy_const_of_psychrometer(psychrometer, atmos_pres): + """ + Calculate the psychrometric constant for different types of + psychrometer at a given atmospheric pressure. + + Based on FAO equation 16 in Allen et al (1998). + + :param psychrometer: Integer between 1 and 3 which denotes type of + psychrometer: + 1. ventilated (Asmann or aspirated type) psychrometer with + an air movement of approximately 5 m/s + 2. natural ventilated psychrometer with an air movement + of approximately 1 m/s + 3. non ventilated psychrometer installed indoors + :param atmos_pres: Atmospheric pressure [kPa]. Can be estimated using + ``atm_pressure()``. + :return: Psychrometric constant [kPa degC-1]. + :rtype: float + """ + # Select coefficient based on type of ventilation of the wet bulb + if psychrometer == 1: + psy_coeff = 0.000662 + elif psychrometer == 2: + psy_coeff = 0.000800 + elif psychrometer == 3: + psy_coeff = 0.001200 + else: + raise ValueError( + 'psychrometer should be in range 1 to 3: {0!r}'.format(psychrometer)) + + return psy_coeff * atmos_pres + + +def rh_from_avp_svp(avp, svp): + """ + Calculate relative humidity as the ratio of actual vapour pressure + to saturation vapour pressure at the same temperature. + + See Allen et al (1998), page 67 for details. + + :param avp: Actual vapour pressure [units do not matter so long as they + are the same as for *svp*]. Can be estimated using functions whose + name begins with 'avp_from'. + :param svp: Saturated vapour pressure [units do not matter so long as they + are the same as for *avp*]. Can be estimated using ``svp_from_t()``. + :return: Relative humidity [%]. + :rtype: float + """ + return 100.0 * avp / svp + + +def sol_dec(day_of_year): + """ + Calculate solar declination from day of the year. + + Based on FAO equation 24 in Allen et al (1998). + + :param day_of_year: Day of year integer between 1 and 365 or 366). + :return: solar declination [radians] + :rtype: float + """ + _check_doy(day_of_year) + return 0.409 * math.sin(((2.0 * math.pi / 365.0) * day_of_year - 1.39)) + + +def sol_rad_from_sun_hours(daylight_hours, sunshine_hours, et_rad): + """ + Calculate incoming solar (or shortwave) radiation, *Rs* (radiation hitting + a horizontal plane after scattering by the atmosphere) from relative + sunshine duration. + + If measured radiation data are not available this method is preferable + to calculating solar radiation from temperature. If a monthly mean is + required then divide the monthly number of sunshine hours by number of + days in the month and ensure that *et_rad* and *daylight_hours* was + calculated using the day of the year that corresponds to the middle of + the month. + + Based on equations 34 and 35 in Allen et al (1998). + + :param dl_hours: Number of daylight hours [hours]. Can be calculated + using ``daylight_hours()``. + :param sunshine_hours: Sunshine duration [hours]. + :param et_rad: Extraterrestrial radiation [MJ m-2 day-1]. Can be + estimated using ``et_rad()``. + :return: Incoming solar (or shortwave) radiation [MJ m-2 day-1] + :rtype: float + """ + _check_day_hours(sunshine_hours, 'sun_hours') + _check_day_hours(daylight_hours, 'daylight_hours') + + # 0.5 and 0.25 are default values of regression constants (Angstrom values) + # recommended by FAO when calibrated values are unavailable. + return (0.5 * sunshine_hours / daylight_hours + 0.25) * et_rad + + +def sol_rad_from_t(et_rad, cs_rad, tmin, tmax, coastal): + """ + Estimate incoming solar (or shortwave) radiation, *Rs*, (radiation hitting + a horizontal plane after scattering by the atmosphere) from min and max + temperature together with an empirical adjustment coefficient for + 'interior' and 'coastal' regions. + + The formula is based on equation 50 in Allen et al (1998) which is the + Hargreaves radiation formula (Hargreaves and Samani, 1982, 1985). This + method should be used only when solar radiation or sunshine hours data are + not available. It is only recommended for locations where it is not + possible to use radiation data from a regional station (either because + climate conditions are heterogeneous or data are lacking). + + **NOTE**: this method is not suitable for island locations due to the + moderating effects of the surrounding water. + + :param et_rad: Extraterrestrial radiation [MJ m-2 day-1]. Can be + estimated using ``et_rad()``. + :param cs_rad: Clear sky radiation [MJ m-2 day-1]. Can be estimated + using ``cs_rad()``. + :param tmin: Daily minimum temperature [deg C]. + :param tmax: Daily maximum temperature [deg C]. + :param coastal: ``True`` if site is a coastal location, situated on or + adjacent to coast of a large land mass and where air masses are + influenced by a nearby water body, ``False`` if interior location + where land mass dominates and air masses are not strongly influenced + by a large water body. + :return: Incoming solar (or shortwave) radiation (Rs) [MJ m-2 day-1]. + :rtype: float + """ + # Determine value of adjustment coefficient [deg C-0.5] for + # coastal/interior locations + if coastal: + adj = 0.19 + else: + adj = 0.16 + + sol_rad = adj * math.sqrt(tmax - tmin) * et_rad + + # The solar radiation value is constrained by the clear sky radiation + return min(sol_rad, cs_rad) + + +def sol_rad_island(et_rad): + """ + Estimate incoming solar (or shortwave) radiation, *Rs* (radiation hitting + a horizontal plane after scattering by the atmosphere) for an island + location. + + An island is defined as a land mass with width perpendicular to the + coastline <= 20 km. Use this method only if radiation data from + elsewhere on the island is not available. + + **NOTE**: This method is only applicable for low altitudes (0-100 m) + and monthly calculations. + + Based on FAO equation 51 in Allen et al (1998). + + :param et_rad: Extraterrestrial radiation [MJ m-2 day-1]. Can be + estimated using ``et_rad()``. + :return: Incoming solar (or shortwave) radiation [MJ m-2 day-1]. + :rtype: float + """ + return (0.7 * et_rad) - 4.0 + + +def sunset_hour_angle(latitude, sol_dec): + """ + Calculate sunset hour angle (*Ws*) from latitude and solar + declination. + + Based on FAO equation 25 in Allen et al (1998). + + :param latitude: Latitude [radians]. Note: *latitude* should be negative + if it in the southern hemisphere, positive if in the northern + hemisphere. + :param sol_dec: Solar declination [radians]. Can be calculated using + ``sol_dec()``. + :return: Sunset hour angle [radians]. + :rtype: float + """ + _check_latitude_rad(latitude) + _check_sol_dec_rad(sol_dec) + + cos_sha = -math.tan(latitude) * math.tan(sol_dec) + # If tmp is >= 1 there is no sunset, i.e. 24 hours of daylight + # If tmp is <= 1 there is no sunrise, i.e. 24 hours of darkness + # See http://www.itacanet.org/the-sun-as-a-source-of-energy/ + # part-3-calculating-solar-angles/ + # Domain of acos is -1 <= x <= 1 radians (this is not mentioned in FAO-56!) + return math.acos(min(max(cos_sha, -1.0), 1.0)) + + +def svp_from_t(t): + """ + Estimate saturation vapour pressure (*es*) from air temperature. + + Based on equations 11 and 12 in Allen et al (1998). + + :param t: Temperature [deg C] + :return: Saturation vapour pressure [kPa] + :rtype: float + """ + return 0.6108 * math.exp((17.27 * t) / (t + 237.3)) + + +def wind_speed_2m(ws, z): + """ + Convert wind speed measured at different heights above the soil + surface to wind speed at 2 m above the surface, assuming a short grass + surface. + + Based on FAO equation 47 in Allen et al (1998). + + :param ws: Measured wind speed [m s-1] + :param z: Height of wind measurement above ground surface [m] + :return: Wind speed at 2 m above the surface [m s-1] + :rtype: float + """ + return ws * (4.87 / math.log((67.8 * z) - 5.42)) diff --git a/pyeto/thornthwaite.py b/pyeto/thornthwaite.py new file mode 100644 index 0000000..18e7f7c --- /dev/null +++ b/pyeto/thornthwaite.py @@ -0,0 +1,119 @@ +""" +Calculate potential evapotranspiration using the Thornthwaite (1948 method) + +:copyright: (c) 2015 by Mark Richards. +:license: BSD 3-Clause, see LICENSE.txt for more details. + +References +---------- +Thornthwaite CW (1948) An approach toward a rational classification of + climate. Geographical Review, 38, 55-94. +""" + +import calendar + +from . import fao +from ._check import check_latitude_rad as _check_latitude_rad + +_MONTHDAYS = (31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31) +_LEAP_MONTHDAYS = (31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31) + + +def thornthwaite(monthly_t, monthly_mean_dlh, year=None): + """ + Estimate monthly potential evapotranspiration (PET) using the + Thornthwaite (1948) method. + + Thornthwaite equation: + + *PET* = 1.6 (*L*/12) (*N*/30) (10*Ta* / *I*)***a* + + where: + + * *Ta* is the mean daily air temperature [deg C, if negative use 0] of the + month being calculated + * *N* is the number of days in the month being calculated + * *L* is the mean day length [hours] of the month being calculated + * *a* = (6.75 x 10-7)*I***3 - (7.71 x 10-5)*I***2 + (1.792 x 10-2)*I* + 0.49239 + * *I* is a heat index which depends on the 12 monthly mean temperatures and + is calculated as the sum of (*Tai* / 5)**1.514 for each month, where + Tai is the air temperature for each month in the year + + :param monthly_t: Iterable containing mean daily air temperature for each + month of the year [deg C]. + :param monthly_mean_dlh: Iterable containing mean daily daylight + hours for each month of the year (hours]. These can be calculated + using ``monthly_mean_daylight_hours()``. + :param year: Year for which PET is required. The only effect of year is + to change the number of days in February to 29 if it is a leap year. + If it is left as the default (None), then the year is assumed not to + be a leap year. + :return: Estimated monthly potential evaporation of each month of the year + [mm/month] + :rtype: List of floats + """ + if len(monthly_t) != 12: + raise ValueError( + 'monthly_t should be length 12 but is length {0}.' + .format(len(monthly_t))) + if len(monthly_mean_dlh) != 12: + raise ValueError( + 'monthly_mean_dlh should be length 12 but is length {0}.' + .format(len(monthly_mean_dlh))) + + if year is None or not calendar.isleap(year): + month_days = _MONTHDAYS + else: + month_days = _LEAP_MONTHDAYS + + # Negative temperatures should be set to zero + adj_monthly_t = [t * (t >= 0) for t in monthly_t] + + # Calculate the heat index (I) + I = 0.0 + for Tai in adj_monthly_t: + if Tai / 5.0 > 0.0: + I += (Tai / 5.0) ** 1.514 + + a = (6.75e-07 * I ** 3) - (7.71e-05 * I ** 2) + (1.792e-02 * I) + 0.49239 + + pet = [] + for Ta, L, N in zip(adj_monthly_t, monthly_mean_dlh, month_days): + # Multiply by 10 to convert cm/month --> mm/month + pet.append( + 1.6 * (L / 12.0) * (N / 30.0) * ((10.0 * Ta / I) ** a) * 10.0) + + return pet + + +def monthly_mean_daylight_hours(latitude, year=None): + """ + Calculate mean daylight hours for each month of the year for a given + latitude. + + :param latitude: Latitude [radians] + :param year: Year for the daylight hours are required. The only effect of + *year* is to change the number of days in Feb to 29 if it is a leap + year. If left as the default, None, then a normal (non-leap) year is + assumed. + :return: Mean daily daylight hours of each month of a year [hours] + :rtype: List of floats. + """ + _check_latitude_rad(latitude) + + if year is None or not calendar.isleap(year): + month_days = _MONTHDAYS + else: + month_days = _LEAP_MONTHDAYS + monthly_mean_dlh = [] + doy = 1 # Day of the year + for mdays in month_days: + dlh = 0.0 # Cumulative daylight hours for the month + for daynum in range(1, mdays + 1): + sd = fao.sol_dec(doy) + sha = fao.sunset_hour_angle(latitude, sd) + dlh += fao.daylight_hours(sha) + doy += 1 + # Calc mean daylight hours of the month + monthly_mean_dlh.append(dlh / mdays) + return monthly_mean_dlh diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..f541c9e Binary files /dev/null and b/requirements.txt differ