ulvis.paste.net

Paste Search Dynamic
Recent pastes
idwr
  1. #!/usr/bin/env python
  2. # packages
  3. import math
  4. import numpy as np
  5. import csv
  6. #------------------------------------------------------------
  7. # Distance calculation, degree to km (Haversine method)
  8. def harvesine(lon1, lat1, lon2, lat2):
  9.     rad = math.pi / 180  # degree to radian
  10.     R = 6378.1  # earth average radius at equador (km)
  11.     dlon = (lon2 - lon1) * rad
  12.     dlat = (lat2 - lat1) * rad
  13.     a = (math.sin(dlat / 2)) ** 2 + math.cos(lat1 * rad) * \
  14.         math.cos(lat2 * rad) * (math.sin(dlon / 2)) ** 2
  15.     c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))
  16.     d = R * c
  17.     return(d)
  18. # ------------------------------------------------------------
  19. # Prediction
  20. def idwr(x, y, z, xi, yi):
  21.     lstxyzi = []
  22.     for p in range(len(xi)):
  23.         lstdist = []
  24.         for s in range(len(x)):
  25.             d = (harvesine(x[s], y[s], xi[p], yi[p]))
  26.             lstdist.append(d)
  27.         sumsup = list((1 / np.power(lstdist, 2)))
  28.         suminf = np.sum(sumsup)
  29.         sumsup = np.sum(np.array(sumsup) * np.array(z))
  30.         u = sumsup / suminf
  31.         xyzi = [xi[p], yi[p], u]
  32.         lstxyzi.append(xyzi)
  33.     return(lstxyzi)
  34.  
  35. # know points
  36. #x = [-47.6, -48.9, -48.2, -48.9, -47.6, -48.6]
  37. #y = [-23.4, -24.0, -23.9, -23.1, -22.7, -22.5]
  38. #z = [27.0,  33.4,  34.6,  18.2,   30.8, 42.8]
  39. # running the function
  40. # output
  41. #[[-48.05306, -23.59167, 31.486682779040855]]
  42.  
  43. x = []
  44. y = []
  45. z = []
  46.  
  47.  
  48. with open('hoy.csv') as csv_file:
  49.     csv_reader = csv.reader(csv_file, delimiter=',')
  50.     line_count = 0
  51.     for row in csv_reader:
  52.         x.append(float(row[20]))
  53.         y.append(float(row[19]))
  54.         z.append(float(row[4]))
  55.         line_count += 1
  56.  
  57. def interpola(xi,yi):
  58.      return idwr(x,y,z,xi,yi)[0][2]
  59. xi = [-3.71224722222223]
  60. yi = [40.4238527777779]
  61.  
  62. print interpola(xi,yi)
Parsed in 0.016 seconds