python 站点插值学习

############站点数据插值为网格数据  通过给定的经纬度以及变量值插值为网格数据

############注意插值的网格分辨率，太密会造成数据出现异常值

############根据站点密度设置网格密度  

############   

import pandas as pd

import numpy as np

############读取文件

f = pd.read_csv('/home/mw/input/moyu1828/1957_2012_tm_year_avg.txt',

                 sep='\s+',names=['sta','lat','lon','alt','year','tm'],na_values=-99.90)

########names 对于每一列命名       na_VALUES 设置缺测值

#########由于有的站点不是每年都记录数据

#########如果提取多年，写循环读入

data = f.loc[f.year==2006]                  ##########首先索引出2006年数据

data = data.dropna(axis=0,how='any')

#####data.shape(看数据的行列数？数据形态)

########剃除缺测值

########通过使用dropna进行操作axis=0（0是行，1是列）；；；；how='any'这一行or行含有任意一个缺测值就踢除

#############提取数据纬度、经度、年平均温度、站号

lat = data.lat.values/100          ###########正常经纬度

lon = data.lon.values/100

tm = data.tm.values

sta = data.sta.values

print(data.agg(['min', 'max']))           ###########通过agg函数看min/max

##############cressman 插值

from metpy.interpolate import inverse_distance_to_grid

lon_grid = np.arange(75.0, 134.0, 1.0)  #lon   加一个步长

lat_grid = np.arange(16.0, 54.0, 1.0)   #lat

###print

lon_gridmesh, lat_gridmesh = np.meshgrid(lon_grid, lat_grid)############lon*lat

tm_grid = inverse_distance_to_grid(lon,lat,tm,lon_gridmesh,lat_gridmesh,r=15,min_neighbors=3)

######help（） 帮助看所需函数需要的变量

####插值tm数据(lon,lat, 插值前的经纬度坐标 tm,lon_gridmesh,lat_gridmesh,r=15 ###搜索半径内的站点数

#####          min_neighbors=3 最小插值点数 )

import matplotlib.pyplot as plt

import cartopy.crs as ccrs

import cartopy.feature as cfeature

import cartopy.mpl.ticker as cticker

fig = plt.figure(figsize=(12,8))

ax1 = fig.add_axes([0.1, 0.1, 0.8, 0.4],projection = ccrs.PlateCarree(central_longitude=90))

leftlon, rightlon, lowerlat, upperlat = (0,180,0,90)

img_extent = [leftlon, rightlon, lowerlat, upperlat]

ax1.set_extent(img_extent, crs=ccrs.PlateCarree())

# ax1.add_feature(cfeature.COASTLINE.with_scale('50m')) 

# ax1.add_feature(cfeature.LAKES, alpha=0.5)

ax1.set_xticks(np.arange(60,150+30,30), crs=ccrs.PlateCarree())

ax1.set_yticks(np.arange(0,60+30,30), crs=ccrs.PlateCarree())

lon_formatter = cticker.LongitudeFormatter()

lat_formatter = cticker.LatitudeFormatter()

ax1.xaxis.set_major_formatter(lon_formatter)

ax1.yaxis.set_major_formatter(lat_formatter)

c1 = ax1.scatter(lon,lat,s=2, zorder=0,c=tm,transform=ccrs.PlateCarree(), cmap=plt.cm.Reds)

position=fig.add_axes([0.35, 0.02,  0.35, 0.025])

fig.colorbar(c1,cax=position,orientation='horizontal',format='%d',)