Radiosonde Location Map of Picked Sites


In [1]:
%matplotlib inline

Visakhpatnam to Kandahar


In [136]:
#############
# Count number of strings in list that are the same
#
#
#############

#station_list_pick=[43150, 42867, 43014, 42339, 40990, 40948, 41780, 43003, 42182, 42369, 43353, 43295, 43279, 43369, 42867, 42339, 43128, 43285, 43150, 43185, 43371, 43346, 42492, 43192, 42379, 43014, 42027]
station_list_pick=[43150, 42867, 43014, 42339, 40990, 40948]

from collections import defaultdict
import os
import numpy as np

import re

#import matplotlib
#matplotlib.use('Agg') # Must be before importing matplotlib.pyplot or pylab!

import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

station_data=np.load('/nfs/a90/eepdw/Data/Observations/Radiosonde_downloaded_from_NOAA_GUAN/Embrace_Period_India_Station_and_sounding_Info_measured_derived.npy')

#match_bad = re.compile(r'9999')
station_list='/nfs/a90/eepdw/Data/Observations/Radiosonde_downloaded_from_NOAA_GUAN/igra-stations.txt'



st_lat_c=[]
st_lon_c=[]
st_nam_c=[]
st_cnt_c=[]



for i, c in enumerate(station_data):
 if c[3]:
  
  match_station= re.compile(r"%s" % c[3][0][0][0:5])
  if match_station.search(str(station_list_pick)) is not None:    
       #print st_lat_c
       st_lat_c.append(c[1])
       st_lon_c.append(c[2])
       st_nam_c.append(c[0].lower().title())
        
       #st_cnt_c[i]=str()
       #st_namcnt_c[i]=st_nam_c[i] + ': ' + st_cnt_c[i]
       #st_namcnt_c.append(st_nam_c[i])
    
       #print station_data
#PLOT TIME AND DATE vS FREQ ALL STATIONS
#print st_nam_c
#print st_lon_c
# create figure and axes instances

slope, offset=np.polyfit(np.array(st_lon_c, dtype=float), np.array(st_lat_c, dtype=float), 1)
y_line=slope*np.array(st_lon_c, dtype=float)+offset

lon_1=np.array(st_lon_c, dtype=float)[-1]
lon_2=np.array(st_lon_c, dtype=float)[0]
lat_1=y_line[-1]
lat_2=y_line[0]

initial_bearing,compass_bearing=calculate_initial_compass_bearing((lat_1,lon_1), (lat_2,lon_2))
maths_bearing = compass_to_maths_bearing(initial_bearing)
print 'Compass bearing - %s' % compass_bearing
print 'atan2 bearing - %s' % initial_bearing
print 'Maths bearing - %s' % maths_bearing
print 'In terms of wind rotations, using the maths bearing will rotate the u-winds to be in line with the transect'
print 'I want the u winds (x-axis) to be perpendicular to the transect, so will use the the maths bearing-90,'
print 'which is the same as the negative of the atan2 bearing - %s' % (-initial_bearing)

fig = plt.figure(figsize=(16,16))
ax = fig.add_axes([0.1,0.1,0.8,0.8])
# create polar stereographic Basemap instance.
#m = Basemap(projection='ste',lon_0=lon_mid,lat_0=lat_mid,lat_ts=lat_mid,\

m = Basemap(projection='cyl',\
            llcrnrlat=-10.,urcrnrlat=40.,\
            llcrnrlon=60.,urcrnrlon=100.,\
            rsphere=6371200.,resolution='l',area_thresh=10000)
# draw coastlines, state and country boundaries, edge of map.
m.drawcoastlines()
# draw parallels.
parallels = np.arange(0.,90,10.)
m.drawparallels(parallels,labels=[1,0,0,0],fontsize=10)
# draw meridians
meridians = np.arange(0.,360.,10.)
m.drawmeridians(meridians,labels=[0,0,0,1],fontsize=10)


x, y = m(st_lon_c,st_lat_c)
m.scatter(x,y,3,marker='o',color='red')
   
for i,j,s in zip(x, y, st_nam_c):
 
 plt.text(i, j, '$%s$' % s, fontsize=14)

plt.title('Position and number of soundings in August and September 2011')
#plt.savefig('/nfs/a90/eepdw/Figures/Observation_Plots/sounding_station_map_igra_guan_favourites.png',  format='png', bbox_inches='tight')

m.plot(x, y_line)
plt.show()

for x in sorted(st_nam_c):
    print x


Compass bearing - 319.946441716
atan2 bearing - -40.0535582838
Maths bearing - 130.053558284
In terms of wind rotations, using the maths bearing will rotate the u-winds to be in line with the transect
I want the u winds (x-axis) to be perpendicular to the transect, so will use the the maths bearing-90,
which is the same as the negative of the atan2 bearing - 40.0535582838
Aurangabad/Chikalthana
Jodhpur
Kabul_Airport
Kandahar_Airport
Nagpur_Sonegaon
Visakhapatnam

Malabar Coast


In [137]:
#############
# Count number of strings in list that are the same
#
#
#############

#station_list_pick=[43150, 42867, 43014, 42339, 40990, 40948, 41780, 43003, 42182, 42369, 43353, 43295, 43279, 43369, 42867, 42339, 43128, 43285, 43150, 43185, 43371, 43346, 42492, 43192, 42379, 43014, 42027]
#station_list_pick=[43150, 42867, 43014, 42339, 40990, 40948]
station_list_pick=[43371, 43353, 43285,  43192, 43003]
# Thiruvanthanapuram, Cochin, Panambur, Panjim, Bombay
from collections import defaultdict
import os
import numpy as np

import re

#import matplotlib
#matplotlib.use('Agg') # Must be before importing matplotlib.pyplot or pylab!

import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

station_data=np.load('/nfs/a90/eepdw/Data/Observations/Radiosonde_downloaded_from_NOAA_GUAN/Embrace_Period_India_Station_and_sounding_Info_measured_derived.npy')

#match_bad = re.compile(r'9999')
station_list='/nfs/a90/eepdw/Data/Observations/Radiosonde_downloaded_from_NOAA_GUAN/igra-stations.txt'



st_lat_c=[]
st_lon_c=[]
st_nam_c=[]
st_cnt_c=[]



for i, c in enumerate(station_data):
 if c[3]:
  
  match_station= re.compile(r"%s" % c[3][0][0][0:5])
  if match_station.search(str(station_list_pick)) is not None:    
       #print st_lat_c
       st_lat_c.append(c[1])
       st_lon_c.append(c[2])
       st_nam_c.append(c[0].lower().title())
        
       #st_cnt_c[i]=str()
       #st_namcnt_c[i]=st_nam_c[i] + ': ' + st_cnt_c[i]
       #st_namcnt_c.append(st_nam_c[i])
    
       #print station_data
#PLOT TIME AND DATE vS FREQ ALL STATIONS
#print st_nam_c
#print st_lon_c
# create figure and axes instances

slope, offset=np.polyfit(np.array(st_lon_c, dtype=float), np.array(st_lat_c, dtype=float), 1)
y_line=slope*np.array(st_lon_c, dtype=float)+offset

lon_1=np.array(st_lon_c, dtype=float)[-1]
lon_2=np.array(st_lon_c, dtype=float)[0]
lat_1=y_line[-1]
lat_2=y_line[0]

initial_bearing,compass_bearing=calculate_initial_compass_bearing((lat_1,lon_1), (lat_2,lon_2))
maths_bearing = compass_to_maths_bearing(initial_bearing)
print 'Compass bearing - %s' % compass_bearing
print 'atan2 bearing - %s' % initial_bearing
print 'Maths bearing - %s' % maths_bearing
print 'In terms of wind rotations, using the maths bearing will rotate the u-winds to be in line with the transect'
print 'I want the u winds (x-axis) to be perpendicular to the transect, so will use the the maths bearing-90,'
print 'which is the same as the negative of the atan2 bearing - %s' % (-initial_bearing)


fig = plt.figure(figsize=(16,16))
ax = fig.add_axes([0.1,0.1,0.8,0.8])
# create polar stereographic Basemap instance.
#m = Basemap(projection='ste',lon_0=lon_mid,lat_0=lat_mid,lat_ts=lat_mid,\

m = Basemap(projection='cyl',\
            llcrnrlat=-10.,urcrnrlat=40.,\
            llcrnrlon=60.,urcrnrlon=100.,\
            rsphere=6371200.,resolution='l',area_thresh=10000)
# draw coastlines, state and country boundaries, edge of map.
m.drawcoastlines()
# draw parallels.
parallels = np.arange(0.,90,10.)
m.drawparallels(parallels,labels=[1,0,0,0],fontsize=10)
# draw meridians
meridians = np.arange(0.,360.,10.)
m.drawmeridians(meridians,labels=[0,0,0,1],fontsize=10)


x, y = m(st_lon_c,st_lat_c)
m.scatter(x,y,3,marker='o',color='red')
   
for i,j,s in zip(x, y, st_nam_c):
 
 plt.text(i, j, '$%s$' % s, fontsize=14)

plt.title('Position and number of soundings in August and September 2011')
#plt.savefig('/nfs/a90/eepdw/Figures/Observation_Plots/sounding_station_map_igra_guan_favourites.png',  format='png', bbox_inches='tight')

m.plot(x, y_line)
plt.show()

for x in sorted(st_nam_c):
    print x


Compass bearing - 339.401572226
atan2 bearing - -20.5984277743
Maths bearing - 110.598427774
In terms of wind rotations, using the maths bearing will rotate the u-winds to be in line with the transect
I want the u winds (x-axis) to be perpendicular to the transect, so will use the the maths bearing-90,
which is the same as the negative of the atan2 bearing - 20.5984277743
Bombay/Santa_Cruz
Cochin/Willingdon
Panambur
Panjim
Thiruvananthapuram

Bombay to Nagpur


In [138]:
#############
# Count number of strings in list that are the same
#
#
#############

#station_list_pick=[43150, 42867, 43014, 42339, 40990, 40948, 41780, 43003, 42182, 42369, 43353, 43295, 43279, 43369, 42867, 42339, 43128, 43285, 43150, 43185, 43371, 43346, 42492, 43192, 42379, 43014, 42027]
#station_list_pick=[43150, 42867, 43014, 42339, 40990, 40948]
station_list_pick=[43003, 43014, 42867]
# Bombay, Aurangabd, Nagpur
from collections import defaultdict
import os
import numpy as np

import re

#import matplotlib
#matplotlib.use('Agg') # Must be before importing matplotlib.pyplot or pylab!

import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

station_data=np.load('/nfs/a90/eepdw/Data/Observations/Radiosonde_downloaded_from_NOAA_GUAN/Embrace_Period_India_Station_and_sounding_Info_measured_derived.npy')

#match_bad = re.compile(r'9999')
station_list='/nfs/a90/eepdw/Data/Observations/Radiosonde_downloaded_from_NOAA_GUAN/igra-stations.txt'



st_lat_c=[]
st_lon_c=[]
st_nam_c=[]
st_cnt_c=[]



for i, c in enumerate(station_data):
 if c[3]:
  
  match_station= re.compile(r"%s" % c[3][0][0][0:5])
  if match_station.search(str(station_list_pick)) is not None:    
       #print st_lat_c
       st_lat_c.append(c[1])
       st_lon_c.append(c[2])
       st_nam_c.append(c[0].lower().title())
        
       #st_cnt_c[i]=str()
       #st_namcnt_c[i]=st_nam_c[i] + ': ' + st_cnt_c[i]
       #st_namcnt_c.append(st_nam_c[i])
    
       #print station_data
#PLOT TIME AND DATE vS FREQ ALL STATIONS
#print st_nam_c
#print st_lon_c
# create figure and axes instances

slope, offset=np.polyfit(np.array(st_lon_c, dtype=float), np.array(st_lat_c, dtype=float), 1)
y_line=slope*np.array(st_lon_c, dtype=float)+offset

lon_1=np.array(st_lon_c, dtype=float)[-1]
lon_2=np.array(st_lon_c, dtype=float)[0]
lat_1=y_line[-1]
lat_2=y_line[0]

initial_bearing,compass_bearing=calculate_initial_compass_bearing((lat_1,lon_1), (lat_2,lon_2))
maths_bearing = compass_to_maths_bearing(initial_bearing)
print 'Compass bearing - %s' % compass_bearing
print 'atan2 bearing - %s' % initial_bearing
print 'Maths bearing - %s' % maths_bearing
print 'In terms of wind rotations, using the maths bearing will rotate the u-winds to be in line with the transect'
print 'I want the u winds (x-axis) to be perpendicular to the transect, so will use the the maths bearing-90,'
print 'which is the same as the negative of the atan2 bearing - %s' % (-initial_bearing)

fig = plt.figure(figsize=(16,16))
ax = fig.add_axes([0.1,0.1,0.8,0.8])
# create polar stereographic Basemap instance.
#m = Basemap(projection='ste',lon_0=lon_mid,lat_0=lat_mid,lat_ts=lat_mid,\

m = Basemap(projection='cyl',\
            llcrnrlat=-10.,urcrnrlat=40.,\
            llcrnrlon=60.,urcrnrlon=100.,\
            rsphere=6371200.,resolution='l',area_thresh=10000)
# draw coastlines, state and country boundaries, edge of map.
m.drawcoastlines()
# draw parallels.
parallels = np.arange(0.,90,10.)
m.drawparallels(parallels,labels=[1,0,0,0],fontsize=10)
# draw meridians
meridians = np.arange(0.,360.,10.)
m.drawmeridians(meridians,labels=[0,0,0,1],fontsize=10)


x, y = m(st_lon_c,st_lat_c)
m.scatter(x,y,3,marker='o',color='red')
   
for i,j,s in zip(x, y, st_nam_c):
 
 plt.text(i, j, '$%s$' % s, fontsize=14)

plt.title('Position and number of soundings in August and September 2011')
#plt.savefig('/nfs/a90/eepdw/Figures/Observation_Plots/sounding_station_map_igra_guan_favourites.png',  format='png', bbox_inches='tight')

m.plot(x, y_line)
plt.show()

for x in sorted(st_nam_c):
    print x


Compass bearing - 70.4543001112
atan2 bearing - 70.4543001112
Maths bearing - 19.5456998888
In terms of wind rotations, using the maths bearing will rotate the u-winds to be in line with the transect
I want the u winds (x-axis) to be perpendicular to the transect, so will use the the maths bearing-90,
which is the same as the negative of the atan2 bearing - -70.4543001112
Aurangabad/Chikalthana
Bombay/Santa_Cruz
Nagpur_Sonegaon

In [122]:
def calculate_initial_compass_bearing(pointA, pointB):
    """
    Calculates the bearing between two points.
 
    The formulae used is the following:
    θ = atan2(sin(Δlong).cos(lat2),
    cos(lat1).sin(lat2) − sin(lat1).cos(lat2).cos(Δlong))
 
    :Parameters:
    - `pointA: The tuple representing the latitude/longitude for the
    first point. Latitude and longitude must be in decimal degrees
    - `pointB: The tuple representing the latitude/longitude for the
    second point. Latitude and longitude must be in decimal degrees
 
    :Returns:
    The bearing in degrees
 
    :Returns Type:
    float
    """
    if (type(pointA) != tuple) or (type(pointB) != tuple):
        raise TypeError("Only tuples are supported as arguments")
 
    lat1 = math.radians(pointA[0])
    lat2 = math.radians(pointB[0])
 
    diffLong = math.radians(pointB[1] - pointA[1])
 
    y = math.sin(diffLong) * math.cos(lat2)
    x = math.cos(lat1) * math.sin(lat2) - (math.sin(lat1)
    * math.cos(lat2) * math.cos(diffLong))
 
    initial_bearing = math.atan2(y, x)
 
    # Now we have the initial bearing but math.atan2 return values
    # from -180° to + 180° which is not what we want for a compass bearing
    # The solution is to normalize the initial bearing as shown below
    initial_bearing = math.degrees(initial_bearing)
    compass_bearing = (initial_bearing + 360) % 360
 
    return initial_bearing, compass_bearing

def compass_to_maths_bearing(initial_bearing):
    
    '''
    Converts degrees positive clockwise from north (in range +/- 180, or 0-360) bearing to
    positive anticlockwise from east maths bearing
    '''
    if initial_bearing > 180:
        initial_bearing=initial_bearing-360
        
    maths_bearing=(-initial_bearing)+90
    return maths_bearing

In [132]:
import numpy as np
from math import *
import math

slope, offset=np.polyfit(np.array(st_lon_c, dtype=float), np.array(st_lat_c, dtype=float), 1)
y_line=slope*np.array(st_lon_c, dtype=float)+offset

lon_1=np.array(st_lon_c, dtype=float)[-1]
lon_2=np.array(st_lon_c, dtype=float)[0]
lat_1=y_line[-1]
lat_2=y_line[0]

initial_bearing,compass_bearing=calculate_initial_compass_bearing((lat_1,lon_1), (lat_2,lon_2))
maths_bearing = compass_to_maths_bearing(initial_bearing)
print 'Compass bearing - %s' % compass_bearing
print 'atan2 bearing - %s' % initial_bearing
print 'Maths bearing - %s' % maths_bearing
#u_rot = U * cos(brng) + V * sin(brng)
#v_rot = -U * sin(brng) + V * cos(brng)


Compass bearing - 70.4543001112
atan2 bearing - 70.4543001112
Maths bearing - 19.5456998888
Out[132]:
267.2403843294163

θ = atan2( sin Δλ ⋅ cos φ2 , cos φ1 ⋅ sin φ2 − sin φ1 ⋅ cos φ2 ⋅ cos Δλ )

Bearing

In general, your current heading will vary as you follow a great circle path (orthodrome); the final heading will differ from the initial heading by varying degrees according to distance and latitude (if you were to go from say 35°N,45°E (≈ Baghdad) to 35°N,135°E (≈ Osaka), you would start on a heading of 60° and end up on a heading of 120°!).

This formula is for the initial bearing (sometimes referred to as forward azimuth) which if followed in a straight line along a great-circle arc will take you from the start point to the end point:1 Formula: θ = atan2( sin Δλ ⋅ cos φ2 , cos φ1 ⋅ sin φ2 − sin φ1 ⋅ cos φ2 ⋅ cos Δλ ) JavaScript: (all angles in radians)

var y = Math.sin(λ2-λ1) Math.cos(φ2); var x = Math.cos(φ1)Math.sin(φ2) - Math.sin(φ1)Math.cos(φ2)Math.cos(λ2-λ1); var brng = Math.atan2(y, x).toDegrees();

Excel: (all angles in radians) =ATAN2(COS(lat1)SIN(lat2)-SIN(lat1)COS(lat2)COS(lon2-lon1), SIN(lon2-lon1)COS(lat2)) *note that Excel reverses the arguments to ATAN2 – see notes below

Since atan2 returns values in the range -π ... +π (that is, -180° ... +180°), to normalise the result to a compass bearing (in the range 0° ... 360°, with −ve values transformed into the range 180° ... 360°), convert to degrees and then use (θ+360) % 360, where % is (floating point) modulo.

For final bearing, simply take the initial bearing from the end point to the start point and reverse it (using θ = (θ+180) % 360).

N.B - lats and lons have to be converted to radians first