{ "metadata": { "name": "Interpolation Windspeed" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "heading", "level": 4, "metadata": {}, "source": [ "Import modules" ] }, { "cell_type": "code", "collapsed": false, "input": [ "inStations = r'I:\\Documents\\Klusjes\\GIDS Interpolation Air Temperature\\GSOD//NEARNLSTATIONS_MERGED_8_monday.csv'\n", "folderGSOD = r'I:\\Documents\\Klusjes\\GIDS Interpolation Air Temperature\\GSOD2009\\selection_nl_2'\n", "RasterDEM = r'I:\\Documents\\Klusjes\\GIDS Interpolation Air Temperature\\GSOD//AALBERS_25KM_DEM.tif'\n", "outFolder = r'I:\\Documents\\Klusjes\\GIDS Interpolation Air Temperature\\GSOD2009\\outTifs3//'\n", "prefix = 'WDSP'\n", "date_all = [20120723, 20120724, 20120725, 20120726, 20120727, 20120728, 20120729, 20120730]" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 25 }, { "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": 26 }, { "cell_type": "heading", "level": 4, "metadata": {}, "source": [ "Import weather-stations who have data on 20090101" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def GSODfiles(inGSODFolder):\n", " st_wmo = [os.path.join(root, name)\n", " for root, dirs, files in os.walk(inGSODFolder)\n", " for name in files\n", " if name.endswith('.op')]\n", " return st_wmo" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 27 }, { "cell_type": "code", "collapsed": false, "input": [ "def inRaster(fileDEM):\n", " raster = gdal.Open(fileDEM, gdal.GA_ReadOnly)\n", " band = raster.GetRasterBand(1)\n", " dem = band.ReadAsArray()\n", " extent = raster.GetGeoTransform()\n", " return raster, dem, extent" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 28 }, { "cell_type": "code", "collapsed": false, "input": [ "def saveRaster(path, array, datatype=6, formatraster=\"GTiff\"):\n", " # Set Driver\n", " format_ = formatraster #save as format\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 = 6#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(array) #save input array\n", " #outBand.WriteArray(dem)\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", " driver = None\n", " array = None" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 29 }, { "cell_type": "code", "collapsed": false, "input": [ "def IDWKDtree(x,y): \n", " x_y = [(x,y)] # for kd_tree that starts counting at bottom left with 0,0\n", " #print x_y\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=8, eps=0, p=1) # returns distance and index\n", " dist_tree = dist_tree.flatten()\n", " ix_tree = ix_tree.flatten()\n", " df_selection = df.ix[ix_tree.ravel()] # seleclt 8 neighbours\n", " c = df_selection.ix[:,14]\n", " \n", " if dist_tree[0] < 1e-10:\n", " wz = c.ix[ix_tree[0]]\n", " else:\n", " w = 1 / dist_tree**2\n", " w /= np.sum(w)\n", " wz = np.dot(w, c.ix[ix_tree])\n", " \n", " #df_selection = df.ix[ix_tree.ravel()]\n", "\n", " return wz" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 30 }, { "cell_type": "code", "collapsed": false, "input": [ "st_wmo = GSODfiles(folderGSOD)\n", "raster, dem, extent = inRaster(RasterDEM)\n", "for date in date_all:\n", "\n", " df = pd.read_csv(inStations)\n", " #df_select = None\n", " df_select = pd.DataFrame()\n", " for st in st_wmo:\n", " with open(st) as f:\n", " \"\"\"\n", " 0,6 = STNN--- : USAF number\n", " 14,22 = YEARMODA : YearMonthDay\n", " 24,30 = TEMP : Mean temperature in Fahreinheit\n", " 35,41 = DEWP : Mean dew point in Fahreinheit\n", " 57,63 = STP : Mean station pressure in millibars\n", " 78,83 = WDSP : Mean wind speed in knots\n", " 118,123 = PRCP : Total precipitation in inches\n", "\n", " \"\"\"\n", " colspecs_data = [(0,6), (14,22), (24,30), (35,41), (57,63), (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", "\n", " # Set Missing values\n", " df_select.replace({'STP':{9999.9:np.nan}}, inplace=True)\n", " df_select.replace({'TEMP':{9999.9:np.nan}}, inplace=True)\n", " df_select.replace({'DEWP':{9999.9:np.nan}}, inplace=True)\n", " df_select.replace({'WDSP':{999.9:np.nan}}, inplace=True)\n", " df_select.replace({'PRCP':{99.99:0}}, inplace=True)\n", "\n", " #print df_select.head()\n", " # Merge values with stations\n", " df = pd.merge(df, df_select, how='inner', left_on='USAF', right_on='STN---')\n", " #print df.head()\n", "\n", " ##STP\n", " #df = df[pd.notnull(df['STP'])]\n", "\n", " ##TEMP\n", " #df = df[pd.notnull(df['TEMP'])]\n", " #print df.head()\n", "\n", " ##DEWP\n", " #df = df[pd.notnull(df['DEWP'])]\n", " #df.reset_index(drop=True, inplace=True)\n", " #print df.head()\n", "\n", " ##WDSP\n", " df = df[pd.notnull(df['WDSP'])]\n", " df.reset_index(drop=True, inplace=True)\n", "\n", " ##PRCP\n", " #df = df[pd.notnull(df['PRCP'])]\n", "\n", " Longscaled = df.ix[:,8]\n", " #print 'longscaled\\n', Longscaled\n", " Latscaled = df.ix[:,7]\n", " #print '\\nlatscaled\\n', Latscaled\n", " tree = KDTree(zip(Longscaled,Latscaled), leafsize=11)\n", "\n", " tp = np.zeros([dem.shape[1],dem.shape[0]])\n", "\n", " #x = 0\n", " #y = 0\n", " #GIDS(x,y)\n", " for x in range(0,20,1):\n", " for y in range(0,19,1):\n", " tp[x][y] = IDWKDtree(x,y)\n", " tp = tp.T\n", "\n", " #save output as raster\n", " outFilename = outFolder+prefix+str(date)+'.tif'\n", " saveRaster(outFilename, tp)\n", " outFilename = None\n", " tp = None" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 31 }, { "cell_type": "code", "collapsed": false, "input": [ "df.shape" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "pyout", "prompt_number": 32, "text": [ "(93, 16)" ] } ], "prompt_number": 32 }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 32 } ], "metadata": {} } ] }