前段時間NASA的API更新了，修改了一些參數，現在重新上線。

NASA有專門為作物生長模型提供的氣象數據，主要包括輻射、溫度、降雨、風速。其分辨率為0.5度，通過輸入經緯度和時間段獲取數據。
獲取數據的基本函數（可以不用看）
主要修改main函數里的
start_date = dt.date(2000,1,1)
end_date = dt.date(2020,2,2)
latitude=34.5
longitude=112.5
四個參數，輸出的結果是Excel可以直接被WOFOST模型讀取，其中空值用了前后5日的平均值代替。

# 獲取數據的基礎函數 import pandas as pd import numpy as np import requests import datetime as dt from math import expMJ_to_KJ = lambda x: x * 1e3 mm_to_cm = lambda x: x / 10. tdew_to_kpa = lambda x: ea_from_tdew(x) / 10 * 10. to_date = lambda d: d.date()def getnasadata(latitude, longitude, start_date, end_date):server = "https://power.larc.nasa.gov/api/temporal/daily/point"# Variable names in POWER datapower_variables = ["TOA_SW_DWN", "ALLSKY_SFC_SW_DWN", "T2M", "T2M_MIN","T2M_MAX", "T2MDEW", "WS2M", "PRECTOTCORR"]payload = {"request": "execute","parameters": ",".join(power_variables),"latitude": latitude,"longitude": longitude,"start": start_date.strftime("%Y%m%d"),"end": end_date.strftime("%Y%m%d"),"community": "AG","format": "JSON","user": "anonymous"}req = requests.get(server, params=payload)return req.json()def _process_POWER_records(powerdata):"""Process the meteorological records returned by NASA POWER"""fill_value = float(powerdata["header"]["fill_value"])power_variables = ["TOA_SW_DWN", "ALLSKY_SFC_SW_DWN", "T2M", "T2M_MIN","T2M_MAX", "T2MDEW", "WS2M", "PRECTOTCORR"]df_power = {}for varname in power_variables:s = pd.Series(powerdata["properties"]["parameter"][varname])s[s == fill_value] = np.NaNdf_power[varname] = sdf_power = pd.DataFrame(df_power)df_power["DAY"] = pd.to_datetime(df_power.index, format="%Y%m%d")# find all rows with one or more missing values (NaN)ix = df_power.isnull().any(axis=1)# Get all rows without missing valuesdf_power = df_power[~ix]return df_powerdef _estimate_AngstAB(df_power):"""Determine Angstrom A/B parameters from Top-of-Atmosphere (ALLSKY_TOA_SW_DWN) andtop-of-Canopy (ALLSKY_SFC_SW_DWN) radiation values.:param df_power: dataframe with POWER data:return: tuple of Angstrom A/B valuesThe Angstrom A/B parameters are determined by dividing swv_dwn by toa_dwnand taking the 0.05 percentile for Angstrom A and the 0.98 percentile forAngstrom A+B: toa_dwn*(A+B) approaches the upper envelope whiletoa_dwn*A approaches the lower envelope of the records of swv_dwnvalues."""# check if sufficient data is available to make a reasonable estimate:# As a rule of thumb we want to have at least 200 days availableangstA = 0.29angstB = 0.49if len(df_power) < 200:msg = ("Less then 200 days of data available. Reverting to " +"default Angstrom A/B coefficients (%f, %f)")print(msg)return angstA, angstB# calculate relative radiation (swv_dwn/toa_dwn) and percentilesrelative_radiation = df_power.ALLSKY_SFC_SW_DWN / df_power.TOA_SW_DWNix = relative_radiation.notnull()angstrom_a = float(np.percentile(relative_radiation[ix].values, 5))angstrom_ab = float(np.percentile(relative_radiation[ix].values, 98))angstrom_b = angstrom_ab - angstrom_aMIN_A = 0.1MAX_A = 0.4MIN_B = 0.3MAX_B = 0.7MIN_SUM_AB = 0.6MAX_SUM_AB = 0.9A = abs(angstrom_a)B = abs(angstrom_b)SUM_AB = A + Bif A < MIN_A or A > MAX_A or B < MIN_B or B > MAX_B or SUM_AB < MIN_SUM_AB or SUM_AB > MAX_SUM_AB:msg = ("Angstrom A/B values (%f, %f) outside valid range " +"Reverting to default values.")msg = msg % (angstrom_a, angstrom_b)print(msg)return angstA, angstBreturn angstrom_a, angstrom_bdef _POWER_to_PCSE(df_power):# Convert POWER data to a dataframe with PCSE compatible inputsdf_pcse = pd.DataFrame({"DAY": df_power.DAY.apply(to_date),"IRRAD": df_power.ALLSKY_SFC_SW_DWN.apply(MJ_to_KJ),"TMIN": df_power.T2M_MIN,"TMAX": df_power.T2M_MAX,"VAP": df_power.T2MDEW.apply(tdew_to_kpa),"WIND": df_power.WS2M,"RAIN": df_power.PRECTOTCORR})return df_pcsedef ea_from_tdew(tdew):# Raise exception:if (tdew < -95.0 or tdew > 65.0):# Are these reasonable bounds?msg = 'tdew=%g is not in range -95 to +60 deg C' % tdewraise ValueError(msg)tmp = (17.27 * tdew) / (tdew + 237.3)ea = 0.6108 * exp(tmp)return ea# 單點獲取 def main():start_date = dt.date(2000, 1, 1)end_date = dt.date(2020, 2, 2)latitude = 35.5longitude = 112.5powerdata = getnasadata(latitude, longitude, start_date, end_date)if not powerdata:msg = "Failure retrieving POWER data from server. This can be a connection problem with " \"the NASA POWER server, retry again later."print(msg)description = [powerdata["header"]["title"]]elevation = float(powerdata["geometry"]["coordinates"][2])df_power = _process_POWER_records(powerdata)# Determine Angstrom A/B parametersangstA, angstB = _estimate_AngstAB(df_power)# Convert power records to PCSE compatible structuredf_pcse = _POWER_to_PCSE(df_power)# 去除空值，填充空值，加上降雪SNOWDEPTH = [-999 for i in range((end_date - start_date).days + 1)]std_datatime = pd.DataFrame(pd.date_range(start=start_date, end=end_date, freq='D'), columns=['DAY'])data = df_pcse.dropna()data['DAY'] = pd.to_datetime(data['DAY'])result = pd.merge(std_datatime, data, how='left', on='DAY')result['SNOWDEPTH'] = SNOWDEPTH# 找出缺失值索引nan_indexset = set(np.where(np.isnan(result.iloc[:, 1:]))[0])nan_index = list(nan_indexset)miss_value = len(nan_index)nan_index.sort()# 替換空值for i in nan_index:result.loc[i, 'IRRAD'] = result.iloc[i - 5:i + 6, 1].mean()result.loc[i, 'TMIN'] = result.iloc[i - 5:i + 6, 2].mean()result.loc[i, 'TMAX'] = result.iloc[i - 5:i + 6, 3].mean()result.loc[i, 'VAP'] = result.iloc[i - 5:i + 6, 4].mean()result.loc[i, 'WIND'] = result.iloc[i - 5:i + 6, 5].mean()result.loc[i, 'RAIN'] = result.iloc[i - 5:i + 6, 6].mean()# 寫文件表頭excelhead = pd.DataFrame({'DAY': ['Site Characteristics', 'Country', 'Station', 'Description', 'Source', 'Contact','Missing values', 'Longitude', longitude, 'Observed data', 'DAY', 'date'],'IRRAD': [np.NaN, 'China', 'miss_value={}days'.format(miss_value), description,'Meteorology and Air Quality Group, Wageningen University', 'Peter Uithol',-999, 'Latitude', latitude, np.NaN, 'IRRAD', 'kJ/m2/day or hours'],'TMIN': [np.NaN, np.NaN, np.NaN, np.NaN, np.NaN, np.NaN, np.NaN, 'Elevation', elevation,np.NaN, 'TMIN', 'Celsius'],'TMAX': [np.NaN, np.NaN, np.NaN, np.NaN, np.NaN, np.NaN, np.NaN, 'AngstromA', angstA,np.NaN, 'TMAX', 'Celsius'],'VAP': [np.NaN, np.NaN, np.NaN, np.NaN, np.NaN, np.NaN, np.NaN, 'AngstromB', angstB,np.NaN, 'VAP', 'kPa'],'WIND': [np.NaN, np.NaN, np.NaN, np.NaN, np.NaN, np.NaN, np.NaN, 'HasSunshine', False,np.NaN, 'WIND', 'm/sec'],'RAIN': [np.NaN, np.NaN, np.NaN, np.NaN, np.NaN, np.NaN, np.NaN, np.NaN, np.NaN, np.NaN,'RAIN', 'mm'],'SNOWDEPTH': [np.NaN, np.NaN, np.NaN, np.NaN, np.NaN, np.NaN, np.NaN, np.NaN, np.NaN,np.NaN, 'SNOWDEPTH', 'cm']})dataexcel = pd.concat([excelhead, result], axis=0, ignore_index=True)dataexcel.to_excel(r'C:\Users\Administrator\Desktop\NASA天氣文件lat={},lon={}.xlsx'.format(latitude, longitude),index=False, header=None)print('getNASA天氣文件lat={},lon={}.xlsx successful'.format(latitude, longitude))if __name__ == '__main__':main()

參考資料

[1] https://github.com/ajwdewit/pcse/tree/master/pcse

總結

以上是生活随笔為你收集整理的从NASA获取全球气象数据的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。