SharedUTAH.ipynbOpen in CoCalc
Jupyter notebook UTAH.ipynb
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.')
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>]
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