Notebook Source Codes

Share your source codes

Source codes: 131

Plot Modis Airborne Simulator level 1 products using python 3

					
from pyhdf.SD import SD, SDC

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.cm as cm
import pprint

#----------------------------------------------------------------------------------------#

file_name = 'MASL1B_07915_05_20070707_1722_1737_V04.hdf'
file = SD(file_name, SDC.READ)

print(file.info())

#----------------------------------------------------------------------------------------#

datasets_dic = file.datasets()

for idx,sds in enumerate(datasets_dic.keys()):
    print( idx,sds )

#----------------------------------------------------------------------------------------#
# Get Central Wavelength

sds_obj = file.select('Central100%ResponseWavelength') # select sds

central_wavelength_list = sds_obj.get() # get sds data

print(central_wavelength_list)
print(len(central_wavelength_list))

for idx,central_wavelength in enumerate(central_wavelength_list):
	print(idx,central_wavelength)

#----------------------------------------------------------------------------------------#
# Get Calibrated data

sds_obj = file.select('CalibratedData') # select sds

print( sds_obj.info() )

data = sds_obj.get() # get sds data

print( data.shape )

attr = sds_obj.attributes()

pprint.pprint( sds_obj.attributes() )

scale_factor_list = attr['scale_factor']

fill_value = attr['_FillValue'] * 1.0

#print('attr',attr)
print('attr',scale_factor_list[0])

for idx,central_wavelength in enumerate(central_wavelength_list):
	
	cmap = [(0.90,0.90,0.90)] + [(cm.jet (i)) for i in range(1,256)] 
	cmap = mpl.colors.ListedColormap(cmap)
	
	value_max = np.max(data[:,idx,:]*scale_factor_list[idx]) * 1.0
	value_min = np.min(data[:,idx,:]*scale_factor_list[idx]) * 1.0
	
	print(value_max)

	bounds = [fill_value] + [i for i in np.linspace(value_min,value_max,10)] 
	norm = mpl.colors.BoundaryNorm(bounds, cmap.N)

	plt.imshow(data[:,idx,:]*scale_factor_list[idx], 
			   cmap=cmap, norm=norm, interpolation='none',
			   aspect='auto')

	cbar = plt.colorbar(boundaries=bounds, ticks=bounds)
	
	cbar.ax.set_yticklabels(['{:.2f}'.format(x) for x in bounds])
	
	plt.title('Central wavelength: '+str(round(central_wavelength,2))+'\n watts/meter2/steradian/micron')

	plt.savefig("mas_"+"%02d" % (idx+1)+".png", bbox_inches='tight', dpi=100)
	#plt.show()
	plt.close()

Python script to get a month of modis caliop cloudsat collocated indexes files

					
#!/usr/bin/env python

import os
import glob
import numpy as np
import calendar
import re

#----------------------------------------------------------------------------------------#
# user inputs

year = 2008
month = 7

machine = 'nasa' # (nasa,icare)

create_database_directory = True

if machine == 'nasa':
	path_to_myd03 = '/Volumes/Data/Myd03/'
	path_to_caliop = '/Volumes/Data/Caliop/01km_CLay/'
	path_to_cloudsat = '/Volumes/Data/Cloudsat/'
	path_to_results = '/Users/marchant/Desktop/'
	path_to_collocation_code = '/Users/marchant/Desktop/'

if machine == 'icare':
	path_to_myd03 = '/DATA/LIENS/MODIS/MYD03.006/'
	path_to_caliop = '/DATA/LIENS/CALIOP/01kmCLay.v4.10/'
	path_to_cloudsat = '/DATA/LIENS/CLOUDSAT/2B-GEOPROF/'
	path_to_results = '/home/benjamin.hg.marchant/data/modis_caliop_cloudsat_collocated_pixels_files/'
	path_to_collocation_code = '/home/benjamin.hg.marchant/src/'

#----------------------------------------------------------------------------------------#
# Create directory (Optional)

if create_database_directory:

	os.system( 'mkdir ' + path_to_results + str(year) )

	for day in np.arange(calendar.monthrange(year,month)[1]):
		d = './' + str(year) + '/' + str(year) + '_' + "%02d" % month + \
			 '_' + "%02d" % (day+1) + '/'
		os.system('mkdir ' + path_to_results + d)

#----------------------------------------------------------------------------------------#

for day in np.arange(calendar.monthrange(int(year),int(month))[1]):
        
	d = str(year) + '/' + str(year) + '_' + "%02d" % month + '_' + "%02d" % (day+1) + '/'
	
	path_to_myd03_day = path_to_myd03 + d
	path_to_caliop_day = path_to_caliop + d
	path_to_cloudsat_day = path_to_cloudsat + d

	list_of_myd03_files = glob.glob(path_to_myd03_day + '/*hdf' )
	list_of_caliop_files = glob.glob(path_to_caliop_day + '/*hdf' )
	list_of_cloudsat_files = glob.glob(path_to_cloudsat_day + '/*hdf' )        
        
	for modis_granule_name in list_of_myd03_files:        
		
		tmp = modis_granule_name.split('/')[-1]
		
		hm = tmp.split('.')[2]
		
		modis_granule_hour = int(hm[0:2])
		modis_granule_minute = int(hm[2:4])

		#--------------------------------------------------------------------------------#
		# Search caliop math file
		
		find_caliop_match_file = False
		minimum = 200
		for caliop_granule_name in list_of_caliop_files:
			
			caliop_re = re.compile(str(year) + '-' + \
			           "%02d" % month + '-' + "%02d" % (day+1) + "T(.*)-(.*)-(.*).hdf")
		
			obj = caliop_re.search(caliop_granule_name)
		
			caliop_granule_hour = int(obj.group(1))
			caliop_granule_minute = int(obj.group(2))
						
			difference = modis_granule_hour * 60 + modis_granule_minute - \
						(caliop_granule_hour * 60 + caliop_granule_minute)
	
			if (difference < minimum ) and ( difference > 0 ):
					minimum = difference
					tmp_02 = caliop_granule_name
					find_caliop_match_file = True

		#--------------------------------------------------------------------------------#
		# Search cloudsat math file

		find_cloudsat_match_file = False
		min_minute_diff = 99
		for cloudsat_file in list_of_cloudsat_files:
			
			cloudsat_file_name = cloudsat_file.split('/')[-1]
			cloudsat_granule_hour = int(cloudsat_file_name[7:9])
			cloudsat_granule_minute = int(cloudsat_file_name[9:11])
			cloudsat_granule_second = int(cloudsat_file_name[11:13] )
			
			minute_diff = modis_granule_hour * 60 + modis_granule_minute - \
					( cloudsat_granule_hour * 60 + cloudsat_granule_minute )
			
			if (minute_diff < min_minute_diff ) and ( minute_diff > 0 ):
					min_minute_diff = minute_diff
					cloudsat_match_file = cloudsat_file
					find_cloudsat_match_file = True

		#--------------------------------------------------------------------------------#
		
		print(find_caliop_match_file,find_cloudsat_match_file)
		
		if find_caliop_match_file and find_cloudsat_match_file:
			print(modis_granule_name,tmp_02,cloudsat_match_file)

			command = path_to_collocation_code + \
			          'create_modis_caliop_cloudsat_granule_collocated_pixel_file' + \
			          ' ' + tmp_02 + ' ' + \
					  modis_granule_name + ' ' + \
					  cloudsat_match_file
			os.system(command)
			command = 'mv modis_caliop_cloudsat_collocated_indexes_' + \
					  "%02d" % modis_granule_hour + \
					  "%02d" % modis_granule_minute + ' ' + \
					  path_to_results + d
			os.system(command)

			print( 'collocated file created')

	print('Day ' + str(day+1) + ' done')

How to download a month of CloudSat files from ICARE using python

					
#!/usr/bin/env python

import numpy as np
import os
import calendar

from ftplib import FTP

#----------------------------------------------------------------------------------------#
# user inputs

year = 2008
month = 7

path_to_cloudsat = '/DATA/LIENS/CLOUDSAT/2B-GEOPROF/'

icare_username = 'john.doe'
icare_password = '********'

#----------------------------------------------------------------------------------------#
# Create directory 

os.system( 'mkdir ' + str(year) )

for day in np.arange(calendar.monthrange(year,month)[1]):
	d = './' + str(year) + '/' + str(year) + '_' + "%02d" % month + \
		 '_' + "%02d" % (day+1) + '/'
	os.system('mkdir ' + d)
	
#----------------------------------------------------------------------------------------#
# Opening ICARE FTP connection

print( 'Opening FTP connection' )

ftp = FTP('ftp.icare.univ-lille1.fr')
ftp.login(icare_username,icare_password)

for day in np.arange(calendar.monthrange(int(year),int(month))[1]):

	directory = path_to_cloudsat + str(year) + '/' + str(year) + '_' + \
	            "%02d" % month + '_' + "%02d" % (day+1) + '/'	
	
	ftp.cwd(directory)

	yyy = []
	ftp.retrlines('NLST', yyy.append)
	#if len(yyy) != 288:
	#	print 'Probleme: donnees manquantes'
	for j in np.arange(len(yyy)):
		file_name = yyy[j]
		print(file_name)
		ftp.retrbinary('RETR ' + file_name, open(file_name, 'wb').write)
		command = 'mv ' + file_name + ' ' + './' + str(year) + '/' + str(year) + \
		'_' + "%02d" % month + '_' + "%02d" % (day+1) + '/' + file_name 
		os.system(command)

print( 'Closing FTP connection')
print( ftp.close())

Write a python script to find modis cloudsat granule match files

					
#!/usr/bin/env python

# code written with python 3 

import glob

#----------------------------------------------------------------------------------------#
# User inputs

#path_to_cloudsat_files = '/DATA/LIENS/CLOUDSAT/2B-GEOPROF/' # icare
path_to_cloudsat_files = '/Volumes/Data/Cloudsat/' # local

year = 2008
month = 1
day = 1

modis_granule_hour = 14
modis_granule_minute = 20

#----------------------------------------------------------------------------------------#
# Find match file

path_to_cloudsat_files = path_to_cloudsat_files + "%04d" % year + '/' + \
                         "%04d" % year + '_' + "%02d" % month + '_' + "%02d" % day

cloudsat_files_list = glob.glob(path_to_cloudsat_files + '/*hdf')

find_match_file = False
min_minute_diff = 99
for cloudsat_file in cloudsat_files_list:
	cloudsat_file = cloudsat_file.split('/')[-1]
	cloudsat_granule_hour = int(cloudsat_file[7:9])
	cloudsat_granule_minute = int(cloudsat_file[9:11])
	cloudsat_granule_second = int(cloudsat_file[11:13] )
	minute_diff = modis_granule_hour * 60 + modis_granule_minute - \
	        ( cloudsat_granule_hour * 60 + cloudsat_granule_minute )
	if (minute_diff < min_minute_diff ) and ( minute_diff > 0 ):
			min_minute_diff = minute_diff
			cloudsat_match_file = cloudsat_file
			find_match_file = True

if find_match_file: 
	print(cloudsat_match_file)
else:
	print('no file found !')

# returns here: 2008001133124_08930_CS_2B-GEOPROF_GRANULE_P_R04_E02.hdf

Plot CloudSat CPR cloud mask using python 3

					
#!/usr/bin/env python

from pyhdf.SD import SD, SDC 
from pyhdf.HDF import *
from pyhdf.VS import *
from mpl_toolkits.axes_grid1 import make_axes_locatable

import pprint
import matplotlib.pyplot as plt
import numpy as np
import matplotlib as mpl
import matplotlib.cm as cm

#----------------------------------------------------------------------------------------#

file_path = './cloudsat/plot_cpr_vertical_profile_from_2B-GEOPROF_py/'
file_name = '2008025010959_09272_CS_2B-GEOPROF_GRANULE_P_R04_E02.hdf'

#----------------------------------------------------------------------------------------#
# Force Aspect Ratio

def forceAspect(ax,aspect=1):
    im = ax.get_images()
    extent =  im[0].get_extent()
    ax.set_aspect(abs((extent[1]-extent[0])/(extent[3]-extent[2]))/aspect)

#----------------------------------------------------------------------------------------#
# Get CPR_Cloud_mask sds

file = SD(file_path+file_name, SDC.READ)

file_info = file.info()
print(file_info)  # number of sds and metadata

sds_obj = file.select('CPR_Cloud_mask') # select sds

data = sds_obj.get() 

sds_info = sds_obj.info()

print( 'sds info: CPR_Cloud_mask' )
print(data.shape)
print( sds_info )
print( sds_info[0], sds_info[1] )
print( 'sds attributes' )
pprint.pprint( sds_obj.attributes() )

#----------------------------------------------------------------------------------------#

f = plt.figure()
ax = f.add_subplot(111)

cmap = [(0.0,0.0,0.0)] + [(0.75,0.75,0.75)]
cmap = mpl.colors.ListedColormap(cmap)

bounds = [0,20,40]

norm = mpl.colors.BoundaryNorm(bounds, cmap.N)

img = plt.imshow(data.T, extent=[0,35000,-5,25], aspect='auto', 
                 cmap=cmap, norm=norm, interpolation='none')

cbar_bounds = [0,20,40]
cbar_ticks = [10,30]   
cbar_ticks_labels = ['', 'Cloudy']
               
cbar = plt.colorbar(img, fraction=0.01, cmap=cmap, norm=norm, boundaries=cbar_bounds, ticks=cbar_ticks)
cbar.ax.set_yticklabels(cbar_ticks_labels, fontsize=8)

plt.title('CloudSat CPR cloud mask', fontsize=8)

forceAspect(ax,aspect=3)

plt.xticks(fontsize=8)
plt.yticks(fontsize=8)

plt.savefig("cloudsat_cpr_cloud_mask.png", bbox_inches='tight', dpi=100)
plt.show()

Read latitudes and longitudes from a CloudSat 2B GEOPROF granule using fortran

					
program read_latitudes_and_longitudes_from_cloudsat_granule

! Compile:
!
! gfortran read_latitudes_and_longitudes_from_cloudsat_granule.f90 
!          -I/hdf-4.2.4-snap4/hdf4/include 
!          -L/hdf-4.2.4-snap4/hdf4/lib -lmfhdf -ldf -ljpeg -lz 
!          -o read_latitudes_and_longitudes_from_cloudsat_granule
!
! Run:
!
! cloudsat_file_name = 2008008133740_09032_CS_2B-GEOPROF_GRANULE_P_R04_E02.hdf
! read_latitudes_and_longitudes_from_cloudsat_granule cloudsat_file_name
      
implicit none

integer i

!----------------------------------------------------------------------------------------!

character*300 cloudsat_file_name 
character*300 vdata_name 

integer DFACC_READ, FULL_INTERLACE
parameter (DFACC_READ = 1, FULL_INTERLACE = 0)

integer hopen, hclose
integer vfstart, vsffnd, vsfatch, vsfsfld, vsfrd, vsfseek, vsfdtch, vfend, vsfinq

!----------------------------------------------------------------------------------------!

integer status
integer file_id, vdata_id, vdata_ref
integer n_records
character*80 fieldname_lists 
integer rec_num, num_of_records, rec_pos
integer interlace_mode, vdata_size


character*300 vdata_names_info 
character*80 fieldname_lists_info 
integer :: n_records_info, vdata_size_info

real, dimension(:), allocatable :: cloudsat_latitude
real, dimension(:), allocatable :: cloudsat_longitude

!----------------------------------------------------------------------------------------!

cloudsat_file_name = '2008008133740_09032_CS_2B-GEOPROF_GRANULE_P_R04_E02.hdf'

file_id = hopen(cloudsat_file_name, DFACC_READ, 0)

status = vfstart(file_id) 

!----------------------------------------------------------------------------------------!
! Latitude 

vdata_name = 'Latitude'
  
vdata_ref = vsffnd(file_id, vdata_name)

!write(6,*) vdata_ref

if (vdata_ref .eq. 0) stop
vdata_id = vsfatch(file_id, vdata_ref, 'r') 

!write(6,*) vdata_id

!----------------------------------------------------------------------------------------!

fieldname_lists = ' '
n_records = -1 
vdata_size = -1 

status = vsfinq(vdata_id, n_records, interlace_mode, fieldname_lists, &
& vdata_size, vdata_name)

write(6,*) 'n_records', n_records

allocate( cloudsat_latitude(n_records) )
allocate( cloudsat_longitude(n_records) )

!----------------------------------------------------------------------------------------!

n_records = vsfrd(vdata_id, cloudsat_latitude, n_records, FULL_INTERLACE)

do i = 1, 10 !n_records 
write(6,*) cloudsat_latitude(i)
end do

!----------------------------------------------------------------------------------------!
! Longitude 

write(6,*) '----------'

vdata_name = 'Longitude'
  
vdata_ref = vsffnd(file_id, vdata_name)

!write(6,*) vdata_ref

if (vdata_ref .eq. 0) stop
vdata_id = vsfatch(file_id, vdata_ref, 'r') 

!write(6,*) vdata_id

!----------------------------------------------------------------------------------------!

n_records = vsfrd(vdata_id, cloudsat_longitude, n_records, FULL_INTERLACE)

do i = 1, 10 !n_records 
write(6,*) cloudsat_longitude(i)
end do

!----------------------------------------------------------------------------------------!

status = vsfdtch(vdata_id)
status = vfend(file_id)
status = hclose(file_id)

end

Bootstrap formulaire vertical avec deux inputs sur une même ligne

					
<!DOCTYPE html>
<html lang="en">
<head>

  <title>Bootstrap Example</title>

  <meta charset="utf-8">

  <link rel="stylesheet" href="https://maxcdn.bootstrapcdn.com/bootstrap/3.3.7/css/bootstrap.min.css">

  <script src="https://ajax.googleapis.com/ajax/libs/jquery/3.2.1/jquery.min.js"></script>
  <script src="https://maxcdn.bootstrapcdn.com/bootstrap/3.3.7/js/bootstrap.min.js"></script>

</head>
<body>

<div class="container" >

<div class="row">

	<form action="/" method="post">

	<div class="col-md-4">
		<h2>Sign-up</h2> <br>
		<div class="row">
			<div class="form-group col-xs-6">
			<label for="e1">First name</label>
			<input type="text" class="form-control" id="e1" name="first_name">
			</div>
			<div class="form-group col-xs-6">
			<label for="e2">Last name</label>
			<input type="text" class="form-control" id="e2" name="last_name">		
			</div>		
		</div>
		<div class="form-group">
			<label for="e4">Username</label>
			<input type="text" class="form-control" id="e4" name="username">
		</div>	
		<div class="form-group">
			<label for="e3">Email Adress</label>
			<input type="email" class="form-control" id="e3" name="email">
		</div>	
		<div class="row">
			<div class="form-group col-xs-6">
				<label for="e5">Password</label>
				<input type="password" class="form-control" id="e5" name="password1">	
			</div>
			<div class="form-group col-xs-6">
				<label for="e6">Password confirmation</label>
				<input type="password" class="form-control" id="e6" name="password2">			
			</div>		
		</div>
		<br>
		<button type="submit" class="btn btn-primary btn-block">Create</button>
	</div>

	</form>

</div>


</div>



</body>
</html>

How to plot CloudSat trajectory using python and basemap

					
#!/usr/bin/env python

from mpl_toolkits.basemap import Basemap, cm
from pyhdf.SD import SD, SDC 
from pyhdf.HDF import *
from pyhdf.VS import *

import matplotlib.pyplot as plt
import numpy as np

#----------------------------------------------------------------------------------------#

file_path = './inputs/cloudsat/'
file_name = '2008025010959_09272_CS_2B-GEOPROF_GRANULE_P_R04_E02.hdf'

#----------------------------------------------------------------------------------------#
# Read Data

f = HDF(file_path+file_name, SDC.READ) 
vs = f.vstart() 

Latitude = vs.attach('Latitude')
Longitude = vs.attach('Longitude')

lat = Latitude[:]
long = Longitude[:]

Latitude.detach()
Longitude.detach() 
vs.end() 
f.close() 

#----------------------------------------------------------------------------------------#
# Plot Data

fig = plt.figure()

m = Basemap(projection='cyl',llcrnrlat=-90,urcrnrlat=90, 
            llcrnrlon=-180,urcrnrlon=180,resolution='c')

x,y = m(long,lat)

m.drawcoastlines()
m.drawparallels(np.arange(-90.,90.,30.))
m.drawmeridians(np.arange(-180.,180.,30.))

plt.title("Cloudsat Trajectory")

m.scatter(x,y,3,marker='o',color='b')

plt.savefig("cloudsat_trajectory.png", bbox_inches='tight', dpi=200)

plt.show()
plt.close()

Mie code written in python

					
#!/usr/bin/env python

# Author: Benjamin Marchant
# Ref: Hong Du's paper: "Mie-Scattering Calculation, Applied Optics, Vol. 43, Issue 9, pp. 1951-1956 (2004)

import matplotlib.pyplot as plt

import numpy as np
import math

#----------------------------------------------------------------------------------------#
# user inputs

ri = complex(1.33,-0.0) # refractive index

x = 5.0 # Mie Parameter (2*pi*radius/wavelength)

n = 100   

angular_resolution = 1.0

#----------------------------------------------------------------------------------------#
# Calculations of the Riccati-Bessel functions by upward recurrence.

PHI = np.zeros((n+1))

a = math.cos(x) 
PHI[0] = math.sin(x) 
b = PHI[0]
for i in range(1,n+1):
	PHI[i] = ( ( 2 * i - 1 ) / x ) * b  - a
	a = b
	b = PHI[i]

#----------------------------------------------------------------------------------------#

CHI = np.zeros((n+1))

a = math.cos(x)
CHI[0] = a
CHI[1] = math.cos(x) / x + math.sin(x)
b = CHI[1]
for i in range(2,n+1):
	CHI[i] = ( ( 2 * i - 1 ) / x ) * b  - a
	a = b
	b = CHI[i]

#----------------------------------------------------------------------------------------#

ZETA = np.zeros((n+1), dtype=complex)
for i in range(0,n+1):
	ZETA[i] = complex(PHI[i], CHI[i])

#----------------------------------------------------------------------------------------#
# Calculations of the complex functions rn(mx) by upward recurrence.

RN = np.zeros((n+1), dtype=complex)	

a = ri.real * x
b = ri.imag * x	

c1 = math.tan(a) - math.exp(-2.0*b) * math.tan(a)
c2 = 1.0 + math.exp(-2.0*b)
c3 = -1.0 + math.exp(-2.0*b)
c4 = math.tan(a) + math.exp(-2.0*b) * math.tan(a) 

a = c1 * c3 + c2 * c4
b = c2 * c3 - c1 * c4
c = c3 * c3 + c4 * c4

RN[0] = complex(a/c, b/c) 
for i in range(1,n+1):
	a = ( 2.0 * i - 1.0 ) / ( ri * x ) - RN[i-1]
	RN[i] = 1.0 / a	

#----------------------------------------------------------------------------------------#
# Calculations of the Mie coefficients an et bn.

MA = np.zeros((n), dtype=complex)
MB = np.zeros((n), dtype=complex)

for i in range(0,n):
	a = RN[i+1] / ri + (i + 1) * ( 1.0 - 1.0 / ( ri * ri ) ) / x 
	MA[i] = ( a * PHI[i+1] - PHI[i]) / ( a * ZETA[i+1] - ZETA[i] )
	MB[i] = ( RN[i+1] * ri * PHI[i+1] - PHI[i] ) / ( RN[i+1] * ri * ZETA[i+1] - ZETA[i] ) 

#----------------------------------------------------------------------------------------#

def angular_functions(x,angle,n,MA,MB):

	#------------------------------------------------------------------------------------#
	# Calculations of the Mie angular functions.

	PI = np.zeros((n))

	a = math.cos(angle) 
	b = 0.0
	c = 1.0

	PI[0] = 1.0
	for j in range(1,n):
		s = a * c
		t = s - b 
		PI[j] = s + t + t / j  
		b = c
		c = PI[j]

	#------------------------------------------------------------------------------------#
	
	TAU = np.zeros((n))
	
	TAU[0] = a
	for j in range(1,n):
		t = a * PI[j] - PI[j-1]
		TAU[j] = ( j + 1 ) * t - PI[j-1]
		
	#------------------------------------------------------------------------------------#
	# scattering amplitudes

	b = complex(0.0,0.0)
	c = complex(0.0,0.0)
	for i in range(0,n):
		a = ( 2.0 * ( i + 1 ) + 1.0 ) / ( ( i + 1 ) * ( i + 2 ) )
		b = b + a * ( MA[i] * PI[i] + MB[i] * TAU[i] )
		c = c + a * ( MB[i] * PI[i] + MA[i] * TAU[i] )
	sa1 = b
	sa2 = c

	return sa1, sa2

#----------------------------------------------------------------------------------------#
# phase matrix

angle_list = []
p11_list = []
p12_list = []
p33_list = []
p34_list = []

for angle in np.arange(0.0,180.0,angular_resolution):

	#print angle
	
	sa1,sa2 = angular_functions(x,math.radians(angle),n,MA,MB)
	
	p11 = ( sa1 * sa1.conjugate() + sa2 * sa2.conjugate() )
	p12 = ( sa2 * sa2.conjugate() - sa1 * sa1.conjugate() ) 
	p33 = ( sa2 * sa1.conjugate() + sa2 * sa1.conjugate() ) 
	p34 = complex(0,1) * ( sa1 * sa2.conjugate() - sa2 * sa1.conjugate() ) 

	angle_list.append(angle)
	p11_list.append(p11.real)
	p12_list.append(p12.real)
	p33_list.append(p33.real)
	p34_list.append(p34.real)		

#----------------------------------------------------------------------------------------#
# normalization

norm = 0 

for i in range(0,n):
	a = abs(MA[i]) * abs(MA[i]) + abs(MB[i]) * abs(MB[i])
	norm = norm +  ( 2.0 * ( i + 1.0 ) + 1.0 ) * a 

p11_list = p11_list / norm 
p12_list = p12_list / norm 
p33_list = p33_list / norm 
p34_list = p34_list / norm 

#----------------------------------------------------------------------------------------#
# extinction efficiency.

a = 0.0
for i in range(0,n):
	a = a + ( 2 * ( i + 1 ) + 1 ) * ( MA[i] + MB[i] ).real 

qe = ( 2.0 / ( x * x ) ) * a

#----------------------------------------------------------------------------------------#
# asymmetry factor.

a = 0.0
for i in range(0,n-1):
	c1 = MA[i] * MA[i+1].conjugate() + MB[i] * MB[i+1].conjugate()
	c2 = MA[i] * MB[i].conjugate()
	b = ( (i+1) * ( (i+1) + 2.0 ) ) / ( (i+1) + 1.0 )
	c = ( 2.0 * (i+1) + 1.0 ) / ( (i+1) * ( (i+1) + 1.0 ) )
	a = a + b * c1.real + c * c2.real

gfact = ( 2.0 / norm ) * a

#----------------------------------------------------------------------------------------#
# single scattering albedo.

albedo = ( 2.0 / ( x * x ) ) * norm / qe

#----------------------------------------------------------------------------------------#
# Data Visualization

fig = plt.figure()
ax = fig.add_subplot(111)

plt.plot(angle_list,p11_list,'k-')
plt.yscale('log')

plt.grid(True,which="both")

plt.xlabel(r"Scattering Angle $\Theta$ ($^\circ$)")
plt.ylabel(r"$P_{11}$")

plt.text(0.70,0.85,'x = '+str(x), fontsize=12, transform = ax.transAxes)
plt.text(0.70,0.80,'n=('+str(ri.real)+','+str(ri.imag)+')', \
                    fontsize=12, transform = ax.transAxes)     

plt.text(0.70,0.70,'qe = '+str(round(qe,2)), fontsize=12, transform = ax.transAxes)     
plt.text(0.70,0.65,'w0 = '+str(round(albedo,2)), fontsize=12, transform = ax.transAxes)     
plt.text(0.70,0.60,'g = '+str(round(gfact,2)), fontsize=12, transform = ax.transAxes)     

#plt.savefig('p11.png')
plt.show()

#----------------------------------------------------------------------------------------#
# plot p12

fig = plt.figure()
ax = fig.add_subplot(111)

plt.plot(angle_list,p12_list,'k-')
#plt.yscale('log')

plt.grid(True,which="both")

plt.xlabel(r"Scattering Angle $\Theta$ ($^\circ$)")
plt.ylabel(r"$P_{12}$")    

#plt.savefig('p12.png')
plt.show()

#----------------------------------------------------------------------------------------#
# plot p33

fig = plt.figure()
ax = fig.add_subplot(111)

plt.plot(angle_list,p33_list,'k-')
plt.yscale('log')

plt.grid(True,which="both")

plt.xlabel(r"Scattering Angle $\Theta$ ($^\circ$)")
plt.ylabel(r"$P_{33}$")    

#plt.savefig('p33.png')
plt.show()

#----------------------------------------------------------------------------------------#
# plot p34

fig = plt.figure()

ax = fig.add_subplot(111)

plt.plot(angle_list,p34_list,'k-')
#plt.yscale('log')

plt.grid(True,which="both")

plt.xlabel(r"Scattering Angle $\Theta$ ($^\circ$)")
plt.ylabel(r"$P_{34}$")    

#plt.savefig('p34.png')
plt.show()

How to plot polder/parasol L2 cloud cover using python and parasolascii ?

					
from mpl_toolkits.basemap import Basemap

import os
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.cm as cm

#----------------------------------------------------------------------------------------#

filename = 'P3L2TRGB071037KL'

command = 'parasolascii -c ' + filename + ' cc > ' + filename + '.ascii'

os.system(command)

#----------------------------------------------------------------------------------------#

data = np.loadtxt(filename+'.ascii')

#----------------------------------------------------------------------------------------#

polder_medium_resolution_grid = np.zeros((1080,2160))

polder_medium_resolution_grid[:,:] = -0.5

for idx,row in enumerate(data[:,0]):
	polder_medium_resolution_grid[data[idx,0]-1,data[idx,1]-1] = data[idx,2]

#----------------------------------------------------------------------------------------#

m = Basemap(projection='sinu', lon_0=0, resolution='c')

m.drawcoastlines()
#m.fillcontinents(color='coral',lake_color='aqua')
m.drawparallels(np.arange(-90.,120.,30.))
m.drawmeridians(np.arange(0.,420.,60.))
#m.drawmapboundary(fill_color='aqua')

cmap = [(0.90,0.90,0.90)] + [(cm.jet (i)) for i in range(1,256)] 
cmap = mpl.colors.ListedColormap(cmap)

bounds = [-1.0,0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0]
norm = mpl.colors.BoundaryNorm(bounds, cmap.N)

img = m.imshow(np.flipud(polder_medium_resolution_grid), 
               cmap=cmap, norm=norm, interpolation='none')

bounds = [0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0]

cb = m.colorbar(img,"bottom", size="5%", pad='4%', boundaries=bounds, 
                ticks=[0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0])

plt.title("PARASOL L2 Cloud Cover (P3L2TRGB071037KL) \n 2008-01-01 Sinusoidal Projection")
plt.savefig("polder_l2_cloud_cover_sinusoidal_projection_01.png",
            bbox_inches='tight', 
            dpi=400)
plt.show()

#----------------------------------------------------------------------------------------#

About


To add a source code to your notebook just click on the link below. Sharing scientific source code allows to understand and reproduce your results:

Add a code

Search


Source codes can be retrieved using the following form: