Project

Radar

All HTML code to prettify python code generated using hilite.me by Alexander Kojevnikov.
Original Code is written by Freya Lumb with citetations where relevant for formulas.
Original Code was written on a linux operating system hence the file url is'/' not '\\' as it is for windows. In addition I am using Python 2.7 as that is the only one compatitable for pyart which I need to read in the radar files.
Data sourced from Centre for Environmental Data Analysis ( CEDA ), FAAM core and core cloud and Chilbolton_Dataset, CAMRa RHI files.
I have broken the coincident code up into sections in order to explain each section clearly.

Initial imports and description

Matplotlib is Matlab style plotting for python, and is a very userfriendly plotting module which I recommend.
Pyart is the module specifically for reading in radar. More information here
The os module is very useful, the os stands for operating system, which allows the user to interface simply with the fuctions within the script
The numbers module is now moot, I used it to check if an array had just numbers. I had been having difficulties when defining a new array from a previous array and it was bringing along all its gumph - "true", "false" etc. when I just required the numbers. I think this has been fixed now!
netCDF4 module is the module to read in the netCDF4 files (.nc)

The first steps is to important the data files from their location, and simplify the names for all the aircraft and all the radar data. The loops, loop through all the files in the folder to include the complete dataset for the day.

#!/usr/bin/env python2
# -*- coding: utf-8 -*-
"""
Created on Sun Apr  2 10:56:13 2017

@author: Freya I. Lumb
Coincident Radar Script
This script is designed to point to point match 
in situ aircraft measurements of clouds with 
radar measurements of the same cloud.
The script is under progress so at present
only records how many coincident measurements
are made and not plot, or output to a file.
"""
import matplotlib.pyplot as plt  
import pyart 
from matplotlib import rcParams
import os
import numbers
import numpy as np
import netCDF4
from netCDF4 import Dataset
#Import in relevant modules for code and some for future code i.e pyart for plotting
#Radar data tends to be suffixed with "r"
#Aircraft data tends to be suffixed with "a"
filelocationa = '/home/ubo/Documents/MRes Rain and Radar/FAAM Data/03032009'
#location of the aircraft files
filelocationr = '/home/ubo/Documents/MRes Rain and Radar/CAMRa Data/03032009 RHI Data'
#location of the radar files
for filena in os.listdir(filelocationa):
#loops through folder to retrieve all aircraft files
    if filena.endswith(".nc"):
        nameoffilea = filena
        #Example of nameoffile = 'core-cloud-phy_faam_20090527_v500_r0_b449.nc'
        filenamea = filelocationa+'/'+nameoffilea
		#filename of aircraft file for rest of code
        scriptlocationa ='/home/ubo/Documents/MRes Rain and Radar/findcloud'
		#location of script
        figurelocationa ='/home/ubo/Documents/MRes Rain and Radar/FAAM Data/03032009 Fig'
		#location of figure
        aircraft = Dataset(filenamea)
		#Data from aircraft file
        
        for filenr in os.listdir(filelocationr):
		#loops through all the radar files
            if filenr.endswith(".nc"):
                nameoffiler = filenr
                filenamer = filelocationr+'/'+nameoffiler
                #filename of radar files         
                scriptlocationr ='/home/ubo/Documents/MRes Rain and Radar/Scripts'
				#location of radar scripts
                figurelocationr ='/home/ubo/Documents/MRes Rain and Radar/CAMRa Data/03032009 RHI Figures'
				#location of radar figures
                radar = pyart.aux_io.read_camra_rhi(filenamer)
				#Data from radar files

Location of radar and aircraft

longitude_a = aircraft.variables['LON_GIN'][:]
#longitude from gps in decimal degrees east
latitude_a = aircraft.variables['LAT_GIN'][:]
#latitude from gps in decimal degrees north
altitude_a = aircraft.variables['ALT_GIN'][:]
#altitude from gps in m from sea level
time_a = aircraft.variables['Time'][:]
if radar.longitude['data']<0:
#statement just corrects longitude in terms of east not west
	longitude_r = 360 + radar.longitude['data']
	#longitude of radar in decimal degrees east
else:
    longitude_r = radar.longitude['data']
	#longitude of radar in decimal degrees east
lon_ra = np.deg2rad(longitude_r)
#longitude of radar in radians
latitude_r = radar.latitude['data']
#latitude of radar in decimal degrees north
lat_ra = np.deg2rad(latitude_r)
#latitude of radar in radians
altitude_r = radar.altitude['data'] + (9*(np.sin(radar.elevation['data'])))
#height of the radar corrected for angle of beam
time_r = radar.time['data']
#time in seconds

Solid Angle

solid_angle = 1.133*radar.fields['beamwidthV']['data']*radar.fields['beamwidthH']['data']
#Solid angle is making the assumption of a gaussian function of the beam citep(www.atnf.csiro.au/people/Tobias.Westmeier/tools_spherical.php)    
#Please note that as beamwidthV = beamwidthH we have only continued with one direction for the area
volume_rb = (solid_angle/(4*np.pi))*(4/3)*np.pi*(radar.range['data']**3)
#volume of the partial sphere

Constants

R = 6371
#radius of the earth in km
R_ra = 6371 + altitude_r
rng=radar.range['data']/1000
#radar range in km
el_rad=np.deg2rad(radar.elevation['data'])
#elevation of radar in radians
nin_rad=0.5*np.pi
#0.5*pi
bwv = np.deg2rad(radar.fields['beamwidthV']['data'])
#converts beamwidth to radians

Haversine Formula

In the image, you can see the radar beam length "r" and height "h", below is the Earth with radius R. As the radar beam is a solid angle, the view of the image and lines drawn are in reference to the centre of the radar beam. What we decided to do was take the mid point of the beam, and the width horizontally. The reason for this choice was to cover the main beam area to start and if we thought it was necessary add in the full beam area with the vertical cross section at a later date. This is because as we are working in 3d, each part of the radar beam has a different location in space. The diagram shows a scalene triangle between the radar, the end of the radar beam and the center of the Earth. We can apply the cosine formula to calculate the beam's altitude. We can also use the sine formula to work out "c" which is curvature of the Earth's surface for the sector of the radar beam. Once we know "c" we can use the sector of a circle formula to calculate the distance between the radar and the end of the beam. We then use all this, and our co-ordinates of the location of the Radar to calculate the longitude and latitude reference to the end of the radar beam.

Haversine Model (Made by Freya Lumb, drawn in Solid Works courtesy of George Addison)

for i in range(0,len(rng)):
#loops through all the values
    if 0<rng[i]<250:
    #statement gets rid of unrealistic radar range values
        mid_y = np.sqrt((rng[i]**2)+(R_ra**2)-(2*rng[i]*R_ra*np.cos(nin_rad+el_rad)))
			#distance between mid point end of beam and centre of the earth in km
			#print mid_y
        mid_c=np.arcsin((rng[i]*np.sin(nin_rad+el_rad))/mid_y)
			#curvature of Earth for midpoint
        mid_s = mid_c*R
			#distance from radar to midpoint
        mid_a = mid_y-R
			#altitude of midpoint
        azi_ra = np.deg2rad(radar.azimuth['data'])
            #bearing in radians
        mid_lat0 = np.arcsin((np.sin(lat_ra)*np.cos(mid_s/R))+(np.cos(lat_ra)*np.sin(mid_s/R)*np.cos(azi_ra)))
            #mid latitude point of radar
        mid_lon0 = lon_ra + np.arctan2((np.sin(azi_ra)*np.sin(mid_s/R)*np.cos(lat_ra)),(np.cos(mid_s/R)-np.sin(lat_ra)*np.sin(mid_lat0)))
			#mid longitude point of radar
			#==============================================================================
			#Haversine Formula http://www.movable-type.co.uk/scripts/latlong.html
			#http://stackoverflow.com/questions/7222382/get-lat-long-given-current-point-distance-and-bearing
            #http://mathforum.org/library/drmath/view/51879.html
			#============================================================================== 
        mid_lat = np.rad2deg(mid_lat0)
        mid_lon = np.rad2deg(mid_lon0)
			#final mid lat and long points
                                
        max_y = np.sqrt((rng[i]**2)+(R_ra**2)-(2*rng[i]*R_ra*np.cos(nin_rad+el_rad+0.5*bwv)))
            #distance between maxmid point end of beam and centre of the earth in km
        max_c=np.arcsin((rng[i]*np.sin(nin_rad+el_rad+0.5*bwv))/max_y)
            #curvature of Earth for maxmidpoint
        max_s = max_c*R
            #distance from radar to maxmidpoint
        max_a = max_y-R
            #altitude of maxmidpoint
        max_lat0 = np.arcsin((np.sin(lat_ra)*np.cos(max_s/R))+(np.cos(lat_ra)*np.sin(max_s/R)*np.cos(azi_ra)))
            #maxmid latitude point of radar
        max_lon0 = lon_ra + np.arctan2((np.sin(azi_ra)*np.sin(max_s/R)*np.cos(lat_ra)),(np.cos(max_s/R)-np.sin(lat_ra)*np.sin(max_lat0)))
            #maxmid longitude point of radar
        max_lat = np.rad2deg(max_lat0)
        max_lon = np.rad2deg(max_lon0)
            #final maxmid lat and long points
			
        min_y = np.sqrt((rng[i]**2)+(R_ra**2)-(2*rng[i]*R_ra*np.cos(nin_rad+el_rad-0.5*bwv)))
            #distance between mid point end of beam and centre of the earth in km
        min_c=np.arcsin((rng[i]*np.sin(nin_rad+el_rad-0.5*bwv))/min_y)
            #curvature of Earth for minmidpoint
        min_s = min_c*R
            #distance from radar to minmidpoint
        min_a = min_y-R
            #altitude of minmidpoint
        min_lat0 = np.arcsin((np.sin(lat_ra)*np.cos(min_s/R))+(np.cos(lat_ra)*np.sin(min_s/R)*np.cos(azi_ra)))
            #minmid latitude point of radar
        min_lon0 = lon_ra + np.arctan2((np.sin(azi_ra)*np.sin(min_s/R)*np.cos(lat_ra)),(np.cos(min_s/R)-np.sin(lat_ra)*np.sin(min_lat0)))
            #minmid longitude point of radar
        min_lat = np.rad2deg(min_lat0)
        min_lon = np.rad2deg(min_lon0)
            #final minmid lat and long points

Location Triangle

The diagram below shows the location triangle that the beam makes with the surface of the Earth. I then applied this as a projection on a 2D plane to create an area enclosed by y=mx+c equation lines between each latitude and longitude point of which the beam lies. If the aircraft satifies all these equations, it is within the radar beam, if not, the data is no co-incident, therefore is not used.

Location Triangle (Made by Freya Lumb, drawn in Solid Works courtesy of George Addison)

delta_lat = max_lat0-min_lat0
    #difference between latitude points in radians
delta_lon = max_lon0 - min_lon0
    #difference between longitude points in radians
	
alpha = ((np.sin(delta_lat/2))**2)+(np.cos(min_lat0)*np.cos(max_lat0)*((np.sin(delta_lon/2))**2))
beta = 2*np.arctan2(np.sqrt(alpha),np.sqrt(1-alpha))
delta_d = 0.5*R*beta
                                
hmax_lat0 = np.arcsin((np.sin(mid_lat0)*np.cos(delta_d/R))+(np.cos(mid_lat0)*np.sin(delta_d/R)*np.cos(azi_ra+(np.pi/2))))
    #minmid latitude point of radar
hmax_lon0 = mid_lon0 + np.arctan2((np.sin(azi_ra+(np.pi/2))*np.sin(delta_d/R)*np.cos(mid_lat0)),(np.cos(delta_d/R)-np.sin(mid_lat0)*np.sin(hmax_lat0)))
    #minmid longitude point of radar
hmax_lat = np.rad2deg(hmax_lat0)
hmax_lon = np.rad2deg(hmax_lon0)
    #final minmid lat and long points 
                            
hmin_lat0 = np.arcsin((np.sin(mid_lat0)*np.cos(delta_d/R))+(np.cos(mid_lat0)*np.sin(delta_d/R)*np.cos(azi_ra-(np.pi/2))))
    #minmid latitude point of radar
hmin_lon0 = mid_lon0 + np.arctan2((np.sin(azi_ra-(np.pi/2))*np.sin(delta_d/R)*np.cos(mid_lat0)),(np.cos(delta_d/R)-np.sin(mid_lat0)*np.sin(hmin_lat0)))
    minmid longitude point of radar
hmin_lat = np.rad2deg(hmin_lat0)
hmin_lon = np.rad2deg(hmin_lon0)
#final minmid lat and long points 

#pologonshape, x=lon, y=lat, this isn't a perfect solution as 
m1=(latitude_r-hmin_lat)/(longitude_r-hmin_lon)
c1 = latitude_r - m1*(longitude_r)
    #y1=(m1*x1)+c1
m2=(latitude_r-hmax_lat)/(longitude_r-hmax_lon)
c2=latitude_r - m2*(longitude_r)
    #y2=(m2*x2)+c2
m3=(hmin_lat-min_lat)/(hmin_lon-min_lon)
c3=hmin_lat - m3*(hmin_lon)
    #y3=(m3*x3)+c3
m4=(hmax_lat-min_lat)/(hmax_lon-min_lon)
c4=hmin_lat - m4*(hmin_lon)
    #y4=(m4*x4)+c4

Applying the location triangle to satisfy co-incidence for time and location

Theory of the Location Triangle for Co-incident

for tb in range(0,len(timeia)):
	y1 = m1[timeir[tb]]*longitude_a[timeia[tb]]+c1[timeir[tb]]
	y2 = m2[timeir[tb]]*longitude_a[timeia[tb]]+c2[timeir[tb]]
	y3 = m3[timeir[tb]]*longitude_a[timeia[tb]]+c3[timeir[tb]]
	y4 = m4[timeir[tb]]*longitude_a[timeia[tb]]+c4[timeir[tb]]
	#location triangle equations
	if latitude_a[timeia[tb]].all() >= y1.all():
		
		if latitude_a[timeia[tb]].all() <= y2.all():
					  
			if latitude_a[timeia[tb]].all() >= y3.all():

				if latitude_a[timeia[tb]].all() >= y4.all():
				#check to satisfy all the location triangle equations
					iwctimea.append(timeia[tb])
					iwctimer.append(timeir[tb])
					for tim in range(0,len(iwctimea)):
						#update variables
						expr = radar.fields['ZED_H']['data'][iwctimer[tim]]
						zedhr.append(expr)
						#updated ZED_H for the co-incident times
						expat = aircraft['NV_TCW_U'][iwctimea[tim]]
						twca.append(expat)
						#updated TWC for the co-incident times
						expal = aircraft['NV_LWC_U'][iwctimea[tim]]
						lwca.append(expal)
						#updated LWC for the co-incident times
						iwca.append(twca[tim]-lwca[tim])
						#updated IWC for the co-incident times
						newtimes.append(tim)

Edited section of CAMRa RHI file

This code was put in for all the variable names to account for change in variable names or additional variables in the radar files over the years

try:
  ncvars['variablename']
#action (i.e print "works")
except KeyError:
  print "Variable Name does not exist"

Neely Rain Gauges

"""
Created on Wed Apr 19 10:45:23 2017

@author: Freya I. Lumb
Neely rain gauges
"""

import numpy as np

#All we need is the following mathematics everything needs to be in radians if not already
#lat0,lon0 is the coordinates in radians of the radar
#latrg,lonrg is the coordinates in radians of your rain gauge
R = 6371
#radius of the earth in km
dlat = latrg - lat0
dlon = lonrg - lon0

w = (np.sin(0.5*dlat)*np.sin(0.5*dlat)) + np.cos(lat0)*np.cos(latrg)*np.sin(0.5*dlon)*np.sin(0.5*dlon)
c = 2*np.arctan2(np.sqrt(w),np.sqrt(1-w))
#angle of curvature of Earth between radar and rain gauge (Haversine Formula)
S = R*c
#S is the sector length between the radar and the rain gauge
#el0 is the elevation of the radar beam - the beam width (BW0) to get the bottom section
j = np.pi - (el0+(0.5*np.pi)-(0.5*BW0))-c

r = (np.sin(c))*((R+h)/np.sin(j))
#using the sine rule, r is the range of the radar to the lat and lon of the rain gauge. 
#h is the height of the radar.

Multi RHI Plotting

#!/usr/bin/env python2
# -*- coding: utf-8 -*-
"""
Created on Sun Oct 23 19:16:30 2016

@author: Freya I. Lumb 
"""
import matplotlib.pyplot as plt  
import pyart 
import os

filelocation = '/home/ubo/Documents/MRes Rain and Radar/CAMRa Data/12012017 RHI Data'
#file location, insert own directory path
for filen in os.listdir(filelocation):
#for each file in the directory loop
    if filen.endswith(".nc"):
        #only use the .nc files
        
        nameoffile = filen
        filename = filelocation+'/'+nameoffile

        scriptlocation ='/home/ubo/Documents/MRes Rain and Radar/Scripts'
        #location of this code
        figurelocation ='/home/ubo/Documents/MRes Rain and Radar/CAMRa Data/12012017 RHI Figures'
        #location to put in figures
        radar = pyart.aux_io.read_camra_rhi(filename)
        #simple display object; quickest and easiest object to make plots from
        display = pyart.graph.RadarDisplay(radar)
        written()
        savefig()


def savefig():
    #Saving Figure
    # Multiple RHI Plots  and Adding A simple Correctoin to the Data
    # Define Sweep so you can change them all easily
    sweep=0  #Only 1 sweep in CAMRa files so this may be superfluous
    
    #Make Figure Really Big; Here I plan to Do 2 columns of many plots
    fig=plt.figure(figsize=(14,6))
    
    #Add First plot as axis 1 by adding a subplot
    ax1=fig.add_subplot(311)
    #add_subplot(1 row of Figures, 2 columns of Figure, This is the first plot
    # Starts in the top left)
    
     # Define Labels
    xlabel = 'Distance from Radar [km]'
    ylabel = 'Distance from Radar [km]'
    colorbar_label = 'Horizontal Reflectivity [dBZ]'
    
    #Define Colormap Lot of Options
    # We shoudl define a certain set
    plot_cmap='pyart_Theodore16' #'pyart_NWSRef' 'pyart_RRate11'
    
    display.plot_rhi('ZED_H',ax=ax1, sweep=sweep, vmin=-32, vmax=64.,colorbar_orient='vertical',
                     cmap=plot_cmap, axislabels=('', ylabel))
    display.plot_grid_lines(ax=ax1, col='k', ls=':')
    display.set_limits(xlim= [0,180], ylim=[0, 10])
    
    ax2=fig.add_subplot(312)
    #add_subplot(1 row1 of Figures, 2 columns of Figure, This is the second plot
    # top right)
    display.plot_rhi('ZDR',ax=ax2, sweep=sweep, vmin=-2, vmax=4.,colorbar_orient='vertical',title='ZDR',
                     cmap='pyart_RefDiff', axislabels=('',ylabel))
    display.plot_grid_lines(ax=ax2, col='k', ls=':')
    display.set_limits(xlim= [0,180], ylim=[0, 10])
    
    ax3=fig.add_subplot(313)
    display.plot_rhi('VEL_HV',ax=ax3, sweep=sweep,vmin=-1.0, vmax=1.0,colorbar_orient='vertical',title='VEL_HV',
                     cmap='pyart_RefDiff', axislabels=(xlabel, ylabel))
    display.plot_grid_lines(ax=ax3, col='k', ls=':')
    display.set_limits(xlim= [0,180], ylim=[0, 10])
    
    #Define File Name
    #Use normal String Manipulation 
    Figurename=figurelocation+'/'+nameoffile+'_pyart.png'
    
    fig.tight_layout() # Turn off to to plot faster
    
    
    #Use Matplotlip Save command
    #plt.show()
    plt.savefig(Figurename,dpi=200)
    plt.close(fig)
savefig()

HTML Menu Colour Change

I have written the entire computer project in HTML using Macromedia Dreamweaver 8. The first month of this degree I spent a great deal of time learning html and learned that having a style script ".css" meant I could have a much more effiency over the look of my website. One thing I was really keen on was to be able to change the colour and look of the menu icons, a simple task, yes? Wrong. At first, the hover would stay a changed colour when clicked so I had to reset the colour when visited back to its original colour. It took me a long time to figure out the hover and the link icon colours. The way that html works is you have your main website page background, then you have a padding between that and any menus or your main frame text, and you have padding between each menu and main frame text and that is different to a margin!

Script for Changing the Colour of the Menu Item when Hovering

Tools of the Trade

To comment in Python you use a '#'
Or to block comment a section a series of punctuation: """
Python Set-Up should always be as follows:

""""""
Author Name: Freya I. Lumb
Date of Last Edit: 21.05.2017
Title of Code: Volume of the Earth
Description of Code/Rationale: Example Script
Acknowledgements: None
""""""" 
Import Modules
import numpy as np
Define Global Variables 
R=6371
#Radius of the Earth in km
#Global variables are defined for the entire script, useful for constants such as the radius of the Earth, Avogadro's number, dielectric constant etc. 
#However these are only constant for the script you are in, so if you call a function from one script into another and it uses a global variable it does not work.
#Global variables should be used for constants and not if you are going to redefine them as this is bad coding and presents difficulties in your script.<br />
#There is a lot of debate on the internet how best to use global variables and functions, as long as they are used correctly it is fine. 
Define Functions
#Define separate functions
def function():
#do stuff
       V=(0.75)*np.pi*(R**3)
       return V
#Return functions and any variables required
print V
End Script

Panoply is the tool I use to look at NetCDF files but can be used for other datasets. It has versions for Mac, Windows and Linux.

Future Work

I have written some more code which isn't currently finished to save the new variables and timestamps to a file, this file will then be read by a prewritten code to plot the graphs; (ZED_H v. IWC), (IWC_aircraft v. IWC_radar), (ZDR v. IWC) etc.
There was some calibration issues for the flight I was originally working on but luckily a new flight was done on 17th May 2017, which I was a part of and we were able to capture the aircraft in some of the RHIs which will enable me to check my co-incident script. This has delayed issues. In addition I have been having trouble with metadata attached to my new variables causing them to be bigger arrays than intended, this is an issue I am still working on.
Once I have got the whole of the code working, I will be able to plot up and then apply statistics! Future work would be making the script usuable and faster, at the moment it is one ugly long piece of code, ideally I will split it up into functions but there is an issue with global variables from other files for this, so will need some though and collaboration to sort out. After that is done, the next thing is to make it usuable for other radar, aircraft users, this would be helped by splitting into functions, but also we can do it one of two ways, make the read in applicable for multiple files, or have separate read ins and then ask the question to the user what radar are you using? "X-Band" use the Xband read in file. Would you like ZED_H vs ZDR plotted? Yes, no etc.