# python_nc_code

**Repository Path**: xhfan/python_nc_code

## Basic Information

- **Project Name**: python_nc_code
- **Description**: 利用python+cdo处理nc文件
- **Primary Language**: Unknown
- **License**: Not specified
- **Default Branch**: master
- **Homepage**: None
- **GVP Project**: No

## Statistics

- **Stars**: 1
- **Forks**: 0
- **Created**: 2023-12-24
- **Last Updated**: 2025-02-20

## Categories & Tags

**Categories**: Uncategorized

**Tags**: None

## README

### 1. 概述

- 根据气温、降水、风速等nc格式数据识别极端天气气候事件。
- 数据为1981-2020年的nc格式的数据，分别以198101.nc，198102.nc，……202012.nc存储在
  2m_temperature,2m_dewpoint_temperature, total_precipitation, 10m_u_component_of_wind, 10m_v_component_of_wind
- 极端天气气候事件包括 **[p1:种植-拔节; p2:拔节-抽穗; p3:抽穗-成熟]**：

  + (a) 小于0℃温度之和、小时数之和 **[frez_t_p1, frez_t_p2]**
  + (b) 大于0℃温度之和、小时数之和 **[GDD_total_p1,GDD_total_p2,GDD_total_p3]**。详细地：
    - (b1) 0-5℃温度之和、小时数之和 **[GDDlow_t_p1]**
    - (b2) 5-10℃温度之和、小时数之和 **[GDDmed_t_p1]**
    - (b3) 大于10℃温度之和、小时数之和 **[GDDhigh_t_p1]**
    - (b4) 0-10℃温度之和、小时数之和 **[GDDlow_t_p2]**
    - (b5) 10-20℃温度之和、小时数之和 **[GDDmed_t_p2]**
    - (b6) 大于20℃温度之和、小时数之和 **[GDDhigh_t_p2]**
    - (b7) 0-20℃温度之和、小时数之和 **[GDDlow_t_p3]**
    - (b8) 20-30℃温度之和、小时数之和 **[GDDmed_t_p3]**
    - (b9) 大于30℃温度之和、小时数之和 **[GDDhigh_t_p3]**
  + (c) 高于30℃温度之和 **[EDD_t_p2,EDD_t_p2]**
    - (c1) 30-31℃温度之和、小时数之和 **[t30_31_p2]**
    - (c2) 31-32℃温度之和、小时数之和 **[t31_32_p2]**
    - (c3) 32-33℃温度之和、小时数之和 **[t32_33_p2]**
    - (c4) 33-35℃温度之和、小时数之和 **[t33_35_p2]**
    - (c5) 35-36℃温度之和、小时数之和 **[t35_36_p2]**

  - (d) 降雨量之和 **[tp_total_p1,tp_total_p2,tp_total_p3]**
    - (d1) 0-2cm降雨量之和、小时数之和 **[tp_2_p1]**
    - (d2) 2-4cm降雨量之和、小时数之和 **[tp2_4_p1]**
    - (d3) 4-8cm降雨量之和、小时数之和 **[tp4_8_p1]**
    - (d4) 8-20cm降雨量之和、小时数之和 **[tp8_20_p1]**
    - (d5) 20-50cm降雨量之和、小时数之和 **[tp20_50_p1]**
    - (d6) 高于50cm降雨量之和、小时数之和 **[tp50_p1]**
- 复合极端天气气候事件包括：

  - 干热风日数和小时数 **[HDW_h,HDW_d]**
    - (a) 温度T(31-32℃)&相对湿度RH≤30%&风速U≥3 的小时数、日数 **[hdw1_h,hdw1_d]**
    - (b) 温度T(32-33℃)&相对湿度RH≤30%&风速U≥3 的小时数、日数 **[hdw2_h,hdw2_d]**
    - (c) 温度T(33-35℃)&相对湿度RH≤25%&风速U≥3 的小时数、日数 **[hdw3_h,hdw3_d]**
    - (d) 温度T(35-36℃)&相对湿度RH≤25%&风速U≥3 的小时数、日数 **[hdw4_h,hdw4_d]**
    - (e) 温度T(＞36℃)&相对湿度RH≤25%&风速U≥3 的小时数、日数 **[hdw5_h,hdw5_d]**

- <font color='red'><b>计算极端天气气候事件的方法主要有Python+ArcPy方法和Python+cdo方法。~~需要注意的是在计算极端事件时，需要首先生成StartEndDates_lat_***_lon_***文件。~~</b></font>

### 2. python方法

> 主要思路：从nc文件中按小时提取每个波段，然后按条件对单个波段按要求处理之后合并为日数据数据文件，接着利用待提取的经纬度列表对日数据进行
> 数据的提取，结果存入CSV文件中以备进一步处理。<br>另：文件夹命名规则建议：**文件夹名称“###_hourly_bands_tif“表示从nc文件中按要求的不同的取值范围（bands）导出单日的每小时（hourly）的数据（tif格式）；
> 而文件夹名称“###_hourly_count_bands_tif“表示从nc文件中按要求统计在单日的每小时（hourly）的tif格式的数据中，识别并标注每个不同的取值范围（bands）的
> 数据值是否存在（用0或1标注），从而在“###_dailly_count_bands_tif”文件夹中存储每日的小时数，**

详细步骤如下：

- 极端天气气候事件计算：
  - （1）通过【01_01_temperature_houly_to_GDD_bands.py】文件，从2m_temperature读取温度的小时数据，生成与FDD和GDD相关的bands小时数据（或统计是否在此范围中的1、0数据），结果分别存放在temperature_hourly_band_tif（或temperature_hourly_count_band_tif）文件夹下的FDD_0，GDD_0等子文件夹中（文件夹根据具体的band自行定义）。如：temperature_hourly_band_tif/GDD_step3_20_30表示的是统计每个小时温度介于20度至30度之间的格点。 **若是，则保留原值，若非，则为空：结果保存在 “###_hourly_bands_tif”中；若是，则记为1，否则记为空，则结果保存在"###_hourly_count_bands_tif"中。** 【上述极端天气和气候事件(a)和(b)】
  - （2）通过【01_02_temperature_houly_to_temperature_bands.py】文件，从2m_temperature读取温度的小时数据，生成不同bands的温度之和。结果保存在temperature_hourly_band_tif（或temperature_hourly_count_band_tif）文件夹下的自定义子文件夹中。如：temperature_hourly_band_tif/t30_#表示的是统计每个小时温度大于30度的格点。 **若是，则保留原值，若非，则为空：结果保存在 “###_hourly_bands_tif”中；若是，则记为1，否则记为空，则结果保存在"###_hourly_count_bands_tif"中。**【上述极端天气和气候事件(c)】
  - （3）通过【03_01_prcp_houly_to_bands.py】文件，从total_precipitation文件夹中读取降水每小时数据，生成不同bands的降水之和。文件夹命名规则如（1）和（2）。【上述极端天气和气候事件(d)】
  - （4）通过【02_hourlyTOdaily.py】对（1）~（3）的结果文件夹进行小时数据至日数据的合并，结果存放入形如“###_dailly_bands_tif”和“###_dailly_count_bands_tif”文件夹下的不同band子文件夹中。
  - （5）利用【06_01_sow_harvest_date_from_weatherIndex.R】文件，统计每个格点位置的理论种植或成熟的日期。利用【06_sowToharvest_dateToday.py】文件将日期转为当年距离当前1月1日的天数。将结果保存在StarandEnd_lat_###_lon_###.csv文件中。
  - （6）利用【08_indices_based_different_phase_modified_version.py】文件，根据StarandEnd_lat_###_lon_###.csv中的经纬度从###_daily_bands_tif或###_daily_count_band_tif文件夹中来提取指定经纬度点的值，结果存放在与提取文件夹相对应的###_daily_bands_csv或###_daily_count_bands_csv中的同名子文件夹中。比如：GDD_0日数据的提取，源数据的位置为：temperature_daily_band_tif/GDD_0,提取的结果存放在temperature_daily_band_csv/GDD_0文件夹中。 
  - （7）通过【07_sow_bajie_chousui_havest_GDD.py】计算除种植和成熟以外的，其它两个“拔节（1月1日至成熟日的GDD的35%）”，“抽穗（1月1日至成熟日的GDD的57%）的日期（天数）。
  - （8）根据（5）和（6）计算的三个不同生长阶段，利用【09_cal_event.py】进行每个阶段不同极端天气和气候事件的计算。
- 复合极端天气气候事件计算：

### 3. CDO 方法
>主要思路：通过python调用cdo的cmd命令对下载的nc文件进行数据的筛选及由小时数据生成天的数据
> 然后按照经纬度列表对日数据进行数据的提取，结果存入CSV文件中以备进一步处理
> 与上述方法不同的是：（1）将按照不同区间对小时的数据进行处理和生成相应的日数据同时进行；
> （2）使用CDO命令大大提高了编码效率；（3）CDO的运行效率高于python，从而大大提高了运行速度

详细步骤如下：

- 极端天气气候事件计算：
  - （1）通过【90_cdo_precipitation.py】从total_precipitation文件下读取降水的小时数据，生成不同区间的小时和日bands数据（或统计是否在此范围中的1、0数据）;
  对应存储的文件夹分别为precipitation_hourly_band_nc和precipitation_daily_band_nc（或precipitation_hourly_count_band_nc和
  precipitation_daily_count_band_nc）
  - （2）通过【90_cdo_temperature.py】从2m_temperature文件下读取气温的小时数据，生成不同区间的小时和日bands数据（或统计是否在此范围中的1、0数据）;
  对应存储的文件夹分别为temperature_hourly_band_nc和temperature_daily_band_nc（或temperature_hourly_count_band_nc和
  temperature_daily_count_band_nc）
      >   注意：
      >  + 【Ubuntu版本22.04，cdo版本2.0.4和Ubuntu版本20.04或18.04，cdo版本1.9.3】  这两个版本对应的函数：two_sides_band_stat，max_side_band_stat 和min_side_band_stat这三个函数略有区别（低版本不支持maxc等函数）；
      >  + 大批量数据处理时，尽可能让中间过程文件存入不同的硬盘，提高速度(即：out_hourly_band和out_hourly_count_band 路径的设定)。比如同时处理GDD_0和GDD_step1_1_10时，这二者的这两个文件夹分别存入两个盘，物理上达到这个盘并行运行的目的。
  - （3）通过【91_nc_to_csv.py】将上一步生成的nc文件根据经纬度列表提取数据，并将结果按年保存在csv中。所有的经纬度点一年一个csv文件（表头为：['ID','lat','lon','year','D_001','D_002',……,'D_364','D_365','D_366']）
  - （4）根据【90_cdo_hourly_RH.py,90_cdo_hourly_temp.py和90_cdo_hourly_wind10_to_wind2.py】生成的不同band的小时数据，通过【90_cdo_hourly_【temp_rh_wind】.py】来计算"气温、湿度和风"复合极端时间小时数据。（temp_rh_uv2_hourly_band_nc）
  - （5）根据【90_cdo_dailly_【temp_rh_wind】.py】生成不同band的天数据。（temp_rh_uv2_hourly_band_nc）
  - （6）重复第上述第（3）步，生成规定格式的csv文件
  - （7）调用【91_csv_merge.py】将所有年的csv文件合并为一个csv文件，名称为1981-2020.csv
  - （8）通过【90_cdo_sow_bajie_chousui_havest_GDD.py】计算每个站点每年的bajie和chousui的时间（sow和harvest时间根据【90_cdo_sow_harvest_date_from_weatherIndex.R】获得）
  - （9）通过[92_cal_extremes.py]计算每个阶段不同极端天气和气候事件。
  




### 4. CDO 操作实例

---
> The "clip" function is used to filter values between 10 and 20. Values below 10 and above 20 will be set to 10 and 20, respectively.
> ```
> cdo expr,'myvar=clip(myvar,10,20)' input.nc output.nc
> ```
- ----
> The "if" function is used to set the values greater than 5 to "missing", effectively filtering them out from the output.
> ```
> cdo expr,'myvar=if(myvar>5, myvar, missing)' input.nc output.nc
> ```
----
> This command filters the "temperature" variable in "input.nc" and saves the result in "output.nc," where values greater than 25°C are replaced with NaN.
>> - **setmissval**	Set a new missing value 
>> - **setctomiss**	Set constant to missing value
>> - **setmisstoc**	Set missing value to constant
>> - **setrtomiss**	Set range to missing value
>> - **setvrange**	Set valid range
>> - **setmisstonn**	Set missing value to nearest neighbor
>> - **setmisstodis**	Set missing value to distance-weighted average
>```
> cdo ifthen -gtc,25 input.nc -setmissval,NAN output.nc
> cdo ifthen -ltc,0 input.nc -setmissval,NAN output.nc
> ```
---
>Suppose you have multiple NetCDF files named "file1.nc", "file2.nc", "file3.nc", and so on, and you want to merge them into a single file named "merged.nc".
>```
>cdo merge file1.nc file2.nc file3.nc merged.nc
>```
---BAND---
>Suppose you have a NetCDF file named "input.nc" with two bands (variables) "band1" and "band2", and you want to add these two bands together to create a new variable "sum_band" and save it in a new file "output.nc"
>```
>cdo expr,'sum_band = band1 + band2' input.nc output.nc
>```






import os
import cdo

# Initialize the Cdo class
cdo = cdo.Cdo()

# Set the path to the folder containing the nc files
input_folder = 'path/to/your/folder/'

# Get a list of all nc files in the folder
nc_files = [file for file in os.listdir(input_folder) if file.endswith('.nc')]

# Loop through each file and perform the operations
for nc_file in nc_files:
    input_file = os.path.join(input_folder, nc_file)
    output_file = os.path.join(input_folder, 'filtered_' + nc_file)  # Prefix the output file with 'filtered_'

    # Operation 1: Keep values between 0 and 5, set others to null values
    cdoExpr1 = f'setmisstoc,5 -setrtomiss,0,5 {input_file} {output_file}'
    cdo(cdoExpr1.split())

    # Operation 2: Convert the hour band to the day band
    cdoExpr2 = f'daymean {output_file} {output_file}'  #daysum
	# cdoExpr = f'timselmean,24 {input_file} {output_file}'
    cdo(cdoExpr2.split())
	
	timselsum	Time selection sum
    timselmean	Time selection mean
	
	expression='ifthennc,glt(tp,5) -setrtomiss,0,5 input_file output_file'
	cdoExpr=f'{expresion}{input_file.nc}{out_file}'


import os
import cdo

# Initialize the Cdo class
cdo = cdo.Cdo()
a=[key for key, value in cdo.operators.items()]