Notebook Source Codes

Share your source codes

Source codes: 131

Un réseau de neurones en python sur la reconnaissance optique de caractères

					
# python notebook for Make Your Own Neural Network
# code for a 3-layer neural network, and code for learning the MNIST dataset
# (c) Tariq Rashid, 2016
# license is GPLv2

import numpy
# scipy.special for the sigmoid function expit()
import scipy.special
# library for plotting arrays
import matplotlib.pyplot

# neural network class definition
class neuralNetwork:

	# initialise the neural network
	def __init__(self, inputnodes, hiddennodes, outputnodes, learningrate):
		# set number of nodes in each input, hidden, output layer
		self.inodes = inputnodes
		self.hnodes = hiddennodes
		self.onodes = outputnodes

		# link weight matrices, wih and who
		# weights inside the arrays are w_i_j, where link is from node i to node j in the next layer
		# w11 w21
		# w12 w22 etc 
		self.wih = numpy.random.normal(0.0, pow(self.inodes, -0.5), (self.hnodes, self.inodes))
		self.who = numpy.random.normal(0.0, pow(self.hnodes, -0.5), (self.onodes, self.hnodes))

		# learning rate
		self.lr = learningrate

		# activation function is the sigmoid function
		self.activation_function = lambda x: scipy.special.expit(x)

		pass

	# train the neural network
	def train(self, inputs_list, targets_list):
		# convert inputs list to 2d array
		inputs = numpy.array(inputs_list, ndmin=2).T
		targets = numpy.array(targets_list, ndmin=2).T

		# calculate signals into hidden layer
		hidden_inputs = numpy.dot(self.wih, inputs)
		# calculate the signals emerging from hidden layer
		hidden_outputs = self.activation_function(hidden_inputs)

		# calculate signals into final output layer
		final_inputs = numpy.dot(self.who, hidden_outputs)
		# calculate the signals emerging from final output layer
		final_outputs = self.activation_function(final_inputs)

		# output layer error is the (target - actual)
		output_errors = targets - final_outputs
		# hidden layer error is the output_errors, split by weights, recombined at hidden nodes
		hidden_errors = numpy.dot(self.who.T, output_errors) 

		# update the weights for the links between the hidden and output layers
		self.who += self.lr * numpy.dot((output_errors * final_outputs * (1.0 - final_outputs)), numpy.transpose(hidden_outputs))

		# update the weights for the links between the input and hidden layers
		self.wih += self.lr * numpy.dot((hidden_errors * hidden_outputs * (1.0 - hidden_outputs)), numpy.transpose(inputs))

		pass

	# query the neural network
	def query(self, inputs_list):		
		# convert inputs list to 2d array
		inputs = numpy.array(inputs_list, ndmin=2).T

		# calculate signals into hidden layer
		hidden_inputs = numpy.dot(self.wih, inputs)
		# calculate the signals emerging from hidden layer
		hidden_outputs = self.activation_function(hidden_inputs)

		# calculate signals into final output layer
		final_inputs = numpy.dot(self.who, hidden_outputs)
		# calculate the signals emerging from final output layer
		final_outputs = self.activation_function(final_inputs)

		return final_outputs


# number of input, hidden and output nodes
input_nodes = 784
hidden_nodes = 200
output_nodes = 10

# learning rate
learning_rate = 0.1

# create instance of neural network
n = neuralNetwork(input_nodes,hidden_nodes,output_nodes, learning_rate)

# load the mnist training data CSV file into a list
training_data_file = open("mnist_train_100.csv", 'r')
training_data_list = training_data_file.readlines()
training_data_file.close()

# train the neural network

# epochs is the number of times the training data set is used for training
epochs = 5

for e in range(epochs):
    # go through all records in the training data set
    for record in training_data_list:
        # split the record by the ',' commas
        all_values = record.split(',')
        # scale and shift the inputs
        inputs = (numpy.asfarray(all_values[1:]) / 255.0 * 0.99) + 0.01
        # create the target output values (all 0.01, except the desired label which is 0.99)
        targets = numpy.zeros(output_nodes) + 0.01
        # all_values[0] is the target label for this record
        targets[int(all_values[0])] = 0.99
        n.train(inputs, targets)
        pass
    pass

# load the mnist test data CSV file into a list
test_data_file = open("mnist_test_10.csv", 'r')
test_data_list = test_data_file.readlines()
test_data_file.close()

# test the neural network

# scorecard for how well the network performs, initially empty
scorecard = []

# go through all the records in the test data set
for record in test_data_list:
    # split the record by the ',' commas
    all_values = record.split(',')
    # correct answer is first value
    correct_label = int(all_values[0])
    # scale and shift the inputs
    inputs = (numpy.asfarray(all_values[1:]) / 255.0 * 0.99) + 0.01
    # query the network
    outputs = n.query(inputs)
    # the index of the highest value corresponds to the label
    label = numpy.argmax(outputs)
    # append correct or incorrect to list
    if (label == correct_label):
        # network's answer matches correct answer, add 1 to scorecard
        scorecard.append(1)
    else:
        # network's answer doesn't match correct answer, add 0 to scorecard
        scorecard.append(0)
        pass
    
    pass

# calculate the performance score, the fraction of correct answers
scorecard_array = numpy.asarray(scorecard)
print ("performance = ", scorecard_array.sum() / scorecard_array.size)

Tracer des images avec matplotlib provenant de la base de données MNIST

					
# Sources:
#
# Make Your Own Neural Network (by Tariq Rashid)
# https://github.com/makeyourownneuralnetwork/makeyourownneuralnetwork
# https://pjreddie.com/projects/mnist-in-csv/
# 
# Plot 4 records selected randomly from the mnist_train_100.csv file

from random import randint

import matplotlib.pyplot as plt
import numpy as np

training_data_file = open("mnist_train_100.csv", 'r')
training_data_list = training_data_file.readlines()
training_data_file.close()

for i in range(4):

	record = training_data_list[randint(0,99)]

	all_values = record.split(',')
	
	inputs_list = np.asfarray(all_values[1:])
	
	inputs = np.array(inputs_list, ndmin=1).T
	
	test = np.reshape(inputs, (28, 28))
	
	plt.imshow(test, cmap='Greys')
	
	plt.savefig('mnist_image_{:02d}'.format(i), bbox_inches='tight')
	plt.show()

	plt.close()

Plot modis caliop l2 dataset random sampling with a reservoir with matplotlib

					
#!/usr/bin/env python

from mpl_toolkits.basemap import Basemap, cm
from datetime import datetime
from pylab import *
import numpy as np
import matplotlib.pyplot as plt

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

lats, lons = np.loadtxt("modis_caliop_l2_dataset_unweighted_sampling.txt", unpack=True)

fig = plt.figure()

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

m.drawcoastlines()
#m.drawlsmask(land_color='white',ocean_color='white',lakes=True)
m.drawparallels(np.arange(-90.,90.,30.))
m.drawmeridians(np.arange(-180.,180.,30.))
#m.drawmapboundary(fill_color='aqua')

x, y = m(lons,lats)

m.scatter(x,y,6,marker='o',color='red')

plt.title('Modis caliop l2 dataset \n unweighted_sampling')
#plt.title('Weighted random sampling with a reservoir \n (number of learning data = 1000)')

plt.savefig('modis_caliop_l2_dataset_unweighted_sampling.png',  bbox_inches='tight')
#plt.show()
plt.close()

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

lats, lons = np.loadtxt("modis_caliop_l2_dataset_weighted_sampling.txt", unpack=True)

fig = plt.figure()

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

m.drawcoastlines()
#m.drawlsmask(land_color='white',ocean_color='white',lakes=True)
m.drawparallels(np.arange(-90.,90.,30.))
m.drawmeridians(np.arange(-180.,180.,30.))
#m.drawmapboundary(fill_color='aqua')

x, y = m(lons,lats)

m.scatter(x,y,6,marker='o',color='red')

plt.title('Modis caliop l2 dataset \n weighted_sampling')
#plt.title('Weighted random sampling with a reservoir \n (number of learning data = 1000)')

plt.savefig('modis_caliop_l2_dataset_weighted_sampling.png',  bbox_inches='tight')
#plt.show()
plt.close()

Create modis caliop L2 collocated dataset weighted random sampling with fortran

					
program modis_caliop_baseline

use read_write_hdf_file_variables 
use get_modis_c6_myd03_data
use get_modis_c6_myd06_data
use get_caliop_v4_01kmCLay_data
use get_caliop_v4_05kmCLay_data

implicit none 

integer, parameter :: pr = selected_real_kind(15,3)

integer :: i,j,k,l,m,n
integer :: mod_file
integer :: year,month,day,hour,minute,julian_day
character*300 :: path_calmod
character*100 :: calmod_index,calmod_file_name
character*400 :: input_colocated_file_index,input_colocated_file_name
integer :: ierror, iazm, iszen, iview
integer :: string_len
integer :: ierr ! successful (if == 0) ou failed (si /= 0)
integer :: nb_colocated_data
logical :: found_colocated_data, missing_file
integer, dimension(:,:), allocatable :: modis_index
integer, dimension(:), allocatable :: caliop_index
integer :: AllocateStatus, DeAllocateStatus
character*300 :: input_MYD03, input_MYD35
character*300 :: input_MYD06_c5, input_MYD06_c6, input_MYD06_baseline_c6, input_MYD06_research_c6
character*300 :: input_CAL_01km, input_CAL_05km, input_MYD021KM
character*16 :: file_common_path
integer, dimension(3) :: ShapeArray

integer :: s1,s2,s3,s4,s5,s6,s7,s8

!----- new variables -----!

integer :: nb_total_collocated_data 

real(pr) :: grid_resolution
integer :: grid_bin, idx, jdx, idx_max, jdx_max
integer, dimension(:,:), allocatable :: counts_aggregation

integer :: reservoir_sampling_size, reservoir_sampling_idx

real(pr), dimension(:,:), allocatable :: reservoir_sampling

integer :: lat_idx
integer, dimension(180) :: lat_weight_counts

integer*4 timeArray(3) ! Holds the hour, minute, and second

integer :: reservoir_weigthed_sampling_size
integer :: reservoir_weigthed_sampling_idx
integer :: reservoir_weigthed_sampling_min_k
real(pr), dimension(:,:), allocatable :: reservoir_weigthed_sampling

integer, dimension(180) :: lat_weight

integer :: replace_idx

integer, dimension(1) :: min_idx

real(pr) :: current_k

real(pr) :: ui

!----- new variables (end) -----!

integer, dimension(12) :: dpm_nyear, dpm_byear
data dpm_nyear /31,28,31,30,31,30,31,31,30,31,30,31/
data dpm_byear /31,29,31,30,31,30,31,31,30,31,30,31/

!----------------------------------------------------------------------------------------!
! inputs

year = 2008
month = 07

julian_day = 0
do i = 1, month - 1
julian_day = julian_day + dpm_byear(i)
end do
julian_day = julian_day + day

path_calmod = '/Volumes/Data/Colocated_data/Caliop_modis/01kmCLay/    /    _  _  /'
calmod_index = 'calmod_index_    '
calmod_file_name = 'calmod_file_name_    '

!----------------------------------------------------------------------------------------!
! set-up paths

file_common_path = '    /    _  _  /'
write(file_common_path(1:4),'(i4.4)') year
write(file_common_path(6:9),'(i4.4)') year
write(path_calmod(52:55),'(i4.4)') year
write(path_calmod(57:60),'(i4.4)') year
write(file_common_path(11:12),'(i2.2)') month
write(path_calmod(62:63),'(i2.2)') month

!----------------------------------------------------------------------------------------!
! Initialization

nb_total_collocated_data = 0 

grid_resolution = 1.0

idx_max = int( 180.0 / grid_resolution )
jdx_max = int( 360.0 / grid_resolution ) 

allocate( counts_aggregation( idx_max , jdx_max ) )

counts_aggregation = 0

reservoir_sampling_size = 1000

allocate( reservoir_sampling( reservoir_sampling_size , 2 ) )

reservoir_sampling = 0

lat_weight_counts = 0

call itime(timeArray) ! Generate sequence of random number 

i = rand( timeArray(1) + timeArray(2) + timeArray(3) )

reservoir_weigthed_sampling_size = 1000

allocate( reservoir_weigthed_sampling( reservoir_weigthed_sampling_size , 3 ) )

reservoir_weigthed_sampling_min_k = 9999

open(1,file='./inputs/modis_c6_caliop_weighted_random_sampling_with_a_reservoir/lat_weight.txt')
do i = 1 , 180
read(1,*) lat_weight(i)
end do
close(1)

!----------------------------------------------------------------------------------------!
! start loop

do day = 1, 8 !31 

write(file_common_path(14:15),'(i2.2)') day
write(path_calmod(65:66),'(i2.2)') day

hour = 0
minute = 0

do mod_file = 1, 288

write(calmod_index(14:15),'(i2.2)') hour
write(calmod_index(16:17),'(i2.2)') minute
write(calmod_file_name(18:19),'(i2.2)') hour
write(calmod_file_name(20:21),'(i2.2)') minute

string_len = len( trim(path_calmod) )
input_colocated_file_index = path_calmod(1:string_len)//calmod_index
input_colocated_file_name = path_calmod(1:string_len)//calmod_file_name

found_colocated_data = .TRUE.
open(1,file=input_colocated_file_index,form="formatted",iostat=ierr,status="old")
if (ierr/=0) then
!write(*,*) "Impossible d ouvrir le fichier ",trim(input_colocated_file_index)
!stop
found_colocated_data = .FALSE.
endif
nb_colocated_data = 0
do while (ierr==0)
read(1,*,iostat=ierr)
if (ierr==0) nb_colocated_data = nb_colocated_data + 1
enddo
close(1)

missing_file = .FALSE.
open(1,file=input_colocated_file_name,form="formatted",iostat=ierr,status="old")
if (ierr/=0) then
!write(*,*) "Impossible d ouvrir le fichier ",trim(input_colocated_file_name)
!stop
missing_file = .TRUE.
endif
close(1)

if( found_colocated_data .and. .not. missing_file )then
nb_colocated_data = nb_colocated_data - 2 
! First two lines are the name of caliop and modis files used to do the colocation
if( nb_colocated_data > 0 )then

allocate( modis_index(nb_colocated_data,2), STAT = AllocateStatus )
IF (AllocateStatus /= 0) STOP "*** Not enough memory modis_index***"
allocate( caliop_index(nb_colocated_data), STAT = AllocateStatus )
IF (AllocateStatus /= 0) STOP "*** Not enough memory caliop_index***"

open(1,file=input_colocated_file_index, status="old")
read(1,'(a)') 
read(1,'(a)') 
do i = 1, nb_colocated_data 
read(1,*) modis_index(i,1), modis_index(i,2), caliop_index(i)
end do
close(1)

open(1,file=input_colocated_file_name, status="old")
read(1,'(a)') input_MYD03
read(1,'(a)') input_MYD35
read(1,'(a)') input_MYD06_c5
read(1,'(a)') input_MYD06_c6
read(1,'(a)') input_MYD06_baseline_c6
read(1,'(a)') input_MYD06_research_c6
read(1,'(a)') input_CAL_01km
read(1,'(a)') input_CAL_05km
read(1,'(a)') input_MYD021KM
close(1)

input_MYD03 = '/Volumes/Data/Myd03/'//file_common_path//input_MYD03
!input_MYD35 = '/Volumes/DATA_03/modis/Myd35/'//file_common_path//input_MYD35
!input_MYD06_c5 = '/Volumes/DATA_03/modis/Myd06_c5/'//file_common_path//input_MYD06_c5
input_MYD06_c6 = '/Volumes/DATA_03/modis/Myd06_c6/'//file_common_path//input_MYD06_c6
!input_MYD06_baseline_c6 = '/Volumes/DATA_03/modis/Myd06_baseline_c6/'//file_common_path//input_MYD06_baseline_c6
!input_MYD06_research_c6 = '/Volumes/DATA_03/modis/Myd06_research_c6/'//file_common_path//input_MYD06_research_c6
input_CAL_01km = '/Volumes/Data/Caliop/v4/01km_CLay/'//file_common_path//input_CAL_01km
input_CAL_05km = '/Volumes/Data/Caliop/v4/05km_CLay/'//file_common_path//input_CAL_05km
!input_MYD021KM = '/Volumes/Data/MYD021km/'//file_common_path//input_MYD021KM

CALL read_modis_c6_myd03(input_MYD03)
CALL read_modis_c6_myd06(input_MYD06_c6)
CALL read_caliop_v4_01kmCLay(input_CAL_01km)
CALL read_caliop_v4_05kmCLay(input_CAL_05km)

!----------------------------------------------------------------------------------------!
! start analysis

ShapeArray = shape(modis_c6_myd06_Cloud_Mask_1km)

do i = 1, nb_colocated_data 

j = caliop_index(i)
k = modis_index(i,1)
l = modis_index(i,2)

if( modis_c6_myd03_solarzenith(k,l) * 0.01 < 81.36 .and. modis_c6_myd03_solarzenith(k,l) >= 0.0 )then

nb_total_collocated_data = nb_total_collocated_data + 1

if( ibits(modis_c6_myd06_Cloud_Mask_1km(1,k,l),0,1) == 1 .and. & ! Cloud Mask determined 
&   ibits(modis_c6_myd06_Cloud_Mask_1km(1,k,l),1,2) < 2 )then  ! Cloudy or Probably Cloudy Pixel  
s1 = 1
else
s1 = 2
end if

!write(6,*) s1

idx = int( (modis_c6_myd03_latitude(k,l) + 90.0) / grid_resolution ) + 1
jdx = int( (modis_c6_myd03_longitude(k,l) + 180.0) / grid_resolution ) + 1

if( idx > idx_max) idx = idx_max
if( jdx > jdx_max) jdx = jdx_max

!write(6,*) modis_c6_myd03_latitude(k,l), modis_c6_myd03_longitude(k,l), idx, jdx 

counts_aggregation(idx,jdx) = counts_aggregation(idx,jdx) + 1

lat_idx = int( (modis_c6_myd03_latitude(k,l) + 90.0) ) + 1
if( lat_idx > 180 ) lat_idx = 180

lat_weight_counts(lat_idx) = lat_weight_counts(lat_idx) + 1

!----------------------------------------------------------------------------------------!
! Reservoir Sampling

if( nb_total_collocated_data <= reservoir_sampling_size )then

reservoir_sampling_idx = nb_total_collocated_data

reservoir_sampling(reservoir_sampling_idx,1) = modis_c6_myd03_latitude(k,l)
reservoir_sampling(reservoir_sampling_idx,2) = modis_c6_myd03_longitude(k,l)

else

reservoir_sampling_idx = int( rand(0) * ( nb_total_collocated_data - 1 ) ) + 1

if( reservoir_sampling_idx <= reservoir_sampling_size )then

reservoir_sampling(reservoir_sampling_idx,1) = modis_c6_myd03_latitude(k,l)
reservoir_sampling(reservoir_sampling_idx,2) = modis_c6_myd03_longitude(k,l)

end if

end if

!----------------------------------------------------------------------------------------!
! Weighted random sampling with a reservoir

if( nb_total_collocated_data <= reservoir_weigthed_sampling_size )then

	reservoir_weigthed_sampling_idx = nb_total_collocated_data

	ui = rand(0)

	lat_idx = int( (modis_c6_myd03_latitude(k,l) + 90.0) ) + 1
	if( lat_idx > 180 ) lat_idx = 180

	reservoir_weigthed_sampling(reservoir_weigthed_sampling_idx,1) = modis_c6_myd03_latitude(k,l)
	reservoir_weigthed_sampling(reservoir_weigthed_sampling_idx,2) = modis_c6_myd03_longitude(k,l)
	reservoir_weigthed_sampling(reservoir_weigthed_sampling_idx,3) = ui**( lat_weight(lat_idx) )

	if( reservoir_weigthed_sampling(reservoir_weigthed_sampling_idx,3) < reservoir_weigthed_sampling_min_k ) &
	& reservoir_weigthed_sampling_min_k = reservoir_weigthed_sampling(reservoir_weigthed_sampling_idx,3)

else

	ui = rand(0)

	lat_idx = int( (modis_c6_myd03_latitude(k,l) + 90.0) ) + 1
	if( lat_idx > 180 ) lat_idx = 180

	current_k = ui**( lat_weight(lat_idx) )
	
	if( current_k > reservoir_weigthed_sampling_min_k )then
	
	min_idx = minloc(reservoir_weigthed_sampling(:,3))
	
	replace_idx = min_idx(1)
	
	reservoir_weigthed_sampling(replace_idx,1) = modis_c6_myd03_latitude(k,l)
	reservoir_weigthed_sampling(replace_idx,2) = modis_c6_myd03_longitude(k,l)
	reservoir_weigthed_sampling(replace_idx,3) = current_k
	
	min_idx = minloc(reservoir_weigthed_sampling(:,3))
	replace_idx = min_idx(1)
	reservoir_weigthed_sampling_min_k = reservoir_weigthed_sampling(replace_idx,3)
	
	end if

end if

end if !----- solar zenith

end do !----- nb_colocated_data granule loop

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

deallocate(modis_c6_myd03_latitude, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd03_latitude) ***"
deallocate(modis_c6_myd03_longitude, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd03_longitude) ***"
deallocate(modis_c6_myd03_sensorzenith, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd03_sensorzenith) ***"
deallocate(modis_c6_myd03_sensorazimuth, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd03_sensorazimuth) ***"
deallocate(modis_c6_myd03_solarzenith, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd03_solarzenith) ***"
deallocate(modis_c6_myd03_solarazimuth, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd03_solarazimuth) ***"

deallocate(modis_c6_myd06_Cloud_Phase_Infrared_1km, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Phase_Infrared_1km) baseline***"
deallocate(modis_c6_myd06_Cloud_Phase_Optical_Properties, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Phase_Optical_Properties) baseline***"
deallocate(modis_c6_myd06_Cloud_Mask_1km, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Mask_1km) baseline***"
deallocate(modis_c6_myd06_Quality_Assurance_1km, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Quality_Assurance_1km) baseline***"
deallocate(modis_c6_myd06_Cloud_Optical_Thickness, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Optical_Thickness) baseline***"
deallocate(modis_c6_myd06_Cloud_Optical_Thickness_PCL, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Optical_Thickness_PCL) baseline***"
deallocate(modis_c6_myd06_cloud_top_temperature_1km, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_cloud_top_temperature_1km) baseline***"
deallocate(modis_c6_myd06_Cloud_Top_Temperature, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Top_Temperature) baseline***"
deallocate(modis_c6_myd06_Cloud_Multi_Layer_Flag, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Multi_Layer_Flag) baseline***"
deallocate(modis_c6_myd06_Cloud_Effective_Radius, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Effective_Radius) baseline***"
deallocate(modis_c6_myd06_Cloud_Effective_Radius_16, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Effective_Radius_16) baseline***"
deallocate(modis_c6_myd06_Cloud_Effective_Radius_37, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Effective_Radius_37) baseline***"
deallocate(modis_c6_myd06_Cloud_Top_Height_1km, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Top_Height_1km) baseline***"

deallocate(caliop_v4_01kmCLay_latitude, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_01kmCLay_latitude)***"
deallocate(caliop_v4_01kmCLay_longitude, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_01kmCLay_longitude)***"
deallocate(caliop_v4_01kmCLay_Opacity_Flag, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_01kmCLay_Opacity_Flag)***"
deallocate(caliop_v4_01kmCLay_Number_Layers_Found, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_01kmCLay_Number_Layers_Found)***"
deallocate(caliop_v4_01kmCLay_Feature_Classification_Flags, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_01kmCLay_Feature_Classification_Flags)***"
deallocate(caliop_v4_01kmCLay_Layer_Top_Altitude, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_01kmCLay_Layer_Top_Altitude)***"
deallocate(caliop_v4_01kmCLay_Layer_Base_Altitude, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_01kmCLay_Layer_Base_Altitude)***"

deallocate(caliop_v4_05kmCLay_Number_Layers_Found, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_05kmCLay_Number_Layers_Found)***"
deallocate(caliop_v4_05kmCLay_Feature_Classification_Flags, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_05kmCLay_Feature_Classification_Flags)***"
deallocate(caliop_v4_05kmCLay_Feature_Optical_Depth_532, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_05kmCLay_Feature_Optical_Depth_532)***"
deallocate(caliop_v4_05kmCLay_Layer_Top_Altitude, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_05kmCLay_Layer_Top_Altitude)***"
deallocate(caliop_v4_05kmCLay_Layer_Base_Altitude, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_05kmCLay_Layer_Base_Altitude)***"

deallocate(modis_index, STAT = DeAllocateStatus)
deallocate(caliop_index, STAT = DeAllocateStatus)

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

end if !----- if nb_colocated_data > 0
end if !----- if found colocated file

write(6,*) 'hour, minute: ', hour, minute
minute = minute + 5
if ( minute > 55 ) then
hour = hour + 1
minute = 0
end if
end do !----- Granule loop
write(6,*) 'day: ', day
end do !----- day loop

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

write(6,*) 'nb_total_collocated_data', nb_total_collocated_data

open(1,file='./outputs/modis_c6_caliop_weighted_random_sampling_with_a_reservoir/counts_aggregation.txt')
do idx = 1, idx_max
write(1,*) (counts_aggregation(idx,jdx), jdx = 1, jdx_max)
end do
close(1)

!open(1,file='./outputs/modis_c6_caliop_weighted_random_sampling_with_a_reservoir/lat_weight.txt')
!do idx = 1, 180
!write(1,*) lat_weight_counts(idx)
!end do
!close(1)

open(1,file='./outputs/modis_c6_caliop_weighted_random_sampling_with_a_reservoir/resevoir_sampling.txt')
do idx = 1, reservoir_sampling_size
write(1,*) reservoir_sampling(idx,1), reservoir_sampling(idx,2)
end do
close(1)

open(1,file='./outputs/modis_c6_caliop_weighted_random_sampling_with_a_reservoir/reservoir_weigthed_sampling.txt')
do idx = 1, reservoir_weigthed_sampling_size
write(1,*) reservoir_weigthed_sampling(idx,1), reservoir_weigthed_sampling(idx,2)
end do
close(1)

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

end program modis_caliop_baseline

Create a modis caliop cloudsat granule collocated pixel file using fortran

					
program colocated_modis_caliop

! Compile:
!
! gfortran create_modis_caliop_cloudsat_granule_collocated_pixel_file.f90 
!          -I/hdf-4.2.4-snap4/hdf4/include 
!          -L/hdf-4.2.4-snap4/hdf4/lib -lmfhdf -ldf -ljpeg -lz 
!          -o create_modis_caliop_cloudsat_granule_collocated_pixel_file
!
! Inputs (example): 
!
! myd03_file_name = MYD03.A2008008.1420.006.2012066144951.hdf
! caliop_file_name = CAL_LID_L2_01kmCLay-ValStage1-V3-01.2008-01-08T13-55-27ZD.hdf
! cloudsat_file_name = 2008008133740_09032_CS_2B-GEOPROF_GRANULE_P_R04_E02.hdf
!
! create_modis_caliop_cloudsat_granule_collocated_pixel_file CAL_LID_L2_01kmCLay-ValStage1-V3-01.2008-01-08T13-55-27ZD.hdf MYD03.A2008008.1420.006.2012066144951.hdf 2008008133740_09032_CS_2B-GEOPROF_GRANULE_P_R04_E02.hdf  2008008133740_09032_CS_2B-GEOPROF_GRANULE_P_R04_E02.hdf 
       
implicit none

integer :: i,j,k,l,n

integer, parameter :: pr = selected_real_kind(15,3)

integer :: ierr ! Read File Status: Success (if == 0) or Fail (if /= 0) 
integer :: uin=42 ! logical unita associated to the file 
character*300 caliop_file_name, myd03_file_name, cloudsat_file_name 
character*300 fn
character*2 day, heure, minute
integer :: offset

integer DFACC_READ , DFNT_INT32
parameter (DFACC_READ = 1, DFNT_INT32 = 24)
integer sfstart, sfselect, sfrdata, sfendacc, sfend 
integer sffattr , sfrattr, var_id, attr_id, err_code
real*4 :: scale_factors(20), offsets(20)
integer sd_id, sds_id, sds_index
integer status, n_attrs
integer sffinfo , sfginfo
integer n_datasets, n_file_attrs , index
integer rank, data_type
integer start(3), edges(3), stride(3)

integer X_LENGTH_01 , Y_LENGTH_01 
integer X_LENGTH_01_max , Y_LENGTH_01_max 
parameter ( X_LENGTH_01_max = 1354 , Y_LENGTH_01_max = 2040 ) 
integer :: nb_sds_index 
parameter ( nb_sds_index = 2 )
integer tab_sds_index(nb_sds_index)
integer MAX_NC_NAME , MAX_VAR_DIMS
parameter ( MAX_NC_NAME = 256 , MAX_VAR_DIMS = 32 )
integer dim_sizes( MAX_VAR_DIMS )
character*( MAX_NC_NAME ) name 

!----------------------------------------------------------------------------------------!
! Caliop Data 

real(kind=4), dimension(:,:), allocatable :: caliop_latitude ! data type 5 
real(kind=4), dimension(:,:), allocatable :: caliop_longitude ! data type 5 

real(kind=4), dimension(:,:), allocatable :: Geographic_Coordinate_v3 ! data type 5 

integer :: nb_sds_index_v3 
parameter ( nb_sds_index_v3 = 2 )
integer tab_sds_index_v3( nb_sds_index_v3 )      
data tab_sds_index_v3 /1,2/
integer :: caliop_dim_sizes

integer :: nb_caliop_data_screen
real(pr) :: min_caliop_longitude_screen, max_caliop_longitude_screen
real(pr) :: min_caliop_latitude_screen, max_caliop_latitude_screen
integer, dimension(:), allocatable :: caliop_index_screen
real(kind=4), dimension(:), allocatable :: caliop_latitude_screen ! data type 5
real(kind=4), dimension(:), allocatable :: caliop_longitude_screen ! data type 5 !

!----------------------------------------------------------------------------------------!
! MYD03 Data

real(kind=4), dimension(:,:), allocatable :: mod03_latitude ! data type 5
real(kind=4), dimension(:,:), allocatable :: mod03_longitude ! data type 5

real(kind=4), dimension(:,:,:), allocatable :: Geographic_Coordinate_mod03 ! data type 5

integer :: nb_sds_index_mod03
parameter ( nb_sds_index_mod03 = 2 )
integer tab_sds_index_mod03( nb_sds_index_mod03 )
data tab_sds_index_mod03 /0,1/
integer :: Cell_Along_Swath_1km , Cell_Across_Swath_1km

integer :: nb_modis_data_screen
real(pr) :: min_modis_longitude , max_modis_longitude
real(pr) :: min_modis_latitude , max_modis_latitude
real(pr) :: min_modis_longitude_screen, max_modis_longitude_screen
real(pr) :: min_modis_latitude_screen, max_modis_latitude_screen
integer, dimension(:,:), allocatable :: modis_index_screen
real(kind=4), dimension(:), allocatable :: modis_latitude_screen ! data type 5
real(kind=4), dimension(:), allocatable :: modis_longitude_screen ! data type 5

!----------------------------------------------------------------------------------------!
! CloudSat Data

integer :: nb_cloudsat_pixels

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

character*300 vdata_name

integer FULL_INTERLACE
parameter (FULL_INTERLACE = 0)

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

integer file_id, vdata_id, vdata_ref
integer n_records
character*80 fieldname_lists 
integer interlace_mode, vdata_size

!----------------------------------------------------------------------------------------!
! Collocation Variables

real(pr) :: distance
real(pr) :: distance_pr

real(pr) :: minimum_distance
real(pr) :: distance_min

real(pr) :: earth_radius = 6371.0d0
real(pr) :: diff_lat, diff_long 
real(pr) :: phim
real(pr) :: K1, K2
real(pr) :: tmp_01, tmp_02, tmp_03
integer :: min_j, min_k

logical :: collocation

real(pr) :: cdr
real(pr) :: pi = dacos(-1.0d0)

integer :: i_pr, j_pr

integer :: i_min, j_min

logical :: find_min_distance
logical :: find_first_collocated_pixel
logical :: find_modis_pixel

integer :: tmp_n

integer :: nb_collocated_pixels

integer, dimension(:,:), allocatable :: collocated_pixels 

cdr = pi / 180.0d0

minimum_distance = 1.0d0

!----------------------------------------------------------------------------------------!
! Read Caliop (Longitude & Latitude) 

call getarg(1,caliop_file_name)

sd_id = sfstart(caliop_file_name, DFACC_READ)
status = sffinfo( sd_id , n_datasets , n_file_attrs )

do n = 1 , nb_sds_index_v3

	sds_index = tab_sds_index_v3(n)
	sds_id = sfselect(sd_id, sds_index)
	status = sfginfo( sds_id , name , rank , dim_sizes , data_type , n_attrs )
	
	do i = 1 , rank
	start(i) = 0
	edges(i) = dim_sizes(i)
	stride(i) = 1
	end do

	!---------- caliop_latitude	----------!

	if( sds_index == 1 )then
	allocate( caliop_latitude(1,dim_sizes(2)) )
	status = sfrdata( sds_id, start, stride, edges, caliop_latitude )
	end if

	!---------- caliop_longitude ----------!

	if( sds_index == 2 )then
	allocate( caliop_longitude(1,dim_sizes(2)) )
	status = sfrdata( sds_id, start, stride, edges, caliop_longitude )
	end if

end do

caliop_dim_sizes = dim_sizes(2)

status = sfendacc(sds_id)
status = sfend(sd_id)

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

allocate( collocated_pixels(caliop_dim_sizes,4) )

collocated_pixels = 0

!----------------------------------------------------------------------------------------!
! Read MYD03 (Longitude & Latitude) 

call getarg(2,myd03_file_name)

sd_id = sfstart(myd03_file_name, DFACC_READ)
status = sffinfo( sd_id , n_datasets , n_file_attrs )

do n = 1 , nb_sds_index_mod03

	sds_index = tab_sds_index_mod03(n)
	sds_id = sfselect(sd_id, sds_index)
	status = sfginfo( sds_id , name , rank , dim_sizes , data_type , n_attrs )

	do i = 1 , rank
	start(i) = 0
	edges(i) = dim_sizes(i)
	stride(i) = 1
	end do

	!---------- myd03_longitude ----------!

	if( sds_index == 0 )then
	k = 1
	allocate( mod03_latitude(dim_sizes(1),dim_sizes(2)))
	allocate( Geographic_Coordinate_mod03(dim_sizes(1),dim_sizes(2),2) )
	status = sfrdata( sds_id, start, stride, edges, mod03_latitude )
	end if

	!---------- myd03_longitude ----------!
	
	if( sds_index == 1 )then
	k = 2
	allocate( mod03_longitude(dim_sizes(1),dim_sizes(2)))
	status = sfrdata( sds_id, start, stride, edges, mod03_longitude )
	end if


end do

status = sfendacc(sds_id)
status = sfend(sd_id)

Cell_Across_Swath_1km = dim_sizes(1)
Cell_Along_Swath_1km = dim_sizes(2) 

min_modis_longitude = MINVAL( mod03_longitude )
max_modis_longitude = MAXVAL( mod03_longitude )
min_modis_latitude = MINVAL( mod03_latitude )
max_modis_latitude = MAXVAL( mod03_latitude )

!----------------------------------------------------------------------------------------!
! Read CLOUDSAT (Longitude & Latitude) 

call getarg(3,cloudsat_file_name)

file_id = hopen(cloudsat_file_name, DFACC_READ, 0)

status = vfstart(file_id) 


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

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

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

nb_cloudsat_pixels = n_records

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

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

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)

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

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

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

!write(6,*) 'min_modis_latitude, max_modis_latitude', min_modis_latitude, max_modis_latitude
!write(6,*) 'min_modis_longitude, max_modis_longitude', min_modis_longitude, max_modis_longitude

!write(6,*) shape(caliop_latitude), shape(mod03_longitude)
!write(6,*) Cell_Along_Swath_1km, Cell_Across_Swath_1km

!do i = 1, 10

!write(6,*) caliop_latitude(1,i), caliop_longitude(1,i)

!end do

!write(6,*)

!i = 650
!do j = 1, 10
!	write(6,*) mod03_latitude(i,j), mod03_longitude(i,j)
!end do
!write(6,*)

!----------------------------------------------------------------------------------------!
! Step 1: find the first modis-caliop collocated pixel

offset = len( trim(myd03_file_name) ) - 41 ! myd03 file name has 41 letters

heure = myd03_file_name(offset+16:offset+17)
minute = myd03_file_name(offset+18:offset+19)

write(fn,'("modis_caliop_cloudsat_collocated_indexes_    ")')

write(fn(42:43),'(a)') heure
write(fn(44:45),'(a)') minute

find_first_collocated_pixel = .False.

n = 1

do while( .not. find_first_collocated_pixel .and. n < caliop_dim_sizes )

	if( min_modis_latitude < caliop_latitude(1,n) .AND. &
	&   caliop_latitude(1,n) < max_modis_latitude )then
	if( min_modis_longitude < caliop_longitude(1,n) .AND. &
	&   caliop_longitude(1,n) < max_modis_longitude )then

		i_pr = 1
		j_pr = 1

		!--------------------------------------------------------------------------------!
		! Ellipsoidal Earth projected to a plane 
		
		diff_lat = caliop_latitude(1,n) - mod03_latitude(i_pr,j_pr)
		diff_long = caliop_longitude(1,n) - mod03_longitude(i_pr,j_pr)
	
		phim = ( mod03_latitude(i_pr,j_pr) + caliop_latitude(1,n) ) / 2.0
	
		K1 = 111.13209 - 0.56605 * cos( 2.0 * phim * cdr ) + &
		&    0.00120 * cos( 4.0 * phim * cdr )
	
		K2 = 111.41513 * cos( phim * cdr ) - 0.09455 * cos( 3.0 * phim * cdr ) + &
		&    0.00012 * cos( 5.0 * phim * cdr )
	
		tmp_01 = K1 * diff_lat
		tmp_02 = K2 * diff_long
	
		distance_pr = dsqrt( tmp_01**2 + tmp_02**2 )
		
		!--------------------------------------------------------------------------------!
	
		find_min_distance = .False.
		!do while( .not. find_min_distance  )
		do i = 1, 1354	
				
		!i = i_pr + 1
		j = j_pr
	
		!--------------------------------------------------------------------------------!
		! Ellipsoidal Earth projected to a plane 
		
		diff_lat = caliop_latitude(1,n) - mod03_latitude(i,j)
		diff_long = caliop_longitude(1,n) - mod03_longitude(i,j)
	
		phim = ( mod03_latitude(i,j) + caliop_latitude(1,n) ) / 2.0
	
		K1 = 111.13209 - 0.56605 * cos( 2.0 * phim * cdr ) + &
		&    0.00120 * cos( 4.0 * phim * cdr )
	
		K2 = 111.41513 * cos( phim * cdr ) - 0.09455 * cos( 3.0 * phim * cdr ) + &
		&    0.00012 * cos( 5.0 * phim * cdr )
	
		tmp_01 = K1 * diff_lat
		tmp_02 = K2 * diff_long
	
		distance = dsqrt( tmp_01**2 + tmp_02**2 )
				
		if( distance < distance_pr )then
			i_pr = i
			j_pr = j
			distance_pr = distance
			!write(6,*) i,j, distance	
		else
			find_min_distance = .True.
		end if
	
		!--------------------------------------------------------------------------------!
		
		end do

		!--------------------------------------------------------------------------------!	
				
		if( distance_pr < minimum_distance )then
			!minimum_distance = distance
			min_j = j
			collocation = .true.
			find_first_collocated_pixel = .True.		
		end if

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

	end if
	end if

n = n + 1

end do

!write(6,*) n-1,i_pr,j_pr, distance_pr 
!write(1,*) n-1,i_pr,j_pr, distance_pr 
!write(1,*) i_pr,j_pr, n-1

nb_collocated_pixels = 1

collocated_pixels(1,1) = i_pr
collocated_pixels(1,2) = j_pr
collocated_pixels(1,3) = n-1

!----------------------------------------------------------------------------------------!
! Step 2: iterate until no more collocated pixel can be found

tmp_n = 1

find_modis_pixel = .True.

do while( find_modis_pixel )

	!------------------------------------------------------------------------------------!
	! Ellipsoidal Earth projected to a plane 

	diff_lat = caliop_latitude(1,n) - mod03_latitude(i_pr,j_pr)
	diff_long = caliop_longitude(1,n) - mod03_longitude(i_pr,j_pr)

	phim = ( mod03_latitude(i_pr,j_pr) + caliop_latitude(1,n) ) / 2.0

	K1 = 111.13209 - 0.56605 * cos( 2.0 * phim * cdr ) + &
	&    0.00120 * cos( 4.0 * phim * cdr )

	K2 = 111.41513 * cos( phim * cdr ) - 0.09455 * cos( 3.0 * phim * cdr ) + &
	&    0.00012 * cos( 5.0 * phim * cdr )

	tmp_01 = K1 * diff_lat
	tmp_02 = K2 * diff_long

	distance_pr = dsqrt( tmp_01**2 + tmp_02**2 )

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

	find_min_distance = .False.

	do while( .not. find_min_distance  )

		distance_min = 999.9
		i_min = 1
		j_min = 1

		!--------------------------------------------------------------------------------!
		! neighbour 1

		i = i_pr + 1
		j = j_pr 

		if( i > 0 .and. i <=  Cell_Along_Swath_1km )then
		if( j > 0 .and. j <=  Cell_Along_Swath_1km )then

		diff_lat = caliop_latitude(1,n) - mod03_latitude(i,j)
		diff_long = caliop_longitude(1,n) - mod03_longitude(i,j)
	
		phim = ( mod03_latitude(i,j) + caliop_latitude(1,n) ) / 2.0
	
		K1 = 111.13209 - 0.56605 * cos( 2.0 * phim * cdr ) + &
		&    0.00120 * cos( 4.0 * phim * cdr )
	
		K2 = 111.41513 * cos( phim * cdr ) - 0.09455 * cos( 3.0 * phim * cdr ) + &
		&    0.00012 * cos( 5.0 * phim * cdr )
	
		tmp_01 = K1 * diff_lat
		tmp_02 = K2 * diff_long
	
		distance = dsqrt( tmp_01**2 + tmp_02**2 )

		if( distance < distance_min ) then
			distance_min = distance
			i_min = i
			j_min = j
		end if

		end if
		end if

		!--------------------------------------------------------------------------------!
		! neighbour 2

		i = i_pr + 1
		j = j_pr + 1

		if( i > 0 .and. i <=  Cell_Along_Swath_1km )then
		if( j > 0 .and. j <=  Cell_Along_Swath_1km )then
 
 		diff_lat = caliop_latitude(1,n) - mod03_latitude(i,j)
		diff_long = caliop_longitude(1,n) - mod03_longitude(i,j)
	
		phim = ( mod03_latitude(i,j) + caliop_latitude(1,n) ) / 2.0
	
		K1 = 111.13209 - 0.56605 * cos( 2.0 * phim * cdr ) + &
		&    0.00120 * cos( 4.0 * phim * cdr )
	
		K2 = 111.41513 * cos( phim * cdr ) - 0.09455 * cos( 3.0 * phim * cdr ) + &
		&    0.00012 * cos( 5.0 * phim * cdr )
	
		tmp_01 = K1 * diff_lat
		tmp_02 = K2 * diff_long
	
		distance = dsqrt( tmp_01**2 + tmp_02**2 )

		if( distance < distance_min ) then
			distance_min = distance
			i_min = i
			j_min = j
		end if
 
 		end if
		end if
 
		!--------------------------------------------------------------------------------!   
 		! neighbour 3
 		
		i = i_pr 
		j = j_pr + 1

		if( i > 0 .and. i <=  Cell_Along_Swath_1km )then
		if( j > 0 .and. j <=  Cell_Along_Swath_1km )then

		diff_lat = caliop_latitude(1,n) - mod03_latitude(i,j)
		diff_long = caliop_longitude(1,n) - mod03_longitude(i,j)
	
		phim = ( mod03_latitude(i,j) + caliop_latitude(1,n) ) / 2.0
	
		K1 = 111.13209 - 0.56605 * cos( 2.0 * phim * cdr ) + &
		&    0.00120 * cos( 4.0 * phim * cdr )
	
		K2 = 111.41513 * cos( phim * cdr ) - 0.09455 * cos( 3.0 * phim * cdr ) + &
		&    0.00012 * cos( 5.0 * phim * cdr )
	
		tmp_01 = K1 * diff_lat
		tmp_02 = K2 * diff_long
	
		distance = dsqrt( tmp_01**2 + tmp_02**2 )

		if( distance < distance_min ) then
			distance_min = distance
			i_min = i
			j_min = j
		end if

		end if
		end if

		!--------------------------------------------------------------------------------!
		! neighbour 4
		
		i = i_pr - 1
		j = j_pr + 1

		if( i > 0 .and. i <=  Cell_Along_Swath_1km )then
		if( j > 0 .and. j <=  Cell_Along_Swath_1km )then

		diff_lat = caliop_latitude(1,n) - mod03_latitude(i,j)
		diff_long = caliop_longitude(1,n) - mod03_longitude(i,j)
	
		phim = ( mod03_latitude(i,j) + caliop_latitude(1,n) ) / 2.0
	
		K1 = 111.13209 - 0.56605 * cos( 2.0 * phim * cdr ) + &
		&    0.00120 * cos( 4.0 * phim * cdr )
	
		K2 = 111.41513 * cos( phim * cdr ) - 0.09455 * cos( 3.0 * phim * cdr ) + &
		&    0.00012 * cos( 5.0 * phim * cdr )
	
		tmp_01 = K1 * diff_lat
		tmp_02 = K2 * diff_long
	
		distance = dsqrt( tmp_01**2 + tmp_02**2 )

		if( distance < distance_min ) then
			distance_min = distance
			i_min = i
			j_min = j
		end if

		end if
		end if

		!--------------------------------------------------------------------------------!
		! neighbour 5

		i = i_pr - 1 
		j = j_pr		

		if( i > 0 .and. i <=  Cell_Along_Swath_1km )then
		if( j > 0 .and. j <=  Cell_Along_Swath_1km )then

		diff_lat = caliop_latitude(1,n) - mod03_latitude(i,j)
		diff_long = caliop_longitude(1,n) - mod03_longitude(i,j)
	
		phim = ( mod03_latitude(i,j) + caliop_latitude(1,n) ) / 2.0
	
		K1 = 111.13209 - 0.56605 * cos( 2.0 * phim * cdr ) + &
		&    0.00120 * cos( 4.0 * phim * cdr )
	
		K2 = 111.41513 * cos( phim * cdr ) - 0.09455 * cos( 3.0 * phim * cdr ) + &
		&    0.00012 * cos( 5.0 * phim * cdr )
	
		tmp_01 = K1 * diff_lat
		tmp_02 = K2 * diff_long
	
		distance = dsqrt( tmp_01**2 + tmp_02**2 )

		if( distance < distance_min ) then
			distance_min = distance
			i_min = i
			j_min = j
		end if

		end if
		end if

		!--------------------------------------------------------------------------------!
		! neighbour 6

		i = i_pr - 1 
		j = j_pr - 1

		if( i > 0 .and. i <=  Cell_Along_Swath_1km )then
		if( j > 0 .and. j <=  Cell_Along_Swath_1km )then

		diff_lat = caliop_latitude(1,n) - mod03_latitude(i,j)
		diff_long = caliop_longitude(1,n) - mod03_longitude(i,j)
	
		phim = ( mod03_latitude(i,j) + caliop_latitude(1,n) ) / 2.0
	
		K1 = 111.13209 - 0.56605 * cos( 2.0 * phim * cdr ) + &
		&    0.00120 * cos( 4.0 * phim * cdr )
	
		K2 = 111.41513 * cos( phim * cdr ) - 0.09455 * cos( 3.0 * phim * cdr ) + &
		&    0.00012 * cos( 5.0 * phim * cdr )
	
		tmp_01 = K1 * diff_lat
		tmp_02 = K2 * diff_long
	
		distance = dsqrt( tmp_01**2 + tmp_02**2 )

		if( distance < distance_min ) then
			distance_min = distance
			i_min = i
			j_min = j
		end if

		end if
		end if

		!--------------------------------------------------------------------------------!
		! neighbour 7

		i = i_pr 
		j = j_pr - 1

		if( i > 0 .and. i <=  Cell_Along_Swath_1km )then
		if( j > 0 .and. j <=  Cell_Along_Swath_1km )then

		diff_lat = caliop_latitude(1,n) - mod03_latitude(i,j)
		diff_long = caliop_longitude(1,n) - mod03_longitude(i,j)
	
		phim = ( mod03_latitude(i,j) + caliop_latitude(1,n) ) / 2.0
	
		K1 = 111.13209 - 0.56605 * cos( 2.0 * phim * cdr ) + &
		&    0.00120 * cos( 4.0 * phim * cdr )
	
		K2 = 111.41513 * cos( phim * cdr ) - 0.09455 * cos( 3.0 * phim * cdr ) + &
		&    0.00012 * cos( 5.0 * phim * cdr )
	
		tmp_01 = K1 * diff_lat
		tmp_02 = K2 * diff_long
	
		distance = dsqrt( tmp_01**2 + tmp_02**2 )

		if( distance < distance_min ) then
			distance_min = distance
			i_min = i
			j_min = j
		end if

		end if
		end if

		!--------------------------------------------------------------------------------!
		! neighbour 8
		
		i = i_pr + 1
		j = j_pr - 1		

		if( i > 0 .and. i <=  Cell_Along_Swath_1km )then
		if( j > 0 .and. j <=  Cell_Along_Swath_1km )then

		diff_lat = caliop_latitude(1,n) - mod03_latitude(i,j)
		diff_long = caliop_longitude(1,n) - mod03_longitude(i,j)
	
		phim = ( mod03_latitude(i,j) + caliop_latitude(1,n) ) / 2.0
	
		K1 = 111.13209 - 0.56605 * cos( 2.0 * phim * cdr ) + &
		&    0.00120 * cos( 4.0 * phim * cdr )
	
		K2 = 111.41513 * cos( phim * cdr ) - 0.09455 * cos( 3.0 * phim * cdr ) + &
		&    0.00012 * cos( 5.0 * phim * cdr )
	
		tmp_01 = K1 * diff_lat
		tmp_02 = K2 * diff_long
	
		distance = dsqrt( tmp_01**2 + tmp_02**2 )

		if( distance < distance_min ) then
			distance_min = distance
			i_min = i
			j_min = j
		end if

		end if
		end if

		!write(6,*) distance_pr, i_pr, j_pr
		!write(6,*) distance_min, i_min, j_min
		!write(6,*) '-----------',  n, distance_min < distance_pr

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

		if( distance_min < distance_pr )then
			i_pr = i_min
			j_pr = j_min
			distance_pr = distance_min
		else
			find_min_distance = .True.
		end if

	end do

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

	if( distance_pr < minimum_distance )then
		!minimum_distance = distance
		min_j = j
		collocation = .true.
		find_first_collocated_pixel = .True.	
		!write(6,*) n,i_pr,j_pr, distance_pr, minimum_distance
		!write(1,*) n,i_pr,j_pr, distance_pr
		!write(1,*) i_pr,j_pr,n

		nb_collocated_pixels = nb_collocated_pixels + 1

		collocated_pixels(nb_collocated_pixels,1) = i_pr
		collocated_pixels(nb_collocated_pixels,2) = j_pr
		collocated_pixels(nb_collocated_pixels,3) = n

	end if

	if( distance_pr > 5.0 )then
		find_modis_pixel = .False.
	end if

	!tmp_n = tmp_n + 1
	!if( tmp_n > 10 ) STOP

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

n = n + 1

end do

!----------------------------------------------------------------------------------------!
! Add CloudSat

do i = 1, nb_collocated_pixels

i_pr = collocated_pixels(i,1)
i_pr = collocated_pixels(i,2)
n = collocated_pixels(i,3)

distance_min = 999.9

do j = 1, nb_cloudsat_pixels

		!--------------------------------------------------------------------------------!
		! Ellipsoidal Earth projected to a plane 
		
		diff_lat = caliop_latitude(1,n) - cloudsat_latitude(j)
		diff_long = caliop_longitude(1,n) - cloudsat_longitude(j)
	
		phim = ( cloudsat_latitude(j) + caliop_latitude(1,n) ) / 2.0
	
		K1 = 111.13209 - 0.56605 * cos( 2.0 * phim * cdr ) + &
		&    0.00120 * cos( 4.0 * phim * cdr )
	
		K2 = 111.41513 * cos( phim * cdr ) - 0.09455 * cos( 3.0 * phim * cdr ) + &
		&    0.00012 * cos( 5.0 * phim * cdr )
	
		tmp_01 = K1 * diff_lat
		tmp_02 = K2 * diff_long
	
		distance_pr = dsqrt( tmp_01**2 + tmp_02**2 )

		if( distance_pr < distance_min ) then
			distance_min = distance_pr
			j_min = j
		end if

end do

if( distance_min < minimum_distance )then
	collocated_pixels(i,4) = j_min
end if

end do

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

open(1,file=fn)

write(1,'(a)') caliop_file_name
write(1,'(a)') myd03_file_name
write(1,'(a)') cloudsat_file_name

do i = 1, nb_collocated_pixels
write(1,*) collocated_pixels(i,1),collocated_pixels(i,2),collocated_pixels(i,3),collocated_pixels(i,4)
end do

close(1)

end

Create a fortran module to read CloudSat 2B GEOPROF data

					
MODULE get_cloudsat_2b_geoprof_data

implicit none

!----------------------------------------------------------------------------------------!
! Dynamic Allocation of Arrays

real(kind=4), dimension(:), allocatable :: cloudsat_latitude
real(kind=4), dimension(:), allocatable :: cloudsat_longitude

integer*1, dimension(:,:), allocatable :: cloudsat_cloud_mask

!----------------------------------------------------------------------------------------!
! User Inputs (user_selected_sds)

integer :: nb_user_selected_sds
parameter ( nb_user_selected_sds = 1 )

character*300 selected_sds( nb_user_selected_sds )      

data selected_sds  /'CPR_Cloud_mask'/

!----------------------------------------------------------------------------------------!
! Subroutine read_cloudsat_2b_geoprof

CONTAINS

subroutine read_cloudsat_2b_geoprof(cloudsat_file_name)

use read_write_hdf_file_variables

implicit none

integer :: i,j,k,l,m,n
integer :: AllocateStatus, DeAllocateStatus
integer :: ierr ! etat de reussite (si ==0) ou echec (si /= 0)
character*300 :: cloudsat_file_name

integer sfstart, sfselect, sfrdata, sfendacc, sfend , sffattr , sfrattr, var_id, attr_id, err_code
integer status, n_attrs
integer sffinfo , sfginfo

character*300 vdata_name 
integer file_id, vdata_id, vdata_ref
integer n_records
character*80 fieldname_lists 
integer FULL_INTERLACE
parameter (FULL_INTERLACE = 0)
integer interlace_mode, vdata_size

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

file_id = hopen(cloudsat_file_name, DFACC_READ, 0)

status = vfstart(file_id) 

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

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

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

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

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

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

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

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

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

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

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

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

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

!----------------------------------------------------------------------------------------!
! SDS

sd_id = sfstart(cloudsat_file_name, DFACC_READ)
status = sffinfo( sd_id , n_datasets , n_file_attrs )

do index = 0, n_datasets - 1

	sds_id = sfselect(sd_id, index)
	status = sfginfo( sds_id , name , rank , dim_sizes , data_type , n_attrs )
    sds_name = name(1:len(name))

if ( ANY( selected_sds==sds_name ) ) then

	do i = 1 , rank
	start(i) = 0
	edges(i) = dim_sizes(i)
	stride(i) = 1
	end do
	
	if( sds_name == 'CPR_Cloud_mask' )then
	allocate( cloudsat_cloud_mask(dim_sizes(1),dim_sizes(2)), STAT = AllocateStatus)
	IF (AllocateStatus /= 0) STOP "*** Not enough memory (AllocateStatus cloudsat_cloud_mask)***"
	status = sfrdata( sds_id, start, stride, edges, cloudsat_cloud_mask )
	end if	
	
end if

end do

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

end subroutine read_cloudsat_2b_geoprof

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

END MODULE get_cloudsat_2b_geoprof_data

How to plot a contingency (confusion) table with python 3 and pandas ?

					
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sn
import pandas as pd

conf_arr = np.array([[ 60.36446469,   2.73348519],
                     [  2.73348519,  34.16856492]])

df_cm = pd.DataFrame(conf_arr, 
                  index = ['CALIOP (v4) LIQUID', 'CALIOP (v4) ICE'],
                  columns = ['MODIS LIQUID', 'MODIS ICE'])

fig = plt.figure()

plt.clf()

ax = fig.add_subplot(111)

ax.set_aspect(1)

res = sn.heatmap(df_cm, annot=True, vmin=0.0, vmax=100.0, fmt='.2f')

plt.yticks([0.5,1.5], ['CALIOP (v4) ICE', 'CALIOP (v4) LIQUID'],va='center')

plt.title('''Gaussian Process Cloud Phase Model \n Testing Dataset (July 2008) 
          -- Ocean, no snow-ice \n PAF = 0.94''')

plt.savefig('confusion_table.png', format='png', bbox_inches='tight')

Plot modis caliop collocated dataset pixel count on a global map with python 3

					
#!/usr/bin/env python

from mpl_toolkits.basemap import Basemap, cm

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

#----------------------------------------------------------------------------------------#
# Inputs

data = np.loadtxt("modis_c6_caliop_v4_monthly_l2_counts.txt", skiprows=1)

print(data.shape)

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

global_map_day_grid_counts = np.zeros((180,360))

index = 0 
for i in range(180):
	for j in range(360):
		global_map_day_grid_counts[i,j] = data[index,0]
		index = index + 1

max_value = np.max(global_map_day_grid_counts)
min_value = np.min(global_map_day_grid_counts)

print(min_value,max_value)

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

fig = plt.figure()

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

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

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


bounds = [0]+[i for i in np.linspace(1,max_value,10)]
norm = mpl.colors.BoundaryNorm(bounds, cmap.N)

img = m.imshow(global_map_day_grid_counts, cmap=cmap, norm=norm, interpolation='none')

cbar_bounds = [i for i in np.linspace(0,max_value,10)]
cbar_ticks = [i for i in np.linspace(0,max_value,10)]  

cb = m.colorbar(img, size="5%", pad='4%', boundaries=cbar_bounds, ticks=cbar_ticks)

plt.title("MODIS CALIOP Dataset (Day) Counts \n July 2008")

plt.savefig('modis_caliop_dataset_day_counts.png',  bbox_inches='tight')
#plt.show()



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

global_map_night_grid_counts = np.zeros((180,360))

index = 0 
for i in range(180):
	for j in range(360):
		global_map_night_grid_counts[i,j] = data[index,1]
		index = index + 1

max_value = np.max(global_map_night_grid_counts)
min_value = np.min(global_map_night_grid_counts)

print(min_value,max_value)

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

fig = plt.figure()

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

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

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


bounds = [0]+[i for i in np.linspace(1,max_value,10)]
norm = mpl.colors.BoundaryNorm(bounds, cmap.N)

img = m.imshow(global_map_night_grid_counts, cmap=cmap, norm=norm, interpolation='none')

cbar_bounds = [i for i in np.linspace(0,max_value,10)]
cbar_ticks = [i for i in np.linspace(0,max_value,10)]  

cb = m.colorbar(img, size="5%", pad='4%', boundaries=cbar_bounds, ticks=cbar_ticks)

plt.title("MODIS CALIOP Dataset (Night) Counts \n July 2008")

plt.savefig('modis_caliop_dataset_night_counts.png',  bbox_inches='tight')
#plt.show()

Create gridded pixel count of the collocated modis caliop dataset using fortran

					
program modis_c6_caliop_v4_monthly_l2_counts

! Note:
!
! Compilation
!
! make -f makefile.modis_c6_caliop_v4_monthly_l2_counts
! make clear -f makefile.modis_c6_caliop_v4_monthly_l2_counts
!
! Execution:
! run_modis_c6_caliop_v4_monthly_l2_counts

use read_write_hdf_file_variables 
use get_modis_c6_myd03_data
use get_modis_c6_myd06_data
use get_caliop_v4_01kmCLay_data
use get_caliop_v4_05kmCLay_data

implicit none 

integer, parameter :: pr = selected_real_kind(15,3)

integer :: i,j,k,l,m,n
integer :: mod_file
integer :: year,month,day,hour,minute,julian_day
character*300 :: path_calmod
character*100 :: calmod_index,calmod_file_name
character*400 :: input_colocated_file_index,input_colocated_file_name
integer :: ierror, iazm, iszen, iview
integer :: string_len
integer :: ierr ! successful (if == 0) ou failed (si /= 0)

character*50 :: machine_id = 'nasa_desktop' 
!character*50 :: machine_id = 'icare_server'

character*4 :: hhmm
character*300 :: file_name_collocated_pixels, file_name_modis_linked
character*300 :: path_to_modis_caliop_collocated_pixels, path_to_modis_linked_files
character*300 :: path_to_modis_caliop_collocated_pixels_01, path_to_modis_linked_files_01

integer :: nb_colocated_data

logical :: found_colocated_data, missing_file

integer, dimension(:,:), allocatable :: modis_index
integer, dimension(:), allocatable :: caliop_index

integer :: AllocateStatus, DeAllocateStatus

character*300 :: input_MYD03
character*300 :: input_MYD06_c6
character*300 :: input_CAL_01km
character*300 :: input_CAL_05km

character*16 :: file_common_path
integer, dimension(3) :: ShapeArray

integer :: s1,s2,s3,s4,s5,s6,s7,s8

integer, dimension(12) :: dpm_nyear, dpm_byear
data dpm_nyear /31,28,31,30,31,30,31,31,30,31,30,31/
data dpm_byear /31,29,31,30,31,30,31,31,30,31,30,31/

integer :: cloud_mask_flag
integer :: surface_flag
integer :: day_flag, night_flag
integer :: modis_cloud_phase_flag

integer :: nb_total_collocated_data 

real(pr) :: grid_resolution
integer :: grid_bin, idx, jdx, idx_max, jdx_max
integer, dimension(:,:), allocatable :: global_map_day_grid_counts 
integer, dimension(:,:), allocatable :: global_map_night_grid_counts 

!----------------------------------------------------------------------------------------!
! inputs

year = 2008
month = 07

julian_day = 0
do i = 1, month - 1
julian_day = julian_day + dpm_byear(i)
end do
julian_day = julian_day + day

if( trim(machine_id) == 'nasa_desktop' )then
  path_to_modis_caliop_collocated_pixels = '/Volumes/Data/Colocated_data/Caliop_modis/01kmCLay/'
  path_to_modis_linked_files = '/Volumes/Data/Colocated_data/Caliop_modis/01kmCLay/'
else 
  if( trim(machine_id) == 'icare_server' )then
    path_to_modis_caliop_collocated_pixels = '/home/benjamin.hg.marchant/data/modis_caliop_l2_collocated_pixels/'
    path_to_modis_linked_files = '/home/benjamin.hg.marchant/data/modis_linked_files/'
 else
   STOP
  end if
end if

file_name_collocated_pixels = 'calmod_index_'
file_name_modis_linked = 'calmod_file_name_'

!----------------------------------------------------------------------------------------!
! set-up paths

file_common_path = '    /    _  _  /'

write(file_common_path(1:4),'(i4.4)') year
write(file_common_path(6:9),'(i4.4)') year
write(file_common_path(11:12),'(i2.2)') month

!----------------------------------------------------------------------------------------!
! Initialization

nb_total_collocated_data = 0

grid_resolution = 1.0

idx_max = int( 180.0 / grid_resolution )
jdx_max = int( 360.0 / grid_resolution ) 

allocate( global_map_day_grid_counts( idx_max , jdx_max ) )
allocate( global_map_night_grid_counts( idx_max , jdx_max ) )

global_map_day_grid_counts = 0 
global_map_night_grid_counts = 0 

!----------------------------------------------------------------------------------------!
! start loop

do day = 1, 30

write(file_common_path(14:15),'(i2.2)') day

string_len = len( trim(path_to_modis_caliop_collocated_pixels) )
path_to_modis_caliop_collocated_pixels_01 = path_to_modis_caliop_collocated_pixels(1:string_len)//file_common_path

string_len = len( trim(path_to_modis_linked_files) )
path_to_modis_linked_files_01 = path_to_modis_linked_files(1:string_len)//file_common_path

hour = 0
minute = 0

do mod_file = 1, 288

hhmm = '    '

write(hhmm(1:2),'(i2.2)') hour
write(hhmm(3:4),'(i2.2)') minute

string_len = len( trim(file_name_collocated_pixels) )
input_colocated_file_index = file_name_collocated_pixels(1:string_len)//hhmm

string_len = len( trim(file_name_modis_linked) )
input_colocated_file_name = file_name_modis_linked(1:string_len)//hhmm

string_len = len( trim(path_to_modis_caliop_collocated_pixels_01) )
input_colocated_file_index = path_to_modis_caliop_collocated_pixels_01(1:string_len)//input_colocated_file_index

string_len = len( trim(path_to_modis_linked_files_01) )
input_colocated_file_name = path_to_modis_linked_files_01(1:string_len)//input_colocated_file_name

found_colocated_data = .TRUE.
open(1,file=input_colocated_file_index,form="formatted",iostat=ierr,status="old")
if (ierr/=0) then
!write(*,*) "Impossible d ouvrir le fichier ",trim(input_colocated_file_index)
!stop
found_colocated_data = .FALSE.
endif
nb_colocated_data = 0
do while (ierr==0)
read(1,*,iostat=ierr)
if (ierr==0) nb_colocated_data = nb_colocated_data + 1
enddo
close(1)

missing_file = .FALSE.
open(1,file=input_colocated_file_name,form="formatted",iostat=ierr,status="old")
if (ierr/=0) then
!write(*,*) "Impossible d ouvrir le fichier ",trim(input_colocated_file_name)
!stop
missing_file = .TRUE.
endif
close(1)

if( found_colocated_data .and. .not. missing_file )then
nb_colocated_data = nb_colocated_data - 2 
! First two lines are the name of caliop and modis files used to do the colocation
if( nb_colocated_data > 0 )then

allocate( modis_index(nb_colocated_data,2), STAT = AllocateStatus )
IF (AllocateStatus /= 0) STOP "*** Not enough memory modis_index***"
allocate( caliop_index(nb_colocated_data), STAT = AllocateStatus )
IF (AllocateStatus /= 0) STOP "*** Not enough memory caliop_index***"

open(1,file=input_colocated_file_index, status="old")
read(1,'(a)') 
read(1,'(a)') 
do i = 1, nb_colocated_data 
read(1,*) modis_index(i,1), modis_index(i,2), caliop_index(i)
end do
close(1)

if( trim(machine_id) == 'nasa_desktop' )then
  open(1,file=input_colocated_file_name, status="old")
  read(1,'(a)') input_MYD03
  read(1,'(a)') 
  read(1,'(a)') 
  read(1,'(a)') input_MYD06_c6
  read(1,'(a)') 
  read(1,'(a)') 
  read(1,'(a)') input_CAL_01km
  read(1,'(a)') input_CAL_05km
  read(1,'(a)') 
  close(1)
  input_MYD03 = '/Volumes/Data/Myd03/'//file_common_path//input_MYD03
  input_MYD06_c6 = '/Volumes/DATA_03/modis/Myd06_c6/'//file_common_path//input_MYD06_c6
  input_CAL_01km = '/Volumes/Data/Caliop/v4/01km_CLay/'//file_common_path//input_CAL_01km
  input_CAL_05km = '/Volumes/Data/Caliop/v4/05km_CLay/'//file_common_path//input_CAL_05km
else 
  if( trim(machine_id) == 'icare_server' )then
    open(1,file=input_colocated_file_name, status="old")
    read(1,'(a)') input_MYD03
    read(1,'(a)') 
    read(1,'(a)') 
    read(1,'(a)') input_MYD06_c6
    read(1,'(a)') input_CAL_01km
    read(1,'(a)') input_CAL_05km
    read(1,'(a)') 
    close(1)
    input_MYD03 = '/DATA/LIENS/MODIS/MYD03/'//file_common_path//input_MYD03
    input_MYD06_c6 = '/DATA/LIENS/MODIS/MYD06_L2.006/'//file_common_path//input_MYD06_c6
    input_CAL_01km = '/DATA/LIENS/CALIOP/01kmCLay.v4.10/'//file_common_path//input_CAL_01km
    input_CAL_05km = '/DATA/LIENS/CALIOP/05kmCLay.v4.10/'//file_common_path//input_CAL_05km
 else
   STOP
  end if
end if

CALL read_modis_c6_myd03(input_MYD03)
CALL read_modis_c6_myd06(input_MYD06_c6)
CALL read_caliop_v4_01kmCLay(input_CAL_01km)
CALL read_caliop_v4_05kmCLay(input_CAL_05km)

!----------------------------------------------------------------------------------------!
! start analysis

ShapeArray = shape(modis_c6_myd06_Cloud_Mask_1km)

do i = 1, nb_colocated_data 

nb_total_collocated_data = nb_total_collocated_data + 1

j = caliop_index(i)
k = modis_index(i,1)
l = modis_index(i,2)

if( ibits(modis_c6_myd06_Cloud_Mask_1km(1,k,l),0,1) == 1 .and. & ! Cloud Mask determined 
&   ibits(modis_c6_myd06_Cloud_Mask_1km(1,k,l),1,2) < 2 )then  ! Cloudy or Probably Cloudy Pixel  
cloud_mask_flag = 1
else
cloud_mask_flag = 2
end if

day_flag = 0
if( modis_c6_myd03_solarzenith(k,l) * 0.01 < 80.0 ) day_flag = 1 

night_flag = 0 
if( modis_c6_myd03_solarzenith(k,l) * 0.01 > 100.0 ) night_flag = 1

modis_cloud_phase_flag = 0
if(modis_c6_myd06_Cloud_Phase_Optical_Properties(k,l) == 2) modis_cloud_phase_flag = 1
if(modis_c6_myd06_Cloud_Phase_Optical_Properties(k,l) == 3) modis_cloud_phase_flag = 2

surface_flag = 0 
if( ibits(modis_c6_myd06_Cloud_Mask_1km(1,k,l),6,2) == 0 ) surface_flag = 0 ! Ocean
if( ibits(modis_c6_myd06_Cloud_Mask_1km(1,k,l),6,2) == 3 ) surface_flag = 3 ! Land
if( ibits(modis_c6_myd06_Cloud_Mask_1km(1,k,l),6,2) == 2 ) surface_flag = 2 ! Desert
if( ibits(modis_c6_myd06_Cloud_Mask_1km(1,k,l),6,2) == 1 ) surface_flag = 1 ! Coastal

!----------------------------------------------------------------------------------------!
! Counts

idx = int( (modis_c6_myd03_latitude(k,l) + 90.0) / grid_resolution ) + 1
jdx = int( (modis_c6_myd03_longitude(k,l) + 180.0) / grid_resolution ) + 1

if( idx > idx_max) idx = idx_max
if( jdx > jdx_max) jdx = jdx_max

if(day_flag==1)then
  global_map_day_grid_counts(idx,jdx) = global_map_day_grid_counts(idx,jdx) + 1
 end if
 
if(night_flag==1)then
  global_map_night_grid_counts(idx,jdx) = global_map_night_grid_counts(idx,jdx) + 1
 end if

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

end do !----- nb_colocated_data granule loop

deallocate(modis_c6_myd03_latitude, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd03_latitude) ***"
deallocate(modis_c6_myd03_longitude, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd03_longitude) ***"
deallocate(modis_c6_myd03_sensorzenith, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd03_sensorzenith) ***"
deallocate(modis_c6_myd03_sensorazimuth, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd03_sensorazimuth) ***"
deallocate(modis_c6_myd03_solarzenith, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd03_solarzenith) ***"
deallocate(modis_c6_myd03_solarazimuth, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd03_solarazimuth) ***"

deallocate(modis_c6_myd06_Cloud_Phase_Infrared_1km, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Phase_Infrared_1km) baseline***"
deallocate(modis_c6_myd06_Cloud_Phase_Optical_Properties, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Phase_Optical_Properties) baseline***"
deallocate(modis_c6_myd06_Cloud_Mask_1km, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Mask_1km) baseline***"
deallocate(modis_c6_myd06_Quality_Assurance_1km, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Quality_Assurance_1km) baseline***"
deallocate(modis_c6_myd06_Cloud_Optical_Thickness, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Optical_Thickness) baseline***"
deallocate(modis_c6_myd06_Cloud_Optical_Thickness_PCL, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Optical_Thickness_PCL) baseline***"
deallocate(modis_c6_myd06_cloud_top_temperature_1km, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_cloud_top_temperature_1km) baseline***"
deallocate(modis_c6_myd06_Cloud_Top_Temperature, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Top_Temperature) baseline***"
deallocate(modis_c6_myd06_Cloud_Multi_Layer_Flag, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Multi_Layer_Flag) baseline***"
deallocate(modis_c6_myd06_Cloud_Effective_Radius, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Effective_Radius) baseline***"
deallocate(modis_c6_myd06_Cloud_Effective_Radius_16, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Effective_Radius_16) baseline***"
deallocate(modis_c6_myd06_Cloud_Effective_Radius_37, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Effective_Radius_37) baseline***"
deallocate(modis_c6_myd06_Cloud_Top_Height_1km, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Top_Height_1km) baseline***"

deallocate(caliop_v4_01kmCLay_latitude, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_01kmCLay_latitude)***"
deallocate(caliop_v4_01kmCLay_longitude, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_01kmCLay_longitude)***"
deallocate(caliop_v4_01kmCLay_Opacity_Flag, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_01kmCLay_Opacity_Flag)***"
deallocate(caliop_v4_01kmCLay_Number_Layers_Found, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_01kmCLay_Number_Layers_Found)***"
deallocate(caliop_v4_01kmCLay_Feature_Classification_Flags, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_01kmCLay_Feature_Classification_Flags)***"
deallocate(caliop_v4_01kmCLay_Layer_Top_Altitude, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_01kmCLay_Layer_Top_Altitude)***"
deallocate(caliop_v4_01kmCLay_Layer_Base_Altitude, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_01kmCLay_Layer_Base_Altitude)***"

deallocate(caliop_v4_05kmCLay_Number_Layers_Found, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_05kmCLay_Number_Layers_Found)***"
deallocate(caliop_v4_05kmCLay_Feature_Classification_Flags, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_05kmCLay_Feature_Classification_Flags)***"
deallocate(caliop_v4_05kmCLay_Feature_Optical_Depth_532, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_05kmCLay_Feature_Optical_Depth_532)***"
deallocate(caliop_v4_05kmCLay_Layer_Top_Altitude, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_05kmCLay_Layer_Top_Altitude)***"
deallocate(caliop_v4_05kmCLay_Layer_Base_Altitude, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_05kmCLay_Layer_Base_Altitude)***"

deallocate(modis_index, STAT = DeAllocateStatus)
deallocate(caliop_index, STAT = DeAllocateStatus)

end if !----- if nb_colocated_data > 0
end if !----- if found colocated file

write(6,*) 'hour, minute: ', hour, minute
minute = minute + 5
if ( minute > 55 ) then
hour = hour + 1
minute = 0
end if
end do !----- Granule loop
write(6,*) 'day: ', day
end do !----- day loop

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

open(1,file='./outputs/modis_c6_caliop_v4_monthly_l2_counts_f90/modis_c6_caliop_v4_monthly_l2_counts.txt')
write(1,*) idx_max, jdx_max
do idx = 1, idx_max 
	do jdx = 1, jdx_max 
		write(1,*) global_map_day_grid_counts(idx,jdx), global_map_night_grid_counts(idx,jdx)
	end do
end do
close(1)

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

end program modis_c6_caliop_v4_monthly_l2_counts

Create a RGB image from MODIS Airborne Simulator L1 using python 3

					
#!/usr/bin/env python

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

#----------------------------------------------------------------------------------------#
# inputs

file_name = 'MASL1B_07915_05_20070707_1722_1737_V04.hdf'

file = SD(file_name, SDC.READ)

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

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

data = sds_obj.get() 

attr = sds_obj.attributes()

scale_factor_list = attr['scale_factor']

data_shape = data.shape

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

rgb = np.zeros((data_shape[0], data_shape[2],3))

rgb[:,:,0] = data[:,19,:]*scale_factor_list[19] 
rgb[:,:,1] = data[:,8,:]*scale_factor_list[8]  
rgb[:,:,2] = data[:,0,:]*scale_factor_list[0] 

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

#rgb = rgb / rgb.max()

rgb[:,:,0] = rgb[:,:,0] / rgb[:,:,0].max()
rgb[:,:,1] = rgb[:,:,1] / rgb[:,:,1].max() 
rgb[:,:,2] = rgb[:,:,2] / rgb[:,:,2].max()

rgb[ rgb > 1 ] = 1.0
rgb[ rgb < 0 ] = 0.0

#----------------------------------------------------------------------------------------#
# plot image using matplotlib

fig = plt.figure()

ax = fig.add_subplot(111)

img = plt.imshow(np.fliplr(rgb), interpolation='nearest', aspect=0.2)

plt.title('MODIS Airborne Simulator \n RGB image (bands 1,9,20)', fontsize=10)

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

plt.close()

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: