Process '.bil' RPISM data to netcdf files.

process_prism.py

import xcdat as xc
import xarray as xa
import glob
import zipfile
from osgeo import gdal
import numpy as np
import argparse
import os

def get_args():
    parser = argparse.ArgumentParser()
    parser.add_argument('--month', type=int)
    parser.add_argument('--year', type=int, default=0)
    parser.add_argument('--var', type=str, choices=['tmax', 'tmin', 'ppt'])
    args = vars(parser.parse_args())
    return args

# PRISM lat lon:
lat = [49.9166666666664-0.0416666666667*i for i in range(621)]
lon = [-125+0.0416666666667*i for i in range(1405)]

args = get_args()
month = args['month']
year = args['year']
var = args['var']
if not os.path.exists('netcdf_data/'+var+'/'):
    os.makedirs('netcdf_data/'+var+'/', exist_ok=True)

month_data = []
filepath = glob.glob('original_data/prism.oregonstate.edu/daily/'+var+'/'+str(year)+'/PRISM_'+var+'_stable_4kmD2_'+str(year)+str(month).zfill(2)+'*_bil.zip')
for file in sorted(filepath):
    files = glob.glob('tmp/PRISM*')
    for f in files:
        os.remove(f)
        # os.remove() # remove previous temp files. 
    print(file)
    with zipfile.ZipFile(file, 'r') as zip_ref:
        zip_ref.extractall('tmp/')
    bil_file = glob.glob('tmp/PRISM_'+var+'_stable_4kmD2_*_bil.bil')[0]
    dataset = gdal.Open(bil_file)
    # Read data
    num_bands = dataset.RasterCount
    band = dataset.GetRasterBand(1)
    data = band.ReadAsArray()
    print(data.shape)
    print(num_bands)

    data_flat = data.flatten()
    data_flat[data_flat<-9990] = np.NAN
    data_plot = data_flat.reshape(data.shape)
    data_xa = xa.DataArray(data=data_plot.reshape(1, data.shape[0], data.shape[1]),
                        dims=['time', 'lat', 'lon'], 
                        coords={'time':[np.datetime64('1981-01-01')], 
                                'lat':lat, 'lon':lon},
                        name='tmax'
                        )
    month_data.append(data_xa)
month_data = xa.concat(month_data, dim='time')
month_data.to_netcdf('netcdf_data/'+var+'/'+var+str(year)+str(month).zfill(2)+'.nc')