<jupyter_start><jupyter_text>aws s3 sync s3://naip-analytic/az/2013/100cm/index/ shp/naip_index/2013/ --request-payeraws s3 sync s3://naip-analytic/az/2015/100cm/index/ shp/naip_index/2015/ --request-payeraws s3 sync s3://naip-analytic/az/2017/60cm/index/ shp/naip_index/2017/ --request-payer naip13 = gpd.read_file('shp/naip_index/2013/naip_3_13_3_1_az.shp')naip15 = gpd.read_file('shp/naip_index/2015/naip_3_15_1_1_az.shp')naip17 = gpd.read_file('shp/naip_index/2017/naip_3_17_3_1_az.shp')naip13.to_file('shp/az_naip_2013.geojson', driver='GeoJSON')naip15.to_file('shp/az_naip_2015.geojson', driver='GeoJSON')naip17.to_file('shp/az_naip_2017.geojson', driver='GeoJSON')<jupyter_code>import fiona
import numpy as np
import math
import pandas as pd
import geopandas as gpd
import utm
from shapely.geometry import Point, box, Polygon
from shapely.strtree import STRtree
from shapely.geometry import GeometryCollection
g13 = gpd.read_file('shp/az_naip_2013.geojson')
g15 = gpd.read_file('shp/az_naip_2015.geojson')
g17 = gpd.read_file('shp/az_naip_2017.geojson')
coords = (33.439856, -112.066936)
coords = Point(coords[::-1])
bbox = box(*coords.buffer(0.01, cap_style=3).bounds)
bbox
matches = []
%%time
utm.from_latlon(coords.y, coords.x)
%%time
m13 = g13[g13.geometry.intersects(coords)]
#%%timeit pass
gx, gy = str(int(abs(coords.x))), str(int(coords.y))
gx, gy
%%timeit pass
gx, gy = str(int(abs(coords.x))), str(int(coords.y))
match = g13[(g13.DX == gx) & (g13.DY == gy)]
match[match.geometry.intersects(coords)]
def dd_dms(dd):
    is_positive = (dd >= 0)
    dd = abs(dd)
    m, s = divmod(dd*3600, 60)
    d, m = divmod(m, 60)
    d = d if is_positive else -d
    return (abs(int(d)), int(m), round(s,6))
%%timeit pass
yd, ym, ys = dd_dms(coords.y)
xd, xm, xs = dd_dms(coords.x)
match = g13[(g13.DX.astype(np.int) == xd) & (g13.DY.astype(np.int) == yd)]
match = match[match.geometry.intersects(coords)]
%%timeit pass
dx, dy = abs(int(coords.x)), int(coords.y)
match = g13[(g13.DX.astype(np.int) == dx) & (g13.DY.astype(np.int) == dy)]
match = match[match.geometry.intersects(coords)]
match
yd, ym, ys = dd_dms(coords.y)
xd, xm, xs = dd_dms(coords.x)
((str(yd) + str(ym)), (str(xd) + str(xm)))
def dms_str(d, m, s):
    if d < 10 and d > 0:
        d = '0' + str(d)
    if m < 10:
        m = '0' + str(m)
    return (str(d), str(m), str(s))
dms_str(yd, ym, ys)
dms_str(xd, xm, xs)
match.OLAT.values[0], match.OLONG.values[0]
%matplotlib inline

match.plot(column='QUADRANT')
a
g13.plot('DY')
g13.plot('DX')
g13.plot('SY')
g13.plot('OLONG')
g13.plot('OLAT')
%%time
m15 = g15[g15.geometry.intersects(p)]
%%time
m17 = g17[g17.geometry.intersects(p)]
m13.OLAT
m15
m17.OLAT
fn = match.at[midx, 'FileName']
fn = fn.replace('.tif', '.mrf')
%%timeit pass
usgs_id = match.at[midx, 'USGSID']
quad = usgs_id[:-2]
timeres = ['2017/60cm', '2015/100cm', '2013/100cm']
s3s = ['s3://naip-analytic/az/{}/rgbir/{}'.format(res, '/'.join([quad, fn])) for res in timeres]
import rasterio
import os
os.environ['AWS_REQUEST_PAYER'] = 'requester'
os.environ['AWS_DEFAULT_REGION'] = 'us-west-2'
%%time
with rasterio.open(cog_url, driver='MRF') as cog:
    coimg = xr.open_rasterio(cog, chunks={'x': 512, 'y': 512})
naip = 's3://naip-analytic/az/2017/60cm/rgbir/31111/m_3111139_ne_12_h_20170610.mrf'
for k, v in match.items():
    print(k, v.values[0])
%matplotlib inline
gdf17.plot(column='USGSID')
%matplotlib inline
gdf17.plot(column='APFONAME')
gdf17.plot(column='GNIS')
gdf17.plot(column='DY')
gdf17.plot(column='QKEY')
gdf17.columns<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup