inStations = r'I:\Documents\Klusjes\GIDS Interpolation Air Temperature\GSOD//NEARNLSTATIONS_MERGED_8_monday.csv'
folderGSOD = r'I:\Documents\Klusjes\GIDS Interpolation Air Temperature\GSOD2009\selection_nl_2'
RasterDEM = r'I:\Documents\Klusjes\GIDS Interpolation Air Temperature\GSOD//AALBERS_25KM_DEM.tif'
outFolder = r'I:\Documents\Klusjes\GIDS Interpolation Air Temperature\GSOD2009\outTifs3//'
prefix = 'WDSP'
date_all = [20120723, 20120724, 20120725, 20120726, 20120727, 20120728, 20120729, 20120730]
from scipy import stats
from scipy.spatial.distance import cdist
from scipy.spatial import cKDTree as KDTree
import statsmodels.api as sm
import numpy as np
import pandas as pd
import gdal
import os
from __future__ import division
def GSODfiles(inGSODFolder):
st_wmo = [os.path.join(root, name)
for root, dirs, files in os.walk(inGSODFolder)
for name in files
if name.endswith('.op')]
return st_wmo
def inRaster(fileDEM):
raster = gdal.Open(fileDEM, gdal.GA_ReadOnly)
band = raster.GetRasterBand(1)
dem = band.ReadAsArray()
extent = raster.GetGeoTransform()
return raster, dem, extent
def saveRaster(path, array, datatype=6, formatraster="GTiff"):
# Set Driver
format_ = formatraster #save as format
driver = gdal.GetDriverByName( format_ )
driver.Register()
# Set Metadata for Raster output
cols = raster.RasterXSize
rows = raster.RasterYSize
bands = raster.RasterCount
datatype = 6#band.DataType
# Set Projection for Raster
outDataset = driver.Create(outFilename, cols, rows, bands, datatype)
geoTransform = raster.GetGeoTransform()
outDataset.SetGeoTransform(geoTransform)
proj = raster.GetProjection()
outDataset.SetProjection(proj)
# Write output to band 1 of new Raster
outBand = outDataset.GetRasterBand(1)
outBand.WriteArray(array) #save input array
#outBand.WriteArray(dem)
# Close and finalise newly created Raster
#F_M01 = None
outBand = None
proj = None
geoTransform = None
outDataset = None
driver = None
datatype = None
bands = None
rows = None
cols = None
driver = None
array = None
def IDWKDtree(x,y):
x_y = [(x,y)] # for kd_tree that starts counting at bottom left with 0,0
#print x_y
x_ = (x*extent[1]+extent[0]+extent[1]/2) # longitude aalbers projection (meters)
y_ = (y*extent[5]+extent[3]+extent[5]/2) # latitude aalbers projection (meters)
long_lat = np.array([[x_,y_]])
#print long_lat
dem_1 = dem[y,x] # elevation x_y coordinate
dist_tree, ix_tree = tree.query(x_y, k=8, eps=0, p=1) # returns distance and index
dist_tree = dist_tree.flatten()
ix_tree = ix_tree.flatten()
df_selection = df.ix[ix_tree.ravel()] # seleclt 8 neighbours
c = df_selection.ix[:,14]
if dist_tree[0] < 1e-10:
wz = c.ix[ix_tree[0]]
else:
w = 1 / dist_tree**2
w /= np.sum(w)
wz = np.dot(w, c.ix[ix_tree])
#df_selection = df.ix[ix_tree.ravel()]
return wz
st_wmo = GSODfiles(folderGSOD)
raster, dem, extent = inRaster(RasterDEM)
for date in date_all:
df = pd.read_csv(inStations)
#df_select = None
df_select = pd.DataFrame()
for st in st_wmo:
with open(st) as f:
"""
0,6 = STNN--- : USAF number
14,22 = YEARMODA : YearMonthDay
24,30 = TEMP : Mean temperature in Fahreinheit
35,41 = DEWP : Mean dew point in Fahreinheit
57,63 = STP : Mean station pressure in millibars
78,83 = WDSP : Mean wind speed in knots
118,123 = PRCP : Total precipitation in inches
"""
colspecs_data = [(0,6), (14,22), (24,30), (35,41), (57,63), (78,83), (118,123)]
df_date = pd.read_fwf(f, colspecs=colspecs_data)
for moda in df_date.YEARMODA:
#print moda
if moda == date:
df_select = df_select.append(df_date[df_date.YEARMODA==date])
break
# Set Missing values
df_select.replace({'STP':{9999.9:np.nan}}, inplace=True)
df_select.replace({'TEMP':{9999.9:np.nan}}, inplace=True)
df_select.replace({'DEWP':{9999.9:np.nan}}, inplace=True)
df_select.replace({'WDSP':{999.9:np.nan}}, inplace=True)
df_select.replace({'PRCP':{99.99:0}}, inplace=True)
#print df_select.head()
# Merge values with stations
df = pd.merge(df, df_select, how='inner', left_on='USAF', right_on='STN---')
#print df.head()
##STP
#df = df[pd.notnull(df['STP'])]
##TEMP
#df = df[pd.notnull(df['TEMP'])]
#print df.head()
##DEWP
#df = df[pd.notnull(df['DEWP'])]
#df.reset_index(drop=True, inplace=True)
#print df.head()
##WDSP
df = df[pd.notnull(df['WDSP'])]
df.reset_index(drop=True, inplace=True)
##PRCP
#df = df[pd.notnull(df['PRCP'])]
Longscaled = df.ix[:,8]
#print 'longscaled\n', Longscaled
Latscaled = df.ix[:,7]
#print '\nlatscaled\n', Latscaled
tree = KDTree(zip(Longscaled,Latscaled), leafsize=11)
tp = np.zeros([dem.shape[1],dem.shape[0]])
#x = 0
#y = 0
#GIDS(x,y)
for x in range(0,20,1):
for y in range(0,19,1):
tp[x][y] = IDWKDtree(x,y)
tp = tp.T
#save output as raster
outFilename = outFolder+prefix+str(date)+'.tif'
saveRaster(outFilename, tp)
outFilename = None
tp = None
df.shape
(93, 16)