-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathbatchEDR2geotiff.py
171 lines (129 loc) · 6.76 KB
/
batchEDR2geotiff.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Dec 18 22:58:41 2018
@author: leonidas
"""
import h5py
import numpy as np
from pyresample import geometry
import os,glob
import xarray as xr
import dask.array as da
from satpy import Scene
from satpy.utils import debug_on
import areaSettings
import argparse
import ntpath
def maskByte(byte1, byte2):
# Joint Polar Satellite System (JPSS)-Algorithm Specification Volume II: Data-Dictionary for the Cloud Mask:
# https://jointmission.gsfc.nasa.gov/sciencedocs/2016-12/474-00448-02-11_JPSS-DD-Vol-II-Part-11_0200E.pdf
# Cloud Detection and Confidence Pixel
CloudDetection_ConfidencePixel_mask=(byte1 & 12)>>2 #apply bitwise operations
# export clouds DN=1, noclouds=0
# Confidently Clear = 0
# Probably Clear = 1
# Probably Cloudy = 2
# Confidently Cloudy = 3
# keep only Confidently Clear
# CloudDetection_ConfidencePixel_mask= ma.masked_greater(CloudDetection_ConfidencePixel_mask, 0)
# CloudDetection_ConfidencePixel_mask = ma.filled(CloudDetection_ConfidencePixel_mask, fill_value=1)
# Snow/Ice Surface Pixel
# No Snow/Ice 0
#Snow/Ice 1
SnowIce_mask=(byte1 & 32)>>5 #apply bitwise operations
Fire_mask=(byte2 & 32)>>5 #apply bitwise operations
# HeavyAerosol_mask=(QF2_VIIRSCMEDR & 16)>>4
# Cirrus_mask=(QF2_VIIRSCMEDR & 64)>>6
# CirrusIR_mask=(QF2_VIIRSCMEDR & 128)>>7
#
# Dust_Candidate_mask=(QF4_VIIRSCMEDR & 16)>>4
# Smoke_Candidate_mask=(QF4_VIIRSCMEDR & 32)>>5
# Dust_VolcanicAsh_mask=(QF4_VIIRSCMEDR & 64)>>6
#Combine pixel level flags
mask=np.logical_or.reduce((
CloudDetection_ConfidencePixel_mask,
SnowIce_mask,
Fire_mask
# HeavyAerosol_mask,
# Cirrus_mask,
# CirrusIR_mask,
# Dust_Candidate_mask,
# Smoke_Candidate_mask,
# Dust_VolcanicAsh_mask,
))
return mask
def readhdfDatasets(HDF):
with h5py.File(HDF,'r') as h5_file:
GROUP_DATA='All_Data/VIIRS-CM-EDR_All'
QF1_VIIRSCMEDR=h5_file[GROUP_DATA]['QF1_VIIRSCMEDR'][...]
QF2_VIIRSCMEDR=h5_file['All_Data/VIIRS-CM-EDR_All']['QF2_VIIRSCMEDR'][...]
#QF4_VIIRSCMEDR=h5_file['All_Data/VIIRS-CM-EDR_All']['QF4_VIIRSCMEDR'][...]
GROUP_GEODATA='All_Data/VIIRS-MOD-GEO_All'
lon_data=h5_file[GROUP_GEODATA]['Longitude'][...]
lat_data=h5_file[GROUP_GEODATA]['Latitude'][...]
return QF1_VIIRSCMEDR, QF2_VIIRSCMEDR, lon_data,lat_data
def EDR2Geotiff(HDF, output_dir,areaid, radius):
QF1_VIIRSCMEDR, QF2_VIIRSCMEDR, lon_data, lat_data = readhdfDatasets(HDF)
HDF = ntpath.basename(HDF)#get filename without path
#https://pytroll.slack.com/archives/C06GJFRN0/p1545083373181100
lon_data[QF1_VIIRSCMEDR == 0] = np.nan
lat_data[QF1_VIIRSCMEDR == 0] = np.nan
mask=maskByte(byte1=QF1_VIIRSCMEDR, byte2=QF2_VIIRSCMEDR)
mask=mask.astype(np.uint8)
fill_value=255 #fill_value is parameter of save_datasets(...). satpy sets 255 for pixels not included in my AOI.
swath_def = geometry.SwathDefinition(
xr.DataArray(da.from_array(lon_data, chunks=4096), dims=('y', 'x')),
xr.DataArray(da.from_array(lat_data, chunks=4096), dims=('y', 'x')))
metadata_dict = {'name': 'mask', 'area':swath_def}
scn = Scene()
scn['mask'] = xr.DataArray(
da.from_array(mask, chunks=4096),
attrs=metadata_dict,
dims=('y', 'x')) #https://satpy.readthedocs.io/en/latest/dev_guide/xarray_migration.html#id1
scn.load(["mask"])
proj_scn = scn.resample(areaSettings.getarea(areaid),radius_of_influence=radius)
proj_scn.save_datasets(writer='geotiff',base_dir=output_dir ,file_pattern="{}.{}.{}".format(HDF,"{name}","tif"),enhancement_config=False,
dtype=np.uint8, fill_value=fill_value) #
def batchEDR2Geotiff(input_dir, output_dir,pattern, areaid, debug):
# ********* SETTINGS *******************************************
#input_dir='/media/leonidas/Hitachi/daily_viirs/2017_packed/EDR_CLOUD_MASK'
#output_dir="geotiffs_greece2"
# pattern='GMODO-VICMO_npp_d20170331_t2346150_e2351554_b28114_c20181105104544883294_noac_ops.h5'
#'GMODO-VICMO_npp_d201703*.h5'
#'GMODO-VICMO_npp_d201703*.h5'
#'GMODO-VICMO_npp_d20170328_t0056143_e0101547_b28058_c20181105104553745389_noac_ops.h5'#GMODO-VICMO_npp_d20170328_t0056143_e0101547_b28058_c20181105104553745389_noac_ops.h5'##'GMODO-VICMO_npp_d201703*.h5'#'GMODO-VICMO_npp_d20170328_t0056143_e0101547_b28058_c20181105104553745389_noac_ops.h5'#GMODO-VICMO_npp_d20170313_t0039560_e0045364_b27845_c20181105104612702584_noac_ops.h5
if debug:
debug_on()
curdir=os.getcwd()
os.chdir(input_dir)
#OUTPUT_DIR=os.path.join(input_dir,OUTPUT_DIR_RELATIVE)
output_dir = os.path.abspath(output_dir) #works for abs and relative paths
if not os.path.exists(output_dir):
os.makedirs(output_dir)
hdf_files = glob.glob(os.path.join(input_dir,pattern))
os.chdir(curdir)
for hdf in hdf_files:
EDR2Geotiff(hdf, output_dir, areaid,radius=2000)
def dir_path(string):
string = os.path.abspath(string)
if os.path.isdir(string):
return string
else:
print('Provide a valid directory.')
raise NotADirectoryError(string)
def parse_args():
parser = argparse.ArgumentParser(description='Extracts Cloud Mask from VIIRS EDR hdf file (Specific Quality flags) as binary geotiff.0=Confidently Clear')
parser.add_argument("-i", "--inputdir", help="Directory contains hdf5 files", required=True, type=dir_path)
parser.add_argument("-o", "--outputdir", help="Directory to export geotiff files",required=True, type=dir_path)
parser.add_argument("-p", "--pattern", help="Regex pattern to match hdf files", default='GMODO-VICMO_npp_d*.h5')
parser.add_argument("-a", "--areaid", help="Area_ID",required=True)
parser.add_argument("-d", "--debug", action='store_true')
args = parser.parse_args()
return args
def main():
args = parse_args()
batchEDR2Geotiff(input_dir=args.inputdir, output_dir=args.outputdir, pattern=args.pattern, areaid=args.areaid, debug=args.debug)
if __name__ == '__main__':
main()
#runfile('/media/leonidas/Hitachi/daily_viirs/2017_packed/batchHDF2Geotiff/batchEDR2geotiff.py', wdir='/media/leonidas/Hitachi/daily_viirs/2017_packed/batchHDF2Geotiff' ,args = '-i /media/leonidas/Hitachi/daily_viirs/2017_packed/EDR_CLOUD_MASK -o /media/leonidas/Hitachi/daily_viirs/2017_packed/EDR_CLOUD_MASK/geotiffs_greece2 -p GMODO-VICMO_npp_d201706*.h5 -a greek_grid4')