read_mod_aerosol_at_a_location.py

#!/usr/bin/python
'''
Module: read_mod_aerosol_at_a_location.py
==========================================================================================
Disclaimer: The code is for demonstration purposes only. Users are responsible to check for accuracy and revise to fit their objective.

Author: Justin Roberts-Pierel, 2015 
Organization: NASA ARSET
Purpose: To view info about a variety of SDS from a MODIS HDF4 file (or series of files) both generally and at a specific lat/lon

See the README associated with this module for more information.
==========================================================================================
'''

#import necessary modules
from pyhdf import SD
import numpy as np

#This uses the file "fileList.txt", containing the list of files, in order to read the files
try:
	fileList=open('fileList.txt','r')
except:
	print('Did not find a text file containing file names (perhaps name does not match)')
	sys.exit()

#loops through all files listed in the text file
for FILE_NAME in fileList:
	FILE_NAME=FILE_NAME.strip()
	user_input=input('\nWould you like to process\n' + FILE_NAME + '\n\n(Y/N)')
	if(user_input == 'N' or user_input == 'n'):
		continue
	else:
		if '3K' in FILE_NAME: #then this is a 3km MODIS file
			userInput=int(input('Which SDS would you like to view? (Type the number and press enter) \n(1) Optical_Depth_Land_And_Ocean \n(2) Land_Ocean_Quality_Flag \n(3) Image_Optical_Depth_Land_And_Ocean \n(4) Land_Sea_Flag \n'))
			while userInput not in {1,2,3,4}:#repeats the question if the user does not choose one of the options
				print('Please try again.')
				userInput=int(input('Which SDS would you like to view? (Type the number and press enter) \n(1) Optical_Depth_Land_And_Ocean \n(2) Land_Ocean_Quality_Flag \n(3) Image_Optical_Depth_Land_And_Ocean \n(4) Land_Sea_Flag \n'))
			#Uses a Python dictionary to choose the SDS indicated by the user
			dataFields=dict([(1,'Optical_Depth_Land_And_Ocean'),(2,'Land_Ocean_Quality_Flag'),(3,'Image_Optical_Depth_Land_And_Ocean'),(4,'Land_Sea_Flag')])
		elif 'L2' in FILE_NAME:#Same as above but for 10km MODIS file
			userInput=int(input('Which SDS would you like to view? (Type the number and press enter) \n(1) Deep_Blue_Aerosol_Optical_Depth_550_Land \n(2) AOD_550_Dark_Target_Deep_Blue_Combined \n(3) AOD_550_Dark_Target_Deep_Blue_Combined_QA_Flag \n'))
			while userInput not in {1,2,3}:
				print('Please try again.')
				userInput=int(input('Which SDS would you like to view? (Type the number and press enter) \n(1) Deep_Blue_Aerosol_Optical_Depth_550_Land \n(2) AOD_550_Dark_Target_Deep_Blue_Combined \n(3) AOD_550_Dark_Target_Deep_Blue_Combined_QA_Flag \n'))
			dataFields=dict([(1,'Deep_Blue_Aerosol_Optical_Depth_550_Land'),(2,'AOD_550_Dark_Target_Deep_Blue_Combined'),(3,'AOD_550_Dark_Target_Deep_Blue_Combined_QA_Flag')])
		SDS_NAME=dataFields[int(userInput)] # The name of the sds to read
		try:
			# open the hdf file for reading
			hdf=SD.SD(FILE_NAME)
		except:
			print('Unable to open file: \n' + FILE_NAME + '\n Skipping...')
			continue
		
		# Get lat and lon info
		lat = hdf.select('Latitude')
		latitude = lat[:,:]
		min_lat=latitude.min()
		max_lat=latitude.max()
		lon = hdf.select('Longitude')
		longitude = lon[:,:]
		min_lon=longitude.min()
		max_lon=longitude.max()
		
		#get SDS, or exit program if SDS is not in the file
		try:
			sds=hdf.select(SDS_NAME)
		except:
			print('Sorry, your MODIS hdf file does not contain the SDS:',SDS_NAME,'. Please try again with the correct file type.')
			continue
		#get scale factor and fill value for data field
		attributes=sds.attributes()
		scale_factor=attributes['scale_factor']
		fillvalue=attributes['_FillValue']
		
		#get SDS data
		data=sds.get()
		#Print the range of latitude and longitude found in the file, then ask for a lat and lon
		print('The range of latitude in this file is: ',min_lat,' to ',max_lat, 'degrees \nThe range of longitude in this file is: ',min_lon, ' to ',max_lon,' degrees')
		user_lat=float(input('\nPlease enter the latitude you would like to analyze (Deg. N): '))
		user_lon=float(input('Please enter the longitude you would like to analyze (Deg. E): '))
		#Continues to ask for lat and lon until the user enters valid values
		while user_lat < min_lat or user_lat > max_lat:
			user_lat=float(input('The latitude you entered is out of range. Please enter a valid latitude: '))
		while user_lon < min_lon or user_lon > max_lon:
			user_lon=float(input('The longitude you entered is out of range. Please enter a valid longitude: '))
			
		#calculation to find nearest point in data to entered location (haversine formula)
		R=6371000#radius of the earth in meters
		lat1=np.radians(user_lat)
		lat2=np.radians(latitude)
		delta_lat=np.radians(latitude-user_lat)
		delta_lon=np.radians(longitude-user_lon)
		a=(np.sin(delta_lat/2))*(np.sin(delta_lat/2))+(np.cos(lat1))*(np.cos(lat2))*(np.sin(delta_lon/2))*(np.sin(delta_lon/2))
		c=2*np.arctan2(np.sqrt(a),np.sqrt(1-a))
		d=R*c
		#gets (and then prints) the x,y location of the nearest point in data to entered location, accounting for no data values
		x,y=np.unravel_index(d.argmin(),d.shape)
		print('\nThe nearest pixel to your entered location is at: \nLatitude:',latitude[x,y],' Longitude:',longitude[x,y])
		if data[x,y]==fillvalue:
			print('The value of ',SDS_NAME,'at this pixel is',fillvalue,',(No Value)\n')
		else:
			print('The value of ', SDS_NAME,'at this pixel is ',round(data[x,y]*scale_factor,3))
			
		#calculates mean, median, stdev in a 3x3 grid around nearest point to entered location
		if x < 1:
			x+=1
		if x > data.shape[0]-2:
			x-=2
		if y < 1:
			y+=1
		if y > data.shape[1]-2:
			y-=2
		three_by_three=data[x-1:x+2,y-1:y+2]
		three_by_three=three_by_three.astype(float)
		three_by_three[three_by_three==float(fillvalue)]=np.nan
		nnan=np.count_nonzero(~np.isnan(three_by_three))
		if nnan == 0:
			print ('\nThere are no valid pixels in a 3x3 grid centered at your entered location.')
		else:
			three_by_three=three_by_three*scale_factor
			three_by_three_average=np.nanmean(three_by_three)
			three_by_three_std=np.nanstd(three_by_three)
			three_by_three_median=np.nanmedian(three_by_three)
			if nnan == 1:
				npixels='is'
				mpixels='pixel'
			else:
				npixels='are'
				mpixels='pixels'
			print('\nThere',npixels,nnan,'valid',mpixels,'in a 3x3 grid centered at your entered location.')
			print('\nThe average value in this grid is: ',round(three_by_three_average,3),' \nThe median value in this grid is: ',round(three_by_three_median,3),'\nThe standard deviation in this grid is: ',round(three_by_three_std,3))
		
		#calculates mean, median, stdev in a 5x5 grid around nearest point to entered location
		if x < 2:
			x+=1
		if x > data.shape[0]-3:
			x-=1
		if y < 2:
			y+=1
		if y > data.shape[1]-3:
			y-=1
		five_by_five=data[x-2:x+3,y-2:y+3]
		five_by_five=five_by_five.astype(float)
		five_by_five[five_by_five==float(fillvalue)]=np.nan
		nnan=np.count_nonzero(~np.isnan(five_by_five))
		if nnan == 0:
			print ('There are no valid pixels in a 5x5 grid centered at your entered location. \n')
		else:
			five_by_five=five_by_five*scale_factor
			five_by_five_average=np.nanmean(five_by_five)
			five_by_five_std=np.nanstd(five_by_five)
			five_by_five_median=np.nanmedian(five_by_five)
			if nnan == 1:
				npixels='is'
				mpixels='pixel'
			else:
				npixels='are'
				mpixels='pixels'
			print('\nThere',npixels,nnan,' valid',mpixels,' in a 5x5 grid centered at your entered location. \n')
			print('The average value in this grid is: ',round(five_by_five_average,3),' \nThe median value in this grid is: ',round(five_by_five_median,3),'\nThe standard deviation in this grid is: ',round(five_by_five_std,3))
print('\nAll valid files have been processed')