| Hosted by CoCalc | Download
Kernel: Python 2 (system-wide)
from math import* import matplotlib.pyplot as plt import numpy as np % matplotlib inline from scipy.interpolate import * from scipy.interpolate import interp1d from scipy.integrate import quad
# Latitude and Longitude points for Utah border lat = [39.5014905, 39.682817, 40.108472, 40.544004, 41.001601, 40.996484, 41.001601, 41.001601, 41.402501, 41.75998, 42.005364, 41.9972, 41.9972, 41.9972, 41.993117, 41.775152, 41.532997, 41.149448, 40.652448, 40.200838, 39.5014905, 38.686515, 37.958207, 37.012354, 37.004607, 37.004607, 37.000222, 36.997617, 37.002004, 37.694552, 38.162501, 38.586283, 39.5014904] print(len(lat)) lng = [-109.0448, -109.05304, -109.047546, -109.05304, -109.05304, -109.731445, -110.475769, -111.047058, -111.041565, -111.047058, -111.041565, -111.55242899, -112.61261, -113.431091, -114.035339, -114.032593, -114.038086, -114.071045, -114.062805, -114.051819, -114.040833, -114.046326, -114.035339, -114.035339, -113.065796, -111.55242899, -110.703735, -109.577637, -109.072266, -109.055786, -109.033813, -109.050293, -109.0448] print(len(lng))
33 33
# finding the midpoint mlat = (max(lat)-min(lat))/(2) + min(lat) mlng = (max(lng)-min(lng))/(2) + min(lng) midpoint = (mlat, mlng) print(midpoint) plt.plot(lng,lat)
(39.5014905, -111.55242899999999)
[<matplotlib.lines.Line2D at 0x7fa1e9b00450>]
/ext/sage/sage-8.1/local/lib/python2.7/site-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment. warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')
Image in a Jupyter notebook
xvec = [] yvec = [] thetavec = [] rvec2 = [] for i in range(len(lat)): # determining x coordinate in km from midpoint xpt = (lat[i] - mlat)*111.0 xvec.append(xpt) print('Distance in km between latitudes') print(np.transpose(xvec)) print('') for i in range(len(lng)): ypt = (lng[i]-mlng)*111.0 # determining y coordinate in km from midpoint yvec.append(ypt) print('Distance in km between longitudes') print(np.transpose(yvec)) print('') for i in range(len(yvec)-1): i = i+1 theta = np.arctan2(xvec[i],yvec[i]) # determining theta using inverse tangent of (x,y) cooridnates if theta < 0: theta = theta + np.pi*2.0 thetavec.append(theta) for i in range(len(xvec)-1): i=i+1 rdist = np.sqrt(abs(xvec[i])**2 + abs(yvec[i])**2) # determining r value rvec2.append(rdist) print('Theta') print(np.transpose(thetavec)) print('') print('radius from midpoint to each point') print(np.transpose(rvec2))
Distance in km between latitudes [ 0.00000000e+00 2.01272415e+01 6.73749465e+01 1.15718998e+02 1.66512265e+02 1.65944278e+02 1.66512265e+02 1.66512265e+02 2.11012165e+02 2.50692334e+02 2.77929958e+02 2.77023754e+02 2.77023754e+02 2.77023754e+02 2.76570541e+02 2.52376426e+02 2.25497221e+02 1.82923282e+02 1.27756282e+02 7.76275725e+01 0.00000000e+00 -9.04622805e+01 -1.71304469e+02 -2.76294152e+02 -2.77154069e+02 -2.77154069e+02 -2.77640804e+02 -2.77929959e+02 -2.77443002e+02 -2.00570174e+02 -1.48627835e+02 -1.01588033e+02 -1.11000001e-05] Distance in km between longitudes [ 2.78346819e+02 2.77432179e+02 2.78042013e+02 2.77432179e+02 2.77432179e+02 2.02129224e+02 1.19509260e+02 5.60961810e+01 5.67059040e+01 5.60961810e+01 5.67059040e+01 1.10999930e-06 -1.17680091e+02 -2.08531482e+02 -2.75603010e+02 -2.75298204e+02 -2.75907927e+02 -2.79566376e+02 -2.78651736e+02 -2.77432290e+02 -2.76212844e+02 -2.76822567e+02 -2.75603010e+02 -2.75603010e+02 -1.67983737e+02 1.10999930e-06 9.42050340e+01 2.19201912e+02 2.75298093e+02 2.77127373e+02 2.79566376e+02 2.77737096e+02 2.78346819e+02] Theta [ 0.07242145 0.23773681 0.39516655 0.54055985 0.68740403 0.94827617 1.24584893 1.3082659 1.35065764 1.36952929 1.57079632 1.97249878 2.21605938 2.35444227 2.39960642 2.45639665 2.562193 2.71170907 2.86876268 3.14159265 3.45744085 3.69771617 3.92824312 4.16749452 4.71238898 5.03950476 5.38019839 5.49390667 5.65669874 5.79454957 5.93253072 6.28318527] radius from midpoint to each point [ 278.16131973 286.08870025 300.59857045 323.56598787 261.52194317 204.9604786 175.70747305 218.49872663 256.89186071 283.65581499 277.0237545 300.98299682 346.73843101 390.44626205 373.4741782 356.33436697 334.09322931 306.54275019 288.08803437 276.212844 291.22870359 324.50306632 390.25053142 324.0878177 277.1540685 293.18766038 353.96968805 390.84991888 342.09351845 316.61900097 295.7330263 278.346819 ]
x = thetavec y = rvec2 f2 = interp1d(x,y, kind='cubic') x2 = np.linspace(x[0],x[len(y)-1],1000) y2 = f2(x2) plt.plot(x,y,'r') + plt.plot(x2,y2,'b')
[<matplotlib.lines.Line2D at 0x7fa1e98973d0>, <matplotlib.lines.Line2D at 0x7fa1e9897510>]
Image in a Jupyter notebook
def f(x): return f2 def Mapping(x,a,b): # Mapping the interval from [-1,1] y = ((b-a)/2.0)*(x+1.0) + a return y
def ThreePtQuadAtB(a,b,f=f): # using quadrature to approximate the area # works on [a,b] x1 = Mapping(0,a,b) x2 = Mapping((-0.2)*np.sqrt(15),a,b) x3 = Mapping((0.2)*np.sqrt(15),a,b) w1 = ((b-a)/2.0)*(0.8888888889) w2 = ((b-a)/2.0)*(0.5555555556) w3 = ((b-a)/2.0)*(0.5555555556) Area = w1*f(x1) + w2*f(x2) + w3*f(x3) return Area
Area = 0.0 for i in range(len(x)): A = ThreePtQuadAtB(x[i],x[len(y)-1],f2) # approxiamtion using 3-pt quadrature rule Area = Area + A print(Area)
33027.9680498