Table of Contents

In [ ]:

    

from pysis.isis import campt, camrange
import pvl
from shapely.geometry import box, polygon


    

In [ ]:

    

from pathlib import Path
p = Path('/Volumes/Data/planet4/season2_3_EDRs/')


    

In [ ]:

    

paths = list(p.glob('*.cub'))


    

In [ ]:

    

from planet4.region_data import Inca


    

In [ ]:

    

season2 = [path for path in paths if path.stem[9:] in Inca.season2]


    

In [ ]:

    

season3 = [path for path in paths if path.stem[9:] in Inca.season3]


    

In [ ]:

    

def get_box(path):
    label = pvl.load(str(path) + '.range.lbl.txt')
    minx = label['UniversalGroundRange']['MinimumLongitude']
    maxx = label['UniversalGroundRange']['MaximumLongitude']
    miny = label['LatitudeRange']['MinimumLatitude']
    maxy = label['LatitudeRange']['MaximumLatitude']
    return box(minx, miny, maxx, maxy)


    

In [ ]:

    

[p for p in paths if get_box(p).bounds[0] < 296]


    

In [ ]:

    

boxes = [get_box(p) for p in paths]


    

In [ ]:

    

get_box(paths[0])


    

In [ ]:

    

rightboxes = [b for b in boxes if b.bounds[0] > 295.7]


    

In [ ]:

    

leftboxes = [b for b in boxes if b.bounds[0] < 295.7]


    

In [ ]:

    

len(leftboxes)


    

In [ ]:

    

len(boxes)


    

In [ ]:

    

intersect = rightboxes[0].intersection(rightboxes[1])


    

In [ ]:

    

n = 2
for b in leftboxes[2:]:
    if b.area < 0.1:
        continue
    temp = intersect.intersection(b)
    if temp.area == 0:
        continue
    intersect = intersect.intersection(b)
    n+=1


    

In [ ]:

    

n


    

In [ ]:

    

len(rightboxes)


    

In [ ]:

    

intersect.bounds


    

In [ ]:

    

s2boxes = [get_box(p) for p in season2]
s3boxes = [get_box(p) for p in season3]


    

In [ ]:

    

from shapely.ops import cascaded_union


    

In [ ]:

    

%matplotlib inline


    

In [ ]:

    

import seaborn as sns


    

In [ ]:

    

s2boxes[0].bounds


    

In [ ]:

    

[box.area for box in rightboxes]


    

In [ ]:

    

from matplotlib.patches import Polygon
from itertools import cycle
fig, ax = plt.subplots(figsize=(8, 8))

# for polygon in leftboxes:
#     mpl_poly = Polygon(np.array(polygon.exterior), color='g',
#                        lw=2, alpha=0.1)
#     ax.add_patch(mpl_poly)

for polygon in rightboxes:
    mpl_poly = Polygon(np.array(polygon.exterior), color='g',
                       lw=2, alpha=0.05)
    ax.add_patch(mpl_poly)


ax.add_patch(Polygon(np.array(intersect.exterior), color='r',
                     lw=2))

ax.relim()
ax.autoscale()


    

In [ ]:

    

intersect.bounds


    

In [ ]:

    




    

In [ ]:

    

from shapely import iterops


    

In [ ]:

    

from shapely.examples import intersect


    

In [ ]:

    

intersect??


    

In [ ]:

    

areas = [b.area for b in s2boxes if b.bounds[0] < 296]


    

In [ ]:

    

len(areas)


    

In [ ]:

    

plt.hist(areas)


    

In [ ]:

    

from rtree import index


    

In [ ]:

    

idx = index.Index()


    

In [ ]:

    

for pos, cell in enumerate(boxes):
    idx.insert(pos, cell.bounds)


    

In [ ]:

    

for poly in boxes:
    merged_cells = cascaded_union([boxes[i] for i in idx.intersection(poly.bounds)])
    print(poly.intersection(merged_cells).area)


    

In [ ]: