{ "metadata": { "name": "Nalder Wein Real Data-Function-GSOD AIR TEMP_Iteration" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "heading", "level": 4, "metadata": {}, "source": [ "Import modules" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from scipy import stats\n", "from scipy.spatial.distance import cdist\n", "from scipy.spatial import cKDTree as KDTree\n", "import statsmodels.api as sm\n", "import numpy as np\n", "import pandas as pd\n", "import gdal\n", "import os\n", "from __future__ import division" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 1 }, { "cell_type": "heading", "level": 4, "metadata": {}, "source": [ "Import weather-stations who have data on 20090101" ] }, { "cell_type": "code", "collapsed": false, "input": [ "st_wmo = [os.path.join(root, name)\n", " for root, dirs, files in os.walk(r'I:\\Documents\\Klusjes\\GIDS Interpolation Air Temperature\\GSOD2009\\selection_nl')\n", " for name in files\n", " if name.endswith('.op')]" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 2 }, { "cell_type": "code", "collapsed": false, "input": [ "date_all = [20120501]#, 20120502, 20120503, 20120504, 20120505, 20120506, 20120507, 20120508, 20120509, 20120510, 20120511, 20120512, 20120513, 20120514, 20120515, 20120516, 20120517, 20120518, 20120519, 20120520, 20120521, 20120522, 20120523, 20120524, 20120525, 20120526, 20120527, 20120528, 20120529, 20120530]" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 3 }, { "cell_type": "code", "collapsed": false, "input": [ "file_raster = r'I:\\Documents\\Klusjes\\GIDS Interpolation Air Temperature\\GSOD//AALBERS_25KM_DEM.tif'\n", "raster = gdal.Open(file_raster, gdal.GA_ReadOnly)\n", "band = raster.GetRasterBand(1)\n", "dem = band.ReadAsArray()\n", "#print band.GetNoDataValue()\n", "#print 'Grid to evaluate\\n', dem" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 4 }, { "cell_type": "code", "collapsed": false, "input": [ "def GIDS(x,y): \n", " x_y = [(x,y)] # for kd_tree that starts counting at bottom left with 0,0\n", " \n", " x_ = (x*extent[1]+extent[0]+extent[1]/2) # longitude aalbers projection (meters)\n", " y_ = (y*extent[5]+extent[3]+extent[5]/2) # latitude aalbers projection (meters)\n", " long_lat = np.array([[x_,y_]])\n", " \n", " #print long_lat\n", " \n", " dem_1 = dem[y,x] # elevation x_y coordinate\n", " dist_tree, ix_tree = tree.query(x_y, k=10, eps=0, p=1) # returns distance and index\n", " df_selection = df.ix[ix_tree.ravel()]\n", " \n", " #print 'elevation from x_y -', dem_1\n", " #print '\\n10 nearest neighbours\\n', df_selection\n", " \n", " Longi = df_selection.ix[:,10]\n", " Lati = df_selection.ix[:,11]\n", " hi = df_selection.ix[:,5]\n", " ti = df_selection.ix[:,16]\n", " ti = (ti-32)*(5/9)\n", " \n", " pr_var = zip(Longi,Lati,hi) # combines predictor variables as tuples\n", " y = zip(ti) # dependent variable\n", " X = sm.add_constant(pr_var) # multiple linear regression\n", " \n", " #fit the model\n", " mlr = sm.OLS(y,X).fit()\n", " b1,b2,b3,b0 = mlr.params\n", " #print mlr.summary()\n", " \n", " long_lat_stations = df_selection.as_matrix(columns=['POINT_X','POINT_Y']) \n", " \n", " di = cdist(long_lat_stations, long_lat, 'euclidean') # Returns Eucleadian distance in meters between grid cell and selected weather-stations\n", " #print di\n", " \n", " # prepare datasets as flattened numpy array or as single values\n", " Hi = df_selection.as_matrix(columns=['Elev_scaled']).flatten()\n", " Ti = df_selection.as_matrix(columns=['TEMP_y']).flatten()\n", " Ti = (Ti-32)*(5/9)\n", " di_ = di.flatten()\n", " long_lat_ = long_lat.flatten()\n", " Longi_ = df_selection.as_matrix(columns=['POINT_X',]).flatten()\n", " Lati_ = df_selection.as_matrix(columns=['POINT_Y',]).flatten()\n", " \n", " top = sum( (1/di_)**2 )**-1\n", " long_f = b1*(long_lat_[0]-Longi_)\n", " lat_f = b2*(long_lat_[1]-Lati_)\n", " h_f = b3*(dem_1-Hi)\n", " middle = Ti + long_f + lat_f + h_f\n", " end = (1/di_)**2 \n", " comb = top * sum( middle * end )\n", " #print comb\n", " return comb" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 5 }, { "cell_type": "code", "collapsed": false, "input": [ "for date in date_all: \n", " df = pd.read_csv(r'I:\\Documents\\Klusjes\\GIDS Interpolation Air Temperature\\GSOD2009//20120502.csv')\n", " df_select = None\n", " df_select = pd.DataFrame()\n", " for st in st_wmo:\n", " with open(st) as f: \n", " colspecs_data = [(0,6), (14,22), (24,30), (78,83), (118,123)]\n", " df_date = pd.read_fwf(f, colspecs=colspecs_data)\n", " for moda in df_date.YEARMODA:\n", " #print moda\n", " if moda == date:\n", " df_select = df_select.append(df_date[df_date.YEARMODA==date])\n", " break\n", " df = pd.merge(df, df_select, how='inner', left_on='USAF', right_on='STN---')\n", " #print df.head()\n", " \n", " Longscaled = df.ix[:,12]\n", " Latscaled = df.ix[:,13]\n", " tree = KDTree(zip(Longscaled,Latscaled), leafsize=11)\n", " #tree.data\n", " \n", " file_raster = r'I:\\Documents\\Klusjes\\GIDS Interpolation Air Temperature\\GSOD//AALBERS_25KM_DEM.tif'\n", " raster = gdal.Open(file_raster, gdal.GA_ReadOnly)\n", " band = raster.GetRasterBand(1)\n", " dem = band.ReadAsArray()\n", " #print band.GetNoDataValue()\n", " #print 'Grid to evaluate\\n', dem\n", " \n", " extent = raster.GetGeoTransform()\n", " #print extent\n", " #print raster.GetProjection() # fyi\n", " \n", " tp = np.zeros([dem.shape[0],dem.shape[1]])\n", " tp.shape\n", " \n", " #Double for-loop\n", " for x in range(tp.shape[1]-1):\n", " for y in range(tp.shape[0]-1):\n", " tp[x][y] = GIDS(x,y)\n", " #print 'T predicted\\n', tp\n", " \n", " outFilename = r'I:\\Documents\\Klusjes\\GIDS Interpolation Air Temperature\\GSOD2009\\outTifs2//'+str(date)+'.tif'\n", " #outFilename\n", " \n", " # Set Driver\n", " format = \"GTiff\"\n", " driver = gdal.GetDriverByName( format )\n", " driver.Register()\n", " \n", " # Set Metadata for Raster output\n", " cols = raster.RasterXSize\n", " rows = raster.RasterYSize\n", " bands = raster.RasterCount\n", " datatype = band.DataType\n", " \n", " # Set Projection for Raster\n", " outDataset = driver.Create(outFilename, cols, rows, bands, datatype)\n", " geoTransform = raster.GetGeoTransform()\n", " outDataset.SetGeoTransform(geoTransform)\n", " proj = raster.GetProjection()\n", " outDataset.SetProjection(proj)\n", " \n", " # Write output to band 1 of new Raster\n", " outBand = outDataset.GetRasterBand(1)\n", " outBand.WriteArray(tp,0,0)\n", " \n", " # Close and finalise newly created Raster\n", " #F_M01 = None\n", " outBand = None\n", " proj = None\n", " geoTransform = None\n", " outDataset = None\n", " driver = None\n", " datatype = None\n", " bands = None\n", " rows = None\n", " cols = None\n", " river = None\n", " tp = None" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stderr", "text": [ "D:\\Python27\\lib\\site-packages\\statsmodels\\tools\\tools.py:306: FutureWarning: The default of `prepend` will be changed to True in 0.5.0, use explicit prepend\n", " FutureWarning)\n" ] } ], "prompt_number": 6 }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] } ], "metadata": {} } ] }