python气象信息系统工程 第七章

第7章 气象数据插值

7.1 空间插值

7.1.1 从站点到栅格

1
2
3
import pandas as pd
df = pd.read_csv('/home/mw/input/pythonbook4259/r160.txt', sep='\s+', names=['lat', 'lon', 'pre_ano'])
print(df)
       lat     lon   pre_ano
0    51.72  126.65  0.098331
1    48.77  121.92 -0.167619
2    49.22  119.75 -0.064597
3    50.45  121.70 -0.046360
4    49.17  125.23  0.017577
..     ...     ...       ...
155  47.73   88.08 -0.036620
156  46.73   83.00  0.071207
157  44.43   84.66 -0.205051
158  43.95   81.33  0.170275
159  43.78   87.62 -0.103539

[160 rows x 3 columns]

1
2
3
4
5
6
7
8
9
from pykrige.ok import OrdinaryKriging
krige = OrdinaryKriging(
          df['lon'],
          df['lat'],
          df['pre_ano'],
          variogram_model="linear",
          verbose=False,
          enable_plotting=False,
       )
---------------------------------------------------------------------------

ModuleNotFoundError                       Traceback (most recent call last)

<ipython-input-2-137613763ad2> in <module>
----> 1 from pykrige.ok import OrdinaryKriging
      2 krige = OrdinaryKriging(
      3           df['lon'],
      4           df['lat'],
      5           df['pre_ano'],


ModuleNotFoundError: No module named 'pykrige'

1
2
3
4
5
import numpy as np
lon = np.arange(76.0, 131.5+0.5, 0.5) # np.arange()函数创建的数组不包含参数传入的终点值
lat = np.arange(20.5, 51.5+0.5, 0.5) # 所以在终点值后再加入一个步长以确保最后一个值被包含在结果中
print(lon.shape, lat.shape)

1
2
3
pre_grid, ss = krige.execute("grid", lon, lat)
print(pre_grid)

1
2
3
4
import xarray as xr
pre_da = xr.DataArray(pre_grid, coords=[lat, lon], dims=['lat', 'lon'])
print(pre_da)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
import numpy as np
from metpy.interpolate import inverse_distance_to_grid

lon = np.arange(76.0, 131.5+0.5, 0.5)
lat = np.arange(20.5, 51.5+0.5, 0.5)
lon_grid, lat_grid = np.meshgrid(lon, lat) # 生成交织的二维经纬栅格
# 注意：inverse_distance_to_grid()使用的目标栅格为交织后的二维经纬栅格
pre_grid = inverse_distance_to_grid(
              df['lon'].values,
              df['lat'].values,
              df['pre_ano'].values,
              lon_grid,
              lat_grid,
              r=15,
              min_neighbors=3
          )

print(pre_grid)

[[        nan         nan         nan ... -0.41675891 -0.41730563
  -0.41772805]
 [        nan         nan         nan ... -0.41783573 -0.41883586
  -0.42000272]
 [        nan         nan         nan ... -0.41931729 -0.41991869
  -0.42128663]
 ...
 [-0.04594149 -0.04876271 -0.05241695 ... -0.15684901 -0.15684024
  -0.15693899]
 [-0.04011599 -0.04204267 -0.04612476 ... -0.1527561  -0.15290136
  -0.15305936]
 [-0.03320475 -0.03582687 -0.03862931 ... -0.14872947 -0.14882568
  -0.14891597]]

1
2
3
4
5
6
import xarray as xr

pre_da = xr.DataArray(pre_grid, coords=[lat, lon], dims=['lat', 'lon'])

print(pre_da)

<xarray.DataArray (lat: 63, lon: 112)>
array([[        nan,         nan,         nan, ..., -0.41675891,
        -0.41730563, -0.41772805],
       [        nan,         nan,         nan, ..., -0.41783573,
        -0.41883586, -0.42000272],
       [        nan,         nan,         nan, ..., -0.41931729,
        -0.41991869, -0.42128663],
       ...,
       [-0.04594149, -0.04876271, -0.05241695, ..., -0.15684901,
        -0.15684024, -0.15693899],
       [-0.04011599, -0.04204267, -0.04612476, ..., -0.1527561 ,
        -0.15290136, -0.15305936],
       [-0.03320475, -0.03582687, -0.03862931, ..., -0.14872947,
        -0.14882568, -0.14891597]])
Coordinates:
  * lat      (lat) float64 20.5 21.0 21.5 22.0 22.5 ... 49.5 50.0 50.5 51.0 51.5
  * lon      (lon) float64 76.0 76.5 77.0 77.5 78.0 ... 130.0 130.5 131.0 131.5

7.1.2 从栅格到站点

1
print(pre_da)
<xarray.DataArray (lat: 63, lon: 112)>
array([[        nan,         nan,         nan, ..., -0.41675891,
        -0.41730563, -0.41772805],
       [        nan,         nan,         nan, ..., -0.41783573,
        -0.41883586, -0.42000272],
       [        nan,         nan,         nan, ..., -0.41931729,
        -0.41991869, -0.42128663],
       ...,
       [-0.04594149, -0.04876271, -0.05241695, ..., -0.15684901,
        -0.15684024, -0.15693899],
       [-0.04011599, -0.04204267, -0.04612476, ..., -0.1527561 ,
        -0.15290136, -0.15305936],
       [-0.03320475, -0.03582687, -0.03862931, ..., -0.14872947,
        -0.14882568, -0.14891597]])
Coordinates:
  * lat      (lat) float64 20.5 21.0 21.5 22.0 22.5 ... 49.5 50.0 50.5 51.0 51.5
  * lon      (lon) float64 76.0 76.5 77.0 77.5 78.0 ... 130.0 130.5 131.0 131.5

1
2
pre_sta = pre_da.interp(lon=lon, lat=lat, method='linear')
print(pre_sta)
<xarray.DataArray (lat: 63, lon: 112)>
array([[        nan,         nan,         nan, ..., -0.41675891,
        -0.41730563, -0.41772805],
       [        nan,         nan,         nan, ..., -0.41783573,
        -0.41883586, -0.42000272],
       [        nan,         nan,         nan, ..., -0.41931729,
        -0.41991869, -0.42128663],
       ...,
       [-0.04594149, -0.04876271, -0.05241695, ..., -0.15684901,
        -0.15684024, -0.15693899],
       [-0.04011599, -0.04204267, -0.04612476, ..., -0.1527561 ,
        -0.15290136, -0.15305936],
       [-0.03320475, -0.03582687, -0.03862931, ..., -0.14872947,
        -0.14882568, -0.14891597]])
Coordinates:
  * lon      (lon) float64 76.0 76.5 77.0 77.5 78.0 ... 130.0 130.5 131.0 131.5
  * lat      (lat) float64 20.5 21.0 21.5 22.0 22.5 ... 49.5 50.0 50.5 51.0 51.5

1
2
pre_df = pre_sta.to_dataframe('pre_ano')
print(pre_df)
             pre_ano
lat  lon            
20.5 76.0        NaN
     76.5        NaN
     77.0        NaN
     77.5        NaN
     78.0        NaN
...              ...
51.5 129.5 -0.148536
     130.0 -0.148608
     130.5 -0.148729
     131.0 -0.148826
     131.5 -0.148916

[7056 rows x 1 columns]

7.1.3 从栅格到栅格

1
print(pre_da)
<xarray.DataArray (lat: 63, lon: 112)>
array([[        nan,         nan,         nan, ..., -0.41675891,
        -0.41730563, -0.41772805],
       [        nan,         nan,         nan, ..., -0.41783573,
        -0.41883586, -0.42000272],
       [        nan,         nan,         nan, ..., -0.41931729,
        -0.41991869, -0.42128663],
       ...,
       [-0.04594149, -0.04876271, -0.05241695, ..., -0.15684901,
        -0.15684024, -0.15693899],
       [-0.04011599, -0.04204267, -0.04612476, ..., -0.1527561 ,
        -0.15290136, -0.15305936],
       [-0.03320475, -0.03582687, -0.03862931, ..., -0.14872947,
        -0.14882568, -0.14891597]])
Coordinates:
  * lat      (lat) float64 20.5 21.0 21.5 22.0 22.5 ... 49.5 50.0 50.5 51.0 51.5
  * lon      (lon) float64 76.0 76.5 77.0 77.5 78.0 ... 130.0 130.5 131.0 131.5

1
2
3
4
5
6
pre_grid1 = pre_da.interp(lon=[80, 90, 100, 120], 
                      lat=[25, 30, 35, 40], 
                      method='linear')

print(pre_sta)

<xarray.DataArray (lat: 63, lon: 112)>
array([[        nan,         nan,         nan, ..., -0.41675891,
        -0.41730563, -0.41772805],
       [        nan,         nan,         nan, ..., -0.41783573,
        -0.41883586, -0.42000272],
       [        nan,         nan,         nan, ..., -0.41931729,
        -0.41991869, -0.42128663],
       ...,
       [-0.04594149, -0.04876271, -0.05241695, ..., -0.15684901,
        -0.15684024, -0.15693899],
       [-0.04011599, -0.04204267, -0.04612476, ..., -0.1527561 ,
        -0.15290136, -0.15305936],
       [-0.03320475, -0.03582687, -0.03862931, ..., -0.14872947,
        -0.14882568, -0.14891597]])
Coordinates:
  * lon      (lon) float64 76.0 76.5 77.0 77.5 78.0 ... 130.0 130.5 131.0 131.5
  * lat      (lat) float64 20.5 21.0 21.5 22.0 22.5 ... 49.5 50.0 50.5 51.0 51.5

7.2 时间插值

7.2.1 站点时间内插

1
2
3
4
5
6
7
8
9
10
import pandas as pd
#构造数组
df = pd.DataFrame([[21.7, 983, 0.64],
                  [19.2, 991, 0.75],
                  [13.4, 973, 0.83]],
                 columns=['t', 'p', 'rh']
                )
df['time'] = pd.to_datetime(['2020-02-19', '2020-02-20', '2020-02-22'])
print(df)

      t    p    rh       time
0  21.7  983  0.64 2020-02-19
1  19.2  991  0.75 2020-02-20
2  13.4  973  0.83 2020-02-22

1
2
3
df = df.set_index('time')
print(df)

               t    p    rh
time                       
2020-02-19  21.7  983  0.64
2020-02-20  19.2  991  0.75
2020-02-22  13.4  973  0.83

1
2
3
df_interp = df.resample('1D').interpolate(method='linear')
print(df_interp)

               t      p    rh
time                         
2020-02-19  21.7  983.0  0.64
2020-02-20  19.2  991.0  0.75
2020-02-21  16.3  982.0  0.79
2020-02-22  13.4  973.0  0.83

7.2.2 栅格时间内插

1
2
3
4
5
6
7
import pandas as pd
import xarray as xr
import numpy as np
#构造数组
da = xr.DataArray(np.abs(np.random.randn(3, 3, 3)*5+15),coords=[pd.to_datetime(['2020-02-19', '2020-02-20', '2020-02-22']),[20, 25, 30], [120, 125, 130]],dims=['time', 'lat', 'lon'])
print(da)

<xarray.DataArray (time: 3, lat: 3, lon: 3)>
array([[[25.36833749,  9.835558  ,  9.70410351],
        [16.96870264,  7.33170443, 14.86629016],
        [16.39984725, 16.46497263, 16.3571344 ]],

       [[14.5675103 , 12.17993797,  9.57918095],
        [23.61240778, 14.57624162, 10.1265022 ],
        [12.1024726 , 10.50885673, 11.87143659]],

       [[ 7.60884169, 20.78943309,  7.90046148],
        [18.41886334, 23.20704066, 17.45145241],
        [ 5.89131442, 19.3040832 , 14.7046148 ]]])
Coordinates:
  * time     (time) datetime64[ns] 2020-02-19 2020-02-20 2020-02-22
  * lat      (lat) int64 20 25 30
  * lon      (lon) int64 120 125 130

1
2
3
da_interp = da.interp(time=['2020-02-21'], method='linear')
print(da_interp)

<xarray.DataArray (time: 1, lat: 3, lon: 3)>
array([[[11.08817599, 16.48468553,  8.73982122],
        [21.01563556, 18.89164114, 13.78897731],
        [ 8.99689351, 14.90646997, 13.28802569]]])
Coordinates:
  * lat      (lat) int64 20 25 30
  * lon      (lon) int64 120 125 130
  * time     (time) datetime64[ns] 2020-02-21
/小书匠/python