Commit 0159d967 authored by Coda Phillips's avatar Coda Phillips
Browse files

Add solar and timing

parent 7e11fab7
......@@ -13,4 +13,5 @@ meteosat8/
notebooks/.ipynb_checkpoints/
satzen_cache/
xrit/
.ipynb_checkpoints/
composite_cache*
demo_20210908 - Latest 0.1 degree data
demo_20211115 - Latest 0.05 degree data
......@@ -22,6 +22,12 @@ M11_ROOT = Path('/arcdata/nongoes/meteosat/meteosat11/')
ROOTS = {'g16':G16_ROOT, 'g17':G17_ROOT, 'h8':H8_ROOT, 'm8':M8_ROOT, 'm11':M11_ROOT}
def band_dir_path(dt, sat, band, l1b_dir=L1B_DIR):
band_dir = l1b_dir/dt.strftime('%Y')/dt.strftime('%Y%m')/dt.strftime('%Y%m%d')/dt.strftime('%Y%m%dT%H%M')/sat/f'{band}'
if not band_dir.is_dir():
raise IOError(f"Missing {band_dir}")
return band_dir
def parse_or(fname):
_,band,sat,start,end,create = fname.split('.')[0].split('_')
start = datetime.strptime(start[:-1], 's%Y%j%H%M%S')
......@@ -75,8 +81,11 @@ def filter_fd_goes(files, start, end):
return pd.DataFrame()
longest = candidates.groupby(['start','band']).apply(lambda x: x.sort_index().iloc[-1])
sizes = longest.groupby('start').size()
full = longest.loc[[sizes.index[sizes == 16][0]]]
return full
if (sizes == 16).any():
full = longest.loc[[sizes.index[sizes == 16][0]]]
return full
else:
return []
def filter_fd_him(files, start, end=None):
all_segments = [f'S{i:02}10' for i in range(1,11)]
......
......@@ -30,11 +30,11 @@ def get_sorting(grid_shape):
src_index_nn = np.memmap(index_dir / 'src_index_nn.dat', mode='r', dtype=np.uint32)
dst_index_nn = np.memmap(index_dir / 'dst_index_nn.dat', mode='r', dtype=np.uint32)
satzen = xr.open_dataset(SATZEN_CACHE / f'{sat}_satzen.nc')
nn_satzen = remap_fast_mean(src_index_nn, dst_index_nn, satzen.satzen.values, grid_shape)
nn_satzen = remap_fast_mean(src_index_nn, dst_index_nn, satzen.satellite_zenith.values, grid_shape)
satzens.values[i,...] = nn_satzen
wmo_ids.values[i,np.isfinite(nn_satzen)] = wmo_id
sample_mode.values[i, np.isfinite(nn_satzen)] = 1
ellip_satzen = remap_fast_mean(src_index, dst_index, satzen.satzen.values, grid_shape)
ellip_satzen = remap_fast_mean(src_index, dst_index, satzen.satellite_zenith.values, grid_shape)
mask = np.isfinite(ellip_satzen) & (ellip_satzen < max_satzen)
satzens.values[i,mask] = ellip_satzen[mask]
wmo_ids.values[i,mask] = wmo_id
......
......@@ -96,12 +96,20 @@ def make_geometry(dt_dir):
sat_azi = get_satazi(area)
ds = xr.Dataset()
ds['satzen'] = ['y','x'], sat_zen
ds['satellite_zenith'] = ['y','x'], sat_zen
out = SATZEN_CACHE / f'{sat}_satzen.nc'
bar.set_description(f'saving {out}')
if out.is_file():
out.unlink()
ds.to_netcdf(out, encoding=SATZEN_ENCODING)
ds = xr.Dataset()
ds['satellite_azimuth'] = ['y','x'], sat_azi
out = SATAZI_CACHE / f'{sat}_satazi.nc'
bar.set_description(f'saving {out}')
if out.is_file():
out.unlink()
ds.to_netcdf(out, encoding=SATAZI_ENCODING)
if __name__ == '__main__':
......
......@@ -9,6 +9,7 @@ from pathlib import Path
import warnings
from tqdm import tqdm
import subprocess
from collect_l1b import band_dir_path
def get_index_bands(res):
......@@ -59,7 +60,7 @@ def get_index_fast(area, pc, radius=2e3, nprocs=8):
grid_coords = pc.get_cartesian_coords()
grid_coords_pad = np.pad(grid_coords, ((0,0),(0,0),(0,1))).astype(np.float32)
grid_coords_pad.astype(np.float32).tofile('coord_descent/grid_coords.dat')
subprocess.run(['./main',str(rows),str(cols), str(grid_rows),str(grid_cols)], cwd='coord_descent', capture_output=False)
subprocess.run(['./main',str(rows),str(cols), str(grid_rows),str(grid_cols), f'{radius:.0f}'], cwd='coord_descent', capture_output=False)
sat_idx = np.memmap('coord_descent/src_index.dat', mode='r', dtype=np.uint32)
grid_idx = np.memmap('coord_descent/dst_index.dat', mode='r', dtype=np.uint32)
#s = pd.Series(sat_idx, index=grid_idx)
......@@ -120,7 +121,7 @@ def main(dt, r_sample=2):
for sat,band,res,r_footprint in bar:
prefix = f'{sat} {band}:'
bar.prefix = prefix
input_dir = Path(dt.strftime('l1b/%Y%m%dT%H%M')) / sat / band
input_dir = band_dir_path(dt, sat, band)
assert input_dir.exists(), str(input_dir)
input_files = list(input_dir.glob('*'))
output_dir = Path('index') / sat / band
......
......@@ -38,13 +38,13 @@ def default_attrs():
'satellite_names':';'.join(f'{v}={utils.SAT_NAMES[k]}' for k,v in utils.WMO_IDS.items())
}
attrs['satellite_zenith'] = {
attrs['satellite_zenith_angle'] = {
'long_name':'satellite zenith angle',
'standard_name':'satellite zenith angle',
'units':'degrees'
}
attrs['satellite_azimuth'] = {
attrs['satellite_azimuth_angle'] = {
'long_name':'satellite azimuth angle',
'standard_name':'satellite azimuth angle',
'units':'degrees'
......@@ -137,7 +137,7 @@ def rewrite_nc(f, out_root, dt, lat, lon):
def filename(k, dt):
#return f"{k}_{dt.strftime('%Y%m%dT%H%M')}.nc"
return f"ISCCP-NG_L1g_demo_A1_v1_res_0_10deg__{k}_{dt.strftime('%Y%m%dT%H%M')}.nc"
return f"ISCCP-NG_L1g_demo_A1_v1_res_0_05deg__{k}_{dt.strftime('%Y%m%dT%H%M')}.nc"
def set_latlon(ds, lat, lon):
......@@ -152,8 +152,14 @@ def set_latlon(ds, lat, lon):
def rewrite_nc_general(f, out_root, dt, lat, lon):
out_dir = out_root / dt.strftime('%Y') / dt.strftime('%Y%m') / dt.strftime('%Y%m%d') / dt.strftime('%Y%m%dT%H%M')
out_dir.mkdir(exist_ok=True, parents=True)
ds = xr.open_dataset(f)
k = next(iter(ds.data_vars))
out = out_dir / filename(k, dt)
if out.exists():
return out
grid_shape = ds[k].shape[-2:]
set_latlon(ds,lat,lon)
......@@ -164,14 +170,17 @@ def rewrite_nc_general(f, out_root, dt, lat, lon):
encoding = default_encoding(grid_shape)
ds = add_time(ds, dt, encoding)
out_dir = out_root / dt.strftime('%Y') / dt.strftime('%Y%m') / dt.strftime('%Y%m%d') / dt.strftime('%Y%m%dT%H%M')
out_dir.mkdir(exist_ok=True, parents=True)
out = out_dir / filename(k, dt)
ds.to_netcdf(out, encoding={k:v for k,v in encoding.items() if k in ds})
return out
def rewrite_wmo_id(f, out_root, dt, lat, lon):
out_dir = out_root / dt.strftime('%Y') / dt.strftime('%Y%m') / dt.strftime('%Y%m%d') / dt.strftime('%Y%m%dT%H%M')
out_dir.mkdir(exist_ok=True)
out = out_dir / filename('wmo_id', dt)
if out.exists():
return out
ds = xr.open_dataset(f)
grid_shape = ds['wmo_id'].shape[-2:]
......@@ -189,24 +198,27 @@ def rewrite_wmo_id(f, out_root, dt, lat, lon):
'complevel':1}
ds = add_time(ds, dt, encoding)
out_dir = out_root / dt.strftime('%Y') / dt.strftime('%Y%m') / dt.strftime('%Y%m%d') / dt.strftime('%Y%m%dT%H%M')
out_dir.mkdir(exist_ok=True)
out = out_dir / filename('wmo_id', dt)
ds.to_netcdf(out, encoding={k:v for k,v in encoding.items() if k in ds})
return out
def rewrite_satazi(f, out_root, dt, lat, lon):
out_dir = out_root / dt.strftime('%Y') / dt.strftime('%Y%m') / dt.strftime('%Y%m%d') / dt.strftime('%Y%m%dT%H%M')
out_dir.mkdir(exist_ok=True)
out = out_dir / filename('satellite_azimuth_angle',dt)
if out.exists():
return out
ds = xr.open_dataset(f)
grid_shape = ds['satellite_azimuth'].shape[-2:]
grid_shape = ds['satellite_azimuth_angle'].shape[-2:]
set_latlon(ds,lat,lon)
attrs = default_attrs()
ds['satellite_azimuth'].attrs.update(attrs.get('satellite_azimuth',{}))
ds['satellite_azimuth_angle'].attrs.update(attrs.get('satellite_azimuth_angle',{}))
encoding = default_encoding(grid_shape)
encoding['satellite_azimuth'] = {
encoding['satellite_azimuth_angle'] = {
'zlib':True,
'scale_factor':.125,
'dtype':'i2',
......@@ -217,25 +229,27 @@ def rewrite_satazi(f, out_root, dt, lat, lon):
}
ds = add_time(ds, dt, encoding)
out_dir = out_root / dt.strftime('%Y') / dt.strftime('%Y%m') / dt.strftime('%Y%m%d') / dt.strftime('%Y%m%dT%H%M')
out_dir.mkdir(exist_ok=True)
out = out_dir / filename('satellite_azimuth',dt)
ds.to_netcdf(out, encoding={k:v for k,v in encoding.items() if k in ds})
return out
def rewrite_satzen(f, out_root, dt, lat, lon):
out_dir = out_root / dt.strftime('%Y') / dt.strftime('%Y%m') / dt.strftime('%Y%m%d') / dt.strftime('%Y%m%dT%H%M')
out_dir.mkdir(exist_ok=True)
out = out_dir / filename('satellite_zenith_angle',dt)
if out.exists():
return out
ds = xr.open_dataset(f)
ds = ds.rename(satzen='satellite_zenith')
grid_shape = ds['satellite_zenith'].shape[-2:]
grid_shape = ds['satellite_zenith_angle'].shape[-2:]
set_latlon(ds,lat,lon)
attrs = default_attrs()
ds['satellite_zenith'].attrs.update(attrs.get('satellite_zenith',{}))
ds['satellite_zenith_angle'].attrs.update(attrs.get('satellite_zenith_angle',{}))
encoding = default_encoding(grid_shape)
encoding['satellite_zenith'] = {
encoding['satellite_zenith_angle'] = {
'zlib':True,
'scale_factor':.125,
'dtype':'i2',
......@@ -246,14 +260,17 @@ def rewrite_satzen(f, out_root, dt, lat, lon):
}
ds = add_time(ds, dt, encoding)
out_dir = out_root / dt.strftime('%Y') / dt.strftime('%Y%m') / dt.strftime('%Y%m%d') / dt.strftime('%Y%m%dT%H%M')
out_dir.mkdir(exist_ok=True)
out = out_dir / filename('satellite_zenith',dt)
ds.to_netcdf(out, encoding={k:v for k,v in encoding.items() if k in ds})
return out
def rewrite_pixel_time(f, out_root, dt, lat, lon):
out_dir = out_root / dt.strftime('%Y') / dt.strftime('%Y%m') / dt.strftime('%Y%m%d') / dt.strftime('%Y%m%dT%H%M')
out_dir.mkdir(exist_ok=True)
out = out_dir / filename('pixel_time',dt)
if out.exists():
return out
ds = xr.open_dataset(f)
grid_shape = ds['pixel_time'].shape[-2:]
......@@ -274,9 +291,6 @@ def rewrite_pixel_time(f, out_root, dt, lat, lon):
}
ds = add_time(ds, dt, encoding)
out_dir = out_root / dt.strftime('%Y') / dt.strftime('%Y%m') / dt.strftime('%Y%m%d') / dt.strftime('%Y%m%dT%H%M')
out_dir.mkdir(exist_ok=True)
out = out_dir / filename('pixel_time',dt)
ds.to_netcdf(out, encoding={k:v for k,v in encoding.items() if k in ds})
return out
......@@ -13,6 +13,7 @@ import satpy
import warnings
from pathlib import Path
from utils import spherical_angle_add, ALL_BANDS, AHI_BANDS, ABI_BANDS, MSG_BANDS, remap_fast_rad_mean, remap_fast_mean, remap_with_stats_rad, remap_with_stats, WMO_IDS, ALL_SATS, STATS_BANDS, STATS_FUNCS, BAND_CENTRAL_WAV
WMO_ID = WMO_IDS
from make_index import get_index_bands
from collect_l1b import L1B_DIR
......@@ -21,6 +22,8 @@ import tempfile
import os
import shutil
import sys
COMP_CACHE = Path('composite_cache')
COMP_CACHE.mkdir(exist_ok=True)
......@@ -34,6 +37,11 @@ ENCODING = {
'add_offset':50
}
WMO_IDS = xr.open_dataset(COMP_CACHE / 'wmo_id.nc').wmo_id
SAMPLE_MODE = xr.open_dataset(COMP_CACHE / 'sample_mode.nc').sample_mode
GRID_SHAPE = WMO_IDS.shape
print(GRID_SHAPE)
orig_print = print
def print(*args, flush=False, **kwargs):
orig_print(*args, flush=True, **kwargs)
......@@ -77,7 +85,7 @@ def read_scene(files, reader, bar=None):
try:
area = scene[ds_names[0]].area
except KeyError as e:
print(e.args)
print('Error', e.args)
raise IOError('Problem reading files')
if bar is not None:
bar.set_description(dt.strftime(f'Loading {sat} band {band} %Y%m%dT%H%M'))
......@@ -87,15 +95,14 @@ def read_scene(files, reader, bar=None):
def composite_band(composite, band, index_band, sat, dt, reader, wmo_id, wmo_ids, sample_mode, with_stats=False, bar=None):
wavelength = BAND_CENTRAL_WAV[band]
grid_shape = wmo_ids.shape
if composite is None:
if with_stats:
composite = {k:
xr.DataArray(np.full(grid_shape, np.nan, dtype=np.float32), dims=['layer','latitude','longitude'])
xr.DataArray(np.full(GRID_SHAPE, np.nan, dtype=np.float32), dims=['layer','latitude','longitude'])
for k in STATS_FUNCS}
else:
composite = xr.DataArray(np.full(grid_shape, np.nan, dtype=np.float32), dims=['layer','latitude','longitude'])
band_dir_path(L1B_DIR, dt, sat, band)
composite = xr.DataArray(np.full(GRID_SHAPE, np.nan, dtype=np.float32), dims=['layer','latitude','longitude'])
band_dir = band_dir_path(L1B_DIR, dt, sat, band)
src_index, dst_index, src_index_nn, dst_index_nn = open_index(INDEX, sat, index_band)
files = list(band_dir.glob('*'))
v, area = read_scene(files, reader)
......@@ -117,9 +124,9 @@ def composite_band(composite, band, index_band, sat, dt, reader, wmo_id, wmo_ids
else:
remap = remap_with_stats
out = remap(src_index, dst_index, v, grid_shape[-2:])
out_nn = remap(src_index_nn, dst_index_nn, v, grid_shape[-2:])
scene.unload()
out = remap(src_index, dst_index, v, GRID_SHAPE[-2:])
out_nn = remap(src_index_nn, dst_index_nn, v, GRID_SHAPE[-2:])
#scene.unload()
def do_composite(composite, out_nn, out, do_nn=True):
if bar is not None:
......@@ -142,10 +149,6 @@ def composite_band(composite, band, index_band, sat, dt, reader, wmo_id, wmo_ids
def main(dt, progress=True):
wmo_ids = xr.open_dataset(COMP_CACHE / 'wmo_id.nc').wmo_id
sample_mode = xr.open_dataset(COMP_CACHE / 'sample_mode.nc').sample_mode
grid_shape = wmo_ids.shape
print(grid_shape)
ordered_bands = ['temp_11_00um', *sorted(ALL_BANDS - set(['temp_11_00um']))]
out_dir = COMP_CACHE / dt.strftime('%Y') / dt.strftime('%Y%m') / dt.strftime('%Y%m%d') / dt.strftime('%Y%m%dT%H%M')
out_dir.mkdir(exist_ok=True, parents=True)
......@@ -159,7 +162,7 @@ def main(dt, progress=True):
else:
out_nc = out_dir / f'{band}.nc'
if out_nc.exists():
#print(f'Already have {out_nc}')
print(f'Already have {out_nc}')
continue
start = time.time()
def run(it,bar=None):
......@@ -170,7 +173,7 @@ def main(dt, progress=True):
if band not in attrs['bands']:
continue
res = attrs['res'][band]
wmo_id = WMO_IDS[sat]
wmo_id = WMO_ID[sat]
index_band = get_index_bands(attrs['res'])[res]
tmp_root = Path(tempfile.gettempdir())
tmp = tmp_root / dt.strftime(f'{sat}_{band}_%Y%m%dT%H%M')
......@@ -178,7 +181,7 @@ def main(dt, progress=True):
try:
tempfile.tempdir = str(tmp)
os.environ['TMP'] = str(tmp)
composite = composite_band(composite, band, index_band, sat, dt, reader, wmo_id, wmo_ids, sample_mode, with_stats=with_stats, bar=bar)
composite = composite_band(composite, band, index_band, sat, dt, reader, wmo_id, WMO_IDS, SAMPLE_MODE, with_stats=with_stats, bar=bar)
except IOError:
print(f'Error reading {sat}')
finally:
......@@ -225,6 +228,7 @@ if __name__ == '__main__':
if args.end is not None:
end = pd.to_datetime(args.end)
for dt in pd.date_range(dt, end, freq='30min'):
print(dt)
try:
main(dt)
except Exception as e:
......
from pathlib import Path
import os
os.environ['XRIT_DECOMPRESS_PATH'] = str(Path('xrit/PublicDecompWT/xRITDecompress/xRITDecompress').absolute())
import timing
import netCDF4
import xarray as xr
import numpy as np
from tqdm import tqdm
from datetime import datetime, timedelta
from utils import ALL_SATS, remap_fast_mean
from make_index import get_index_bands
from make_sample import open_index, band_dir_path, read_scene
import satpy
COMP_CACHE = Path('composite_cache/')
INDEX = Path('index')
L1B_DIR = Path('l1b')
ABI_SCAN_DIR = Path('ancil/abi_scan_schedule/')
WMO_IDS = xr.open_dataset(COMP_CACHE / 'wmo_id.nc').wmo_id
SAMPLE_MODE = xr.open_dataset(COMP_CACHE / 'sample_mode.nc').sample_mode
GRID_SHAPE = WMO_IDS.shape
def saveit(composite, out):
fill = netCDF4.default_fillvals['i2']
ds = xr.Dataset()
ds['pixel_time'] = composite.fillna(fill).astype(np.int16)
encoding = {'pixel_time':{'zlib':True,'chunksizes':(1, 1800, 3600), '_FillValue':fill, 'dtype':'i2'}}
ds.to_netcdf(out, encoding=encoding)
def run_one(dt):
out_dir = COMP_CACHE / dt.strftime('%Y') / dt.strftime('%Y%m') / dt.strftime('%Y%m%d') / dt.strftime('%Y%m%dT%H%M')
out_path = out_dir / 'pixel_time.nc'
if out_path.exists():
return
composite = xr.DataArray(np.full(GRID_SHAPE, np.nan, dtype=np.float32), dims=['layer','latitude','longitude'])
for attrs in ALL_SATS[:]:
prefix = (attrs['name'])
_,index_band = max(get_index_bands(attrs['res']).items())
src_index, dst_index, src_index_nn, dst_index_nn = open_index(INDEX, attrs['sat'], index_band)
band_dir = band_dir_path(L1B_DIR, dt, attrs['sat'], 'temp_11_00um')
print(band_dir)
files = list(band_dir.glob('*'))
try:
v, area = read_scene(files, attrs['reader'])
if attrs['reader'] == 'seviri_l1b_hrit':
start_time, line_times = timing.meteosat_get_time_offset(v)
offsets = timing.meteosat_estimate_pixel_time_offsets(line_times)
elif attrs['reader'] == 'ahi_hsd':
start_time, line_times = timing.himawari_line_times(files)
offsets = timing.himawari_estimate_pixel_time_offsets(line_times)
elif attrs['reader'] == 'abi_l1b':
offsets = timing.goes_pixel_time_offset(ABI_SCAN_DIR)
start_time = timing.goes_start_time(files)
adjust = (start_time - dt).total_seconds()
offsets += adjust
out_nn = remap_fast_mean(src_index_nn, dst_index_nn, offsets, GRID_SHAPE[-2:])
for layer in range(composite.shape[0]):
mask = (WMO_IDS[layer].values == attrs['wmo_id'])
composite.values[layer, mask] = out_nn[mask]
except Exception as e:
print(f'Problem reading {attrs["sat"]}')
saveit(composite, out_path)
def get_timing_list():
dts = sorted([datetime.strptime(i.name,'%Y%m%dT%H%M') for i in COMP_CACHE.glob('*/*/*/*')])
with open('date_list.txt','w') as fp:
for dt in dts:
fp.write(dt.strftime('%Y%m%dT%H%M\n'))
def main(task_id, num_tasks):
with open('date_list.txt') as fp:
dts = [datetime.strptime(i.strip(),'%Y%m%dT%H%M') for i in fp]
dts = dts[task_id::num_tasks]
print(f'{len(dts)} tasks')
for i,dt in enumerate(dts,1):
print(f'{i}/{len(dts)}', flush=True)
try:
run_one(dt)
except IOError:
print('problem reading',dt, flush=True)
print(flush=True)
if __name__ == '__main__':
import argparse
parser = argparse.ArgumentParser()
parser.add_argument('task_id',type=int)
parser.add_argument('max_task_id',type=int)
args = parser.parse_args()
main(args.task_id, args.max_task_id+1)
%% Cell type:code id:11319c50 tags:
``` python
```
%% Cell type:code id: tags:
``` python
%pylab inline
```
%% Cell type:code id: tags:
``` python
import xarray as xr
import numpy as np
import pandas as pd
from netCDF4 import Dataset
from tqdm import tqdm
from pathlib import Path
import shutil
from datetime import datetime
```
%% Cell type:code id: tags:
``` python
ROOT = Path('../final')
```
%% Cell type:code id: tags:
``` python
dirs = sorted(ROOT.glob('*/*/*/*'))
```
%% Cell type:code id: tags:
``` python
f = next(dirs[0].glob('*temp_11_00um_20*'))
```
%% Cell type:code id: tags:
``` python
ds = xr.open_dataset(f)
```
%% Cell type:code id: tags:
``` python
f = '/ships19/cloud/scratch/cphillips/abi_scanline_time_luts/ABI-Time_Model_LUTS/ABI-Timeline05B_Mode 6A_20190612-183017.nc'
```
%% Cell type:code id: tags:
``` python
ds = xr.open_dataset(f, mask_and_scale=False, decode_times=False)
```
%% Cell type:code id: tags:
``` python
ds.FD_pixel_times[::8,::8].plot.imshow()
```
%%%% Output: execute_result
<matplotlib.image.AxesImage at 0x7fae930d5d90>
%%%% Output: display_data
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAEHCAYAAACA3BA3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAjmElEQVR4nO3de5hkdX3n8fenm7tCAAdwwgyCiLpouIw8QIIaDCqoCCRZsyRrQgj7zJpFxWfdCBNy0+ch0VwIipd1VlCIGIIXZNYlIpIQL+GO3AYkjAziyMgwCIKgA9P92T/Or2eKme6ec7qrqk/VfF7Pc56qc+r8Tn2rp6e+fX5X2SYiIqKpkbkOICIiBlMSSEREzEgSSEREzEgSSEREzEgSSEREzMg2cx3AbM2bN88v2mefqU/w+Ayu2rBn2kx6so2PNS6ihp/FY+ubvcFY85iav0fD84Hx9c3jGh9r+LOayXs0LDM+1vz3ZHx9s8/R9HyAsYZxrZ/B73vTMutn8F+qTpG1PLPW9h7Nr77RQu3on1Pv57yWZ66yfdxs3q/NBj6BvGifffj2t7455et69ufNLzre7EtOz/6s8VuMrHtqBmWebHT+2E8ebXT++JOPNTp/Ju8x9viPG7/Hz3/8RPMyjzX7Wf380ebv8bNHf9rsPR5r/rv49Nqnm53/aPPfxcd++kyz859tnqQeWdfs/9Qj62aQ0Guc80m+//3GF97EzxnnN5lf69xP8v15s32/Nhv4BBIRszMqzXUIA0XAaN0f2ZAPs0sCiYhoQMB2IzUzSPMbqYGSBBIR0UB1B5K7NkgCiQGgkXQWjBZRgyqsIZcEEhHRQO5ANkoCiYhooFEj+pBLAomIaES5AymSQCIiGhCwbRIIkAQSEdGI0oi+QRJIRERDqcKqJIFEROvM5C/88T6N+k4j+kbDn0BmMoZgJvMvRsRWId14Nxr+BBIR0WW5A6kkgUTM0Ei+RbZKUoO5sIZcEkhERANpA9koCSQiogFlIOEGSSDDLJMQRvRE7kAqSSAREQ1UAwmTQSAJJCKikUYLSg25JJCIiAbSiL5RzyvJJe0q6QuSvivpHkm/LGl3SVdLuq887tZx/hJJKyTdK+nYXscXEdHUqFRrG3b9aGX9MPBV2y8HDgbuAc4CrrF9AHBN2UfSgcDJwCuA44CPSxrtQ4wREbVIMCLV2oZdTxOIpF2A1wIXANh+xvbjwInAReW0i4CTyvMTgUttr7O9ElgBHN7LGCMimhEarbcNu163gbwYeAT4tKSDgVuAM4C9bK8GsL1a0p7l/L2B6zvKryrHnkPSYmAxwMKFC3sXfcSAURp3e0+ZhWBCr6uwtgEWAZ+wfSjwFKW6agqT/atsNsem7aW2D7N92B7z5nUn0oiIGiQY2W601jbsep1AVgGrbN9Q9r9AlVAeljQfoDyu6Ti/85ZiAfBQj2OcPY003yJiMEmMjNbbhl1Pv8ls/wj4gaSXlUPHAHcDy4BTyrFTgCvK82XAyZK2l7QfcABwYy9jjP7S6EjjLaJtNDJSaxt2/RgH8i7gEknbAfcDp1IlrssknQY8CLwNwPZySZdRJZn1wOm2x/oQY0RELUobyAY9TyC2bwMOm+SlY6Y4/xzgnF7GFBExG93sYVWGKtwM/ND28ZJ2B/4J2Bd4APgt24+Vc5cApwFjwLttX9W1QGZg+O+xZiLtGRExFanb1a9nUI2PmzAw4+Ty7RcR0YAEo9uO1Nq2fC0tAN4CfKrj8MCMk8tcWBERjYiR+ncX8yTd3LG/1PbSjv3zgPcBO3ccm9U4uX5KAomIaEKN2kDW2p6sDRhJxwNrbN8i6eh677yZzcbJ9VMSSET03DBNLKhmCWQ6RwEnSHozsAOwi6TPUsbJlbuPVo+TSxtIRERDI6Mjtbbp2F5ie4Htfakax//F9tsZoHFyuQOJiGhAUq0G8ln4IAMyTi4JJAIy4j3qU/d/X2xfC1xbnj/KgIyTSwKJiGgoI9ErSSAREU1o61jro44kkAHiIerJEjGo1IMqrEGVBBIR0VCqsCpJIBERDVS9sIZ/sag6kkAiIpro3kDCgZcEEhHRRNpANkgCiYhoRFvFaoN1JIFERDRQrUg42AmkLFq1JeO2H5/uhCSQiBZLXXsbaRiqsB4q23S/YKPAPtNdJAkkIqIJiZFtB/6r8x7bh053gqTvbOkiA/9TiIjoK4FGB74b7y9345wkkJgVjQz8f6SIRjQEVVi2fw4gaX9gle11ZVGrg4CLbT8+cc50BvunEBHRb4KRkZFa2wD4IjAm6SXABcB+wOfqFs4dSEREQ4N+B9Jh3PZ6Sb8OnGf7/DptHxN6/lOQ9ICkOyXdNrG4vKTdJV0t6b7yuFvH+UskrZB0r6Rjex1fREQTKo3odbYB8Kyk36Za+fAr5di2dQv3K42+zvYhHYvLnwVcY/sA4Jqyj6QDqZZ2fAVwHPBxSalkj4j2UNUGUmcbAKdSNZafY3tlWSr3s3ULz1WKPBE4ujy/iGolrjPL8UttrwNWSloBHA5cNwcxRkRsboimMrF9N/Dujv2VVEvq1tKPBGLga5IMfNL2UmAv26sBbK+WtGc5d2/g+o6yq8qx55C0GFgMsHDhwl7GHhGxmSEYiX4n1XfzpGwfVOc6/UggR9l+qCSJqyV9d5pzJxsVudmHLEloKcCrFi2a8ocQEd23tQ+Ol4ZiLqzjy+Pp5fEfyuN/BZ6ue5GeJxDbD5XHNZIup6qSeljS/HL3MR9YU05fBXTeUiygGm4fEdEag16FZfv7AJKOsn1Ux0tnSfo28IE61+npT0HS8yTtPPEceCNwF7CMqtWf8nhFeb4MOFnS9qUx5wDgxl7GGBHRyHD1wnqepFdP7Ej6FeB5dQv3+hPuBVyuai3vbYDP2f6qpJuAyySdBjwIvA3A9nJJlwF3A+uB022P9TjGGEKDXkcd7SWJkcGfymTCacCFkn6h7D8O/EHdwj1NILbvBw6e5PijwDFTlDkHOKeXcUVEzMagV2FNsH0LcLCkXQDZ/kmT8gNxjxUR0RpD1I1X0vbAbwL7AtuU2iJs12oDSQKZKxqOX8CIrc9Q9MKacAXwE+AWYF3TwkkgERENaIjuQIAFto+baeEkkIiIJoZjQakJ/y7pl2zfOZPCw/FTSHVQRPSLNEzr4Lwa+H1JK6mqsAS4TSPRIyKGSxcSiKSFwMXAC4FxYKntD0vaHfgnqobtB4Dfsv1YKbOEquvtGPBu21fNMow3zaZw/nSPiGhEMDJSb5veeuC9tv8TcCRwepmRvG+zlZcR6bsCby3brhOj1OtIAomIaKKsiV5nm47t1bZvLc+fBO6hmjz2RKpZyimPJ5XnG2YrL7PmTsxWPvOPIp0BXALsWbbPSnpX3fKpworoE23tsxAOCwm22a7u2fMmFtIrlpbJYDe5pPYFDgVuYJazlTd0GnCE7adKHB+iWj7j/DqFk0AiIhpQs3EgazsW0pv8etLzqdYmf4/tJyYG80361pub7WzkompPmTA2xftMKgkkIqIJ0ZVGdABJ21Ilj0tsf6kc7uds5Z8GbigzpUNVXXZB3cJpA4mNRkabbxFbHXXl/4eqW40LgHtsn9vxUt9mKy/veyrwY+Ax4FTb59UtnzuQiIiGujSVyVHA7wJ3SrqtHPtjqiVl+zJbuaQjgeUTjfmSdpZ0hO0b6pRPAomIaELqyt237W8xdXtDv2Yr/wSwqGP/qUmOTSkJJCKiCQltW7sXVtvJ9oaGeNvjkmrnhbSBRMRQGKmxdUfXBhK2wf2S3i1p27KdAdxft/BAfMKIiNaY6IU1HJ1M3gH8CvBDql5eRwCL6xZOFVZERCPDM5mi7TVU06PMSO5AIrZyo2q2BUNThSXppZKukXRX2T9I0p/ULd/+TxgR0SYaQdtsV2sbAP8HWAI8C2D7DhrckaQKKyKiCTEQdxc17WT7xk2mT1lft3ASSLTeEC0fGkNAaIsz7Q6QtZL2p8ypJek/A6vrFk4CiYhoootzYbXA6cBS4OWSfgisBN5et3Bf/rSTNCrpO5K+UvZ3l3S1pPvK424d5y6RtELSvZKO7Ud8ERH1dWcurDawfb/t1wN7AC+3/WrbD9Qt36+6gTOoFkuZ0LcVtyIiuk0jI7W2tpN0hqRdgKeBv5d0q6Q31i3f808oaQHwFuBTHYf7tuJWRERXSTC6bb2t/f7A9hPAG6lWJDyVajLHWvqRIs8D3ke1aPyE56y4RRU4VKtr/aDjvElX3JK0WNLNkm5+ZO3angQdETE5gUbqbe030f3qzcCnbd/ecWyLevoJJR0PrLF9S90ikxzbbMUt20ttH2b7sD3mzZtVjBERTVkjtbYBcIukr1ElkKsk7cxz/9ifVq97YR0FnCDpzcAOwC6SPkt/V9yKiAEzOvWyrlMa82xXd61JDMrdRR2nAYcA99t+WtILqKqxaunpT8H2EtsLbO9L1Tj+L7bfTj9X3Kr5l8KA/dUQEXNGVTtIna2lJL0Qqunbbd9q+/Gy/2gZjb7hnOnU/raUdJSk55Xnb5d0rqQXzSj6qpHmDZLuA95Q9rG9HJhYceurdGHFrYiIrhv8ubCu7MY5TaqwPgEcLOlgqkbxC4CLgV+tU9j2tcC15fmj9G/FrYieUGYW3CpZwiMDPwb7YElPTPO6gOleB5olkPW2LelE4MO2L5B0yhZLRUQMmwGv6rbdlfF1TRLIk5KWUA1zf20Z4DcQHZ0jIrpHA59AuqVJAvkvwO8Ap9n+kaR9gL/pTVjRDcOy6E1E26SzTaV2ArH9I+Dcjv0HqdpAIiK2LkkgQI0EIulJnjuYT2VfgG3v0qPYIiLaRxqIiRLrKFO5r7K9TtLRwEHAxRPderdki2nU9s62d+nYdu58nE3wERGDaIjGlH0RGJP0EqqetfsBn6tbuFFftNKF9zVl9xsTA04iIrYeavsYjybGba+X9OvAebbPl/SduoWbDCQ8A7iEauLDPYFLJL2rcbgREYNsYiqT4ZhM8VlJv001I8hXyrHavWub3IGcBhxh+ykASR8CrgPOb3CNiIgBN1TdeE8F3gGcY3tlmULqs3ULN0kgAjqnFRmjwbS/EdF7I1k/vj+GJIHYvlvSmcA+ZX8lDdYDaZJAPg3cIOnysn8SVaNLDMkvU0TUMBxTmQAg6a3A3wLbAftJOgT4gO0T6pSv/c1n+1yq250fA48Bp9o+r2nAEREDr0uz8Uo6TtK9klZIOqsPkW/qL6hWfX0cwPZtVD2xaqmdRiUdCSy3fWvZ31nSEbZvaBBsbOU0Ohz952Nr1p02kDId1MeoZiRfBdwkaZntu2d98frW2/6Jnpvsai+s0uSn8Angpx37T5VjERFblS6NAzkcWGH7ftvPAJcCJ/Y8+Oe6S9LvAKOSDpB0PvDvdQs3SSCyNy75ZXuc3q9oGBHRPvW78c6TdHPHtrjjKnsDP+jYX1WO9dO7gFcA64B/pJrC/T11CzdJAPdLejcb7zr+B3B/g/KDI43iETEFS7j+aoNrbR82xWuTXaRP6/KWN7OfBs4uW2NNvinfAfwK8EOqTHkEsHjaEhERw8ZmbLzetgWrgIUd+wuAh3oW9yQkHSbpS5JulXTHxFa3fJPZeNdQrWs+VSBLbP9V3etFRAyqLt0m3AQcUAbv/ZDq+/V3unPp2i4B/gi4ExhvWribbRhvA5JAImKoGdjyzUWN61RzUL0TuAoYBS60vXz2V27kEdvLZlq4mwkko9IjYlLDtnx8R3+i2V7nSuDKrlxsZv5c0qeAa6ga0gGw/aU6hbuZQPra+BMRMRe6dQfSEqcCL6eaQHGiCstA3xPIkP2NERExCcPY8CSQg23/0kwLdzOBfL6L14roq0xCGE10qwqrBa6XdOBMR7/X+l8j6XWlq9fysn2hLH+4ge2/nKTcDpJulHR7Kff+cnx3SVdLuq887tZRZkmZF+ZeScfO5ENFRPSKqep66mwD4NXAbeX79g5Jd3a1G6+ktwAfBT4AvJ+qqmoRcKGkd5ZGoKmsA37N9k8lbQt8S9I/A78BXGP7g2UCsbOAMyUdSNWV7RXALwJfl/RS22NTvUFERL8Nzw0Ix82mcJ0qrD8CTrJ9e8ex2yTdTLWY1JQJpEx9MjF/1rZlM9V8L0eX4xcB1wJnluOX2l4HrJS0gmq+mOtqfp7hlhHyEa0w6I3oknax/QTw5GyuUyeBvHCT5AGA7Tsk7bWlwmXGyVuAlwAfs32DpL1sry7XWS1pz3L63sD1HcXnYm6YiIgp2UPRBvI54Hiq72bz3E5QBl5c5yJ1EshTM3ytiqSqfjpE0q7A5ZJeOc3pteaGKROSLQZYuHDhZgUiInpp0Hth2T6+PNZe+2MydRLI/pImG6koamYpANuPS7qWqs7tYUnzy93HfGBNOa3W3DC2lwJLAV61aNGA/1NGxCCpxoEMx9eOpGtsH7OlY1Opk0Cmm5/+b7cQ3B7AsyV57Ai8HvgQsAw4hWrt3VOAK0qRZcDnJJ1L1Yh+AHBjjRgjIvpm0NOHpB2Anaimm9+NjbU/u1B999ZSJ4GstP1g8xABmA9cVNpBRoDLbH9F0nXAZZJOAx6kmkcL28slXQbcDawHTk8PrIhom0FvRAf+O9W6H79I1Q4ykUCeoFolsZY6CeTLVN12kfRF279Z9+K27wAOneT4o8Ckt0i2zwHOqfseEcMsAxzbadBrsGx/GPiwpHfZPn+q8yS9wfbVU71e57ezs2G7dptHRMQwMma85tZ20yWP4kPTvVjnDsRTPI+I2PoYxgZkmHkXTDvHYZ0EcrCkJ8qFdizPJy5s27vMMsCIiIFhBr8Kq4FpP+kWE4jt0e7FEkNnJL8esfUZhOqpfujmbLwREVuFregO5IHpXkwCiYhoYBgGEkr6jelen1iR0Pa05yWBREQ0YMOzgz6XCbx1mtfmZEXCiIitgBkb8DsQ26d24zoZpRQR0cBEFVadre0k7SXpgrJOE5IOLDOE1JIEEhGNjEqNt6FSxoHU2QbAZ4Cr2Dj/1X9QTXFSSxJIREQDw3QHAsyzfRllBV7b64Ha8w+mDSQioqFBbwPp8JSkF1AGDEo6EvhJ3cJJIBERDQxJL6wJ76VaRmN/Sd8G9qDMjl5HEkgMJY1mhHz0hhmY6qktsn2LpF8FXkY1PdW9bGH+q05pA4mIaGjM9ba2K6vELrC93PZdwCHATXXL5w5kjngGPVOGrC9LxEAahpHoHf4K+KqkjwB7A28Gao8RSQKJiGjCMDYESxIC2L5K0juAq4G1wKG2f1S3fKqwIiIaMPDsuGttsyHpbyR9V9Idki6XtGvHa0skrZB0r6RjO46/StKd5bWPSNNXdUj6U+B84LXAXwDXSnpL3RgHPoF4C1tERDf1cRzI1cArbR9ENcBvCVSjxYGTgVcAxwEflzTRa+QTwGLggLIdt4X3mAccbvs6258EjiUDCSMiesRmfLzeNru38dfKwD6A64EF5fmJwKW219leCawADpc0H9ilJAMDFwMnbeE9zrD9s47979t+Q90Y0wYSEdGAadTDap6kmzv2l9peOoO3/QPgn8rzvakSyoRV5diz5fmmxzcj6Tzb75H0f5mkssb2CXWCSgKJiGioQfXUWtuHTfWipK8DL5zkpbNtX1HOORtYD1wyUWyS8z3N8cn8Q3n826liqyMJJKJPRkbTEXsYVHcg3Wlhtf366V6XdApwPHBMqZaC6s5iYcdpC4CHyvEFkxyf7H1vKY//Jmk74OVUH+1e28/Ujb+nbSCSFkr6V0n3SFou6YxyfHdJV0u6rzzu1lFm0t4FERFtMDGVSZ1tNiQdB5wJnGD76Y6XlgEnS9pe0n5UjeU32l4NPCnpyNL76veAK7bwHm8Bvgd8BPgosELSm+rG2Os7kPXAe23fKmln4BZJVwO/D1xj+4OSzgLOAs7cpHfBLwJfl/RS27Vnh4yI6LU+DST8KLA9cHXpjXu97XfYXi7pMuBuqu/Y0zu+I/+Qaor2HYF/Ltt0/g54ne0VAJL2B/5fjXJAjxNIyYiry/MnJd1D1ahzInB0Oe0i4FqqTLuhdwGwUtIK4HDgul7GGRWNZP6oiC1xn1YktP2SaV47BzhnkuM3A69s8DZrJpJHcT+wpm7hvrWBSNoXOBS4AdirJBdsr5a0Zzltqt4FERHtMEQj0YHlkq4ELqNqA3kbcJOk3wCwPe3a6H1JIJKeD3wReI/tJ6YZHFmrF4GkxVSDZVi4cOFmBSIiesUMVQLZAXgY+NWy/wiwO/BWqo86twlE0rZUyeOSjmz2sKT55e5jPhtvmabqXfAcpR/1UoBFixYNzb9kRLSfh+gOxPa0EydKWmL7r6Z6vde9sARcANxj+9yOl5YBp5Tnp7Cxp8CkvQt6GWNERBPGPLN+vNY2BKZdXKrXdyBHAb8L3CnptnLsj4EPApdJOg14kBLkFnoXRERMqdYwm27cOAzRHUgN0/5Ue90L61vTBHDMFGUm7V0QEdEGQ9YGsiXTftCMRI/Yyo3OYHGzrdkwtYHUMHd3IBERw2grSiCfn+7FJJCIiAbGbdYNQQO5pNcB7wJeVg7dA3zU9rUT59j+y+mukQQSrafRLFsT7TLodyBlDqyPAh8A3k9VVbUIuFDSO21fWec6SSAREQ0MSRvIHwEn2b6949htZe2S84EkkIiIXujHXFg99sJNkgcAtu+QtFfdiySBREQ0YDwMdyBPzfC150gCiYhoYEiqsPaXtGyS4wJeXPciSSAREQ0YeGb9wE+QceI0r9Ve5jYJJCKiCQ9FFdZK2w/O9iLpHxkRrTOq5lu/TExlUmdrsS9PPJH0xZleZPjvQJQcGRHdY8P6dieHOjpTbu02j00NfwKJiOiiIZlM0VM8byQJJGKGRjJCfus0HL2wDpb0BNWdyI7lOWXftnepc5EkkIiIBiYWlBpktke7cZ0kkIiIBoZkHEhXJIEMs5Gu/JEREZtwEgiQBBIR0YgN40kgQBJIRERDxoM/mWJXJIFERDRhGBvwRvRuSQKJiGjAgJM/gCSQiIjGUoVVSQKJiGgijegb9HQoraQLJa2RdFfHsd0lXS3pvvK4W8drSyStkHSvpGN7GVvEMNKoGm/RlPF4va0bJP0vSZY0r+PYpN+Vkl4l6c7y2kck9fQfuNdzMXwGOG6TY2cB19g+ALim7CPpQOBk4BWlzMclDcZABo003yJiIFVtIP1JIJIWAm8AHuw4Nt135SeAxcABZdv0+7erevpNZvsbwI83OXwicFF5fhFwUsfxS22vs70SWAEc3sv4IiIaM4yNjdfauuDvgffx3AkPJ/2ulDQf2MX2da4aaS5m4/drT8xFG8hetlcD2F4tac9yfG/g+o7zVpVjMUwyOj6GQIO7i3mSbu7YX2p7aZ2Ckk4Afmj79k1qoqb6rny2PN/0eM+0qRF9srq6Sf+VJC2muk1j4cKFvYwpIuI5bDdpRF9r+7CpXpT0deCFk7x0NvDHwBsnKzZZWNMc75m5SCAPS5pf7j7mA2vK8VVAZzZYADw02QVKBl8KsGjRonSHiIi+6lY3Xtuvn+y4pF8C9gMm7j4WALdKOpypvytXleebHu+ZuWjNXQacUp6fAlzRcfxkSdtL2o+qAejGOYgPRkaabRGxVfF4vW3G17fvtL2n7X1t70uVHBbZ/hFTfFeWpoEnJR1Zel/9Hhu/X3uip3cgkv4ROJqqHnAV8OfAB4HLJJ1G1bPgbQC2l0u6DLgbWA+cbnusl/FFRDQ115MpbuG78g+per/uCPxz2XqmpwnE9m9P8dIxU5x/DnBO7yKKiJglw3if58IqdyGd+5N+V9q+GXhln8JqVSN6RMQAMOOZygRIAomIPhjt7YDovpoYSBhJIBERzTgJZEISSEREQ5lMsZIEEgEo3bGjJtuMd2eakoGXBBIR0VDuQCpJIBERDXk8Q9QgCWSwZBr4iLlnJ4EUSSAREQ2YJJAJSSAREU0YPJYEAkkgERHNeJzx9c/MdRStkAQSEdFQqrAqSSAREQ2kDWSjJJCIiCacO5AJSSARLabR4ZmEcHiY8SQQIAkkIqKZjAPZIAkkIqIB24w/m15YkAQSs6TR0bkOIaLvcgdSSQKJiGgiVVgbJIFERDSSBDIhCSQiooFqSdusBwJJIBERzdiZyqRIAomIaMIZBzKhdQtMSDpO0r2SVkg6a67jiYjoZKrZeOtsw65VdyCSRoGPAW8AVgE3SVpm++65jSwiJrTur85+Sy+sDVqVQIDDgRW27weQdClwIpAEEhEtkQQyoW0JZG/gBx37q4Aj5iiWiIhJJYFU2pZAJps5zpudJC0GFpfddTvttNNdPY1qduYBa+c6iCm0OTZod3xtjg3aHd9cxvai2V7AT6+96plbPzWv5ult/TfoirYlkFXAwo79BcBDm55keymwFEDSzbYP6094zbU5vjbHBu2Or82xQbvja3Nsddg+bq5jaIu2tYfdBBwgaT9J2wEnA8vmOKaIiJhEq+5AbK+X9E7gKmAUuND28jkOKyIiJtGqBAJg+0rgygZFlvYqli5pc3xtjg3aHV+bY4N2x9fm2KIB2Zu1UUdERGxR29pAIiJiQCSBRETEjAx0ApnrebMkXShpjaS7Oo7tLulqSfeVx906XltSYr1X0rF9iG+hpH+VdI+k5ZLOaEuMknaQdKOk20ts729LbB3vNyrpO5K+0sLYHpB0p6TbJN3cwvh2lfQFSd8tv3+/3Kb4oktsD+RG1Uvre8CLge2A24ED+xzDa4FFwF0dx/4aOKs8Pwv4UHl+YIlxe2C/Evtoj+ObDywqz3cG/qPEMecxUg0afX55vi1wA3BkG2LriPF/Ap8DvtLCf9sHgHmbHGtTfBcB/6083w7YtU3xZevONsh3IBvmzbL9DDAxb1bf2P4G8ONNDp9I9Z+H8nhSx/FLba+zvRJYQfUZehnfatu3ludPAvdQTRcz5zG68tOyu23Z3IbYACQtAN4CfKrjcCtim0Yr4pO0C9UfVxcA2H7G9uNtiS+6Z5ATyGTzZu09R7F02sv2aqi+wIE9y/E5jVfSvsChVH/ptyLGUkV0G7AGuNp2a2IDzgPeB3QuPdeW2KBKtl+TdEuZ2qdN8b0YeAT4dKkC/JSk57UovuiSQU4gtebNapE5i1fS84EvAu+x/cR0p05yrGcx2h6zfQjVlDWHS3rlNKf3LTZJxwNrbN9St8gkx3r9b3uU7UXAm4DTJb12mnP7Hd82VFW7n7B9KPAUVZXVVAbt/3IUg5xAas2bNQceljQfoDyuKcfnJF5J21Ilj0tsf6mNMZbqjWuB41oS21HACZIeoKoa/TVJn21JbADYfqg8rgEup6ryaUt8q4BV5Y4S4AtUCaUt8UWXDHICaeu8WcuAU8rzU4ArOo6fLGl7SfsBBwA39jIQSaKqh77H9rltilHSHpJ2Lc93BF4PfLcNsdleYnuB7X2pfq/+xfbb2xAbgKTnSdp54jnwRuCutsRn+0fADyS9rBw6hmpNn1bEF1001634s9mAN1P1LPoecPYcvP8/AquBZ6n+ijoNeAFwDXBfedy94/yzS6z3Am/qQ3yvpqoKuAO4rWxvbkOMwEHAd0psdwF/Vo7PeWybxHk0G3thtSI2qjaG28u2fOJ3vy3xlfc7BLi5/Pt+GditTfFl686WqUwiImJGBrkKKyIi5lASSEREzEgSSEREzEgSSEREzEgSSEREzEgSSEREzEjrlrSNrZOkMeDOjkMnAftSDTa7H9gJeBj4a9tf6Xd8EbG5JJBoi5+5mhdrgzIB5DdtH1/2DwG+LOlntq/Z0gUlbWN7fQ9ijQhShRUDxPZtwAeAd051jqTPSDpX0r8CH5J0iKTrJd0h6XJJu0naU9It5fyDJVnSPmX/e5J2kvQ2SXeVBa++0Y/PFzFokkCiLXYsq+vdJunyac67FXj5Fq71UuD1tt8LXAycafsgqiqyP3c1AeEOZd2K11BNufEaSS+imoX3aeDPgGNtHwycMLuPFjGcUoUVbbFZFdYUJpv6e1Oftz0m6ReAXW3/Wzl+EfD58vzfqWbdfS3wl1QzAQv4Znn928BnJF0GfImI2EzuQGLQHEq1suJ0nqpxnW9S3X28iKqh/mCqySe/AWD7HcCfUE0zfpukF8w04IhhlQQSA0PSQcCfAh+rc77tnwCPSXpNOfS7wMTdyDeAtwP32R6nWpr4zVR3Hkja3/YNtv8MWMtz16uICFKFFe33GknfoerGuwZ4d50eWB1OAf63pJ2ougOfCmD7gWq5FCYayL8FLLD9WNn/G0kHUFVrXUM1dXpEdMh07hERMSOpwoqIiBlJFVYMJElnA2/b5PDnbZ8zF/FEbI1ShRURETOSKqyIiJiRJJCIiJiRJJCIiJiRJJCIiJiR/w+/TguQ84y8YgAAAABJRU5ErkJggg==)
%% Cell type:code id: tags:
``` python
def get_latlons(fname, bounds):
#this code is from here: https://github.com/makerportal/GOES-16-Fixed-Grid-Projection/blob/master/goes16_lat_lon_algorithm.py
#Init file obj
g16nc = Dataset(fname)
#Get projection information
proj_info = g16nc.variables['goes_imager_projection']
lon_origin = proj_info.longitude_of_projection_origin
H = proj_info.perspective_point_height+proj_info.semi_major_axis
r_eq = proj_info.semi_major_axis
r_pol = proj_info.semi_minor_axis
# read ints
lat_rad_1d = g16nc.variables['x'][bounds[0]:bounds[1]]
lon_rad_1d = g16nc.variables['y'][bounds[2]:bounds[3]]
# create meshgrid filled with radian angles
lat_rad,lon_rad = np.meshgrid(lat_rad_1d,lon_rad_1d)
# lat/lon calc routine from satellite radian angle vectors
lambda_0 = (lon_origin*np.pi)/180.0
a_var = np.power(np.sin(lat_rad),2.0) + (np.power(np.cos(lat_rad),2.0)*(np.power(np.cos(lon_rad),2.0)+(((r_eq*r_eq)/(r_pol*r_pol))*np.power(np.sin(lon_rad),2.0))))
b_var = -2.0*H*np.cos(lat_rad)*np.cos(lon_rad)
c_var = (H**2.0)-(r_eq**2.0)
r_s = (-1.0*b_var - np.sqrt((b_var**2)-(4.0*a_var*c_var)))/(2.0*