In [9]:
import geopandas as gpd
import shapely
import pandas as pd
import numpy as np
from math import pi,cos,radians
import pyproj
from shapely.ops import cascaded_union
from shapely.geometry import Point
import shapely.geometry as shpg
import salem
import datetime as dt
import csv
%matplotlib inline
In [10]:
lp_p = 'C:\\Users\\jlandman\\Desktop\\LeBris_Paul\\ElevationChange_GlaciersSupp1km2_LessThan20%DataVoid.shp' # LeBris/Paul shape
foga_p = 'C:\\Users\\jlandman\\Desktop\\DOI-WGMS-FoG-2015-11\\WGMS-FoG-2015-11-A-GENERAL-INFORMATION.csv'
rgi_ak_p = 'C:\\Users\\jlandman\\Desktop\\01_rgi50_Alaska\\01_rgi50_Alaska_selecxt_LeBris_2km.shp'
wgms_subregions_p = 'C:\\Users\\jlandman\\Desktop\\Glacier_Subregions_jhuber\\Glacier_Subregions_jhuber.shp'
In [11]:
lp = gpd.read_file(lp_p)
rgi_ak = gpd.read_file(rgi_ak_p)
lp['ID'] = lp['ID'].astype(int)
foga = pd.read_csv(foga_p, encoding='iso-8859-1')
wgms_subregions = gpd.read_file(wgms_subregions_p)
#wgms_subregions = wgms_subregions[wgms_subregions.WGMS_Code.str.contains('ALA')] #
In [12]:
fogd_templ = pd.DataFrame([], [['POLITICAL_UNIT', 'GLACIER_NAME', 'WGMS_ID', 'YEAR', 'LOWER_BOUND', 'UPPER_BOUND',
'AREA_SURVEY_YEAR', 'AREA_CHANGE', 'AREA_CHANGE_UNCERTAINTY', 'THICKNESS_CHANGE',
'THICKNESS_CHANGE_UNCERTAINTY', 'VOLUME_CHANGE', 'VOLUME_CHANGE_UNCERTAINTY', 'SURVEY_DATE',
'SURVEY_DATE_PLATFORM_METHOD', 'REFERENCE_DATE', 'REFERENCE_DATE_PLATFORM_METHOD',
'INVESTIGATOR', 'SPONSORING_AGENCY', 'REFERENCE', 'REMARKS']])
In [13]:
def coord_unc(x,y):
""" Determine the uncertainity in meters for a given FoG coordinate at its lat/lon (assuming earth was a sphere)
Parameters
-----------
x: longitude in degrees
y: latitude in degrees
Returns
--------
dx, dy: uncertainty in meters
"""
# assert given numbers are floats
# and determine number of decimals given
if isinstance(x, int):
xd = 0
elif isinstance(x, float):
xd = str(x)[::-1].find('.')
else:
raise TypeError('Given x must be float or integer.)')
if isinstance(y, int):
yd = 0
elif isinstance(y, float):
yd = str(y)[::-1].find('.')
else:
raise TypeError('Given y must be float or integer.)')
# maximum error in degrees
xe = 10**(-xd)
ye = 10**(-yd)
# maximum error in meters
dy = ye * ((2.* pi * 6378000.) / 360.)
dx = xe * (2. * pi * cos(radians(y)) * 6378000.) / 360.
print(dx, dy)
return dx, dy
In [14]:
def transform_coords(x,y,in_proj=None,out_proj=None):
assert in_proj is not None,'Input projection may not be None.'
assert out_proj is not None,'Output projection may not be None.'
from_proj = pyproj.Proj(init='epsg:'+str(in_proj))
to_proj = pyproj.Proj(init='epsg:'+str(out_proj))
#x1, y1 = from_proj(x,y)
return pyproj.transform(from_proj, to_proj, x,y)
In [15]:
utm7n = pyproj.Proj(init='epsg:32605') # WGS84 UTM 5N -> LeBris/Paul
latlon = pyproj.Proj(init='epsg:4326') # WGS84 lat/lon -> FoG
In [16]:
lp_extent = lp.total_bounds
In [17]:
# convert UTM5N to WGS84 lat/lon points
x_lp = [lp_extent[0], lp_extent[2]]
y_lp = [lp_extent[1], lp_extent[3]]
x1,y1 = utm7n(x_lp, y_lp)
x_lpll, y_lpll = pyproj.transform(utm7n,latlon,x_lp,y_lp)
In [18]:
foga = foga[(foga.LONGITUDE >= x_lpll[0]) & (foga.LONGITUDE <= x_lpll[1]) & (foga.LATITUDE >= y_lpll[0]) & (foga.LATITUDE <= y_lpll[1])]
In [19]:
# convert WGS84 lat/lon points to UTM5N
xs = foga.LONGITUDE.values
ys = foga.LATITUDE.values
x1,y1 = latlon(xs, ys)
xm, ym = pyproj.transform(latlon,utm7n,xs,ys)
In [ ]:
buffers = [Point(x,y).buffer(coord_unc(transform_coords(x,y,32605,4326)[0],transform_coords(x,y,32605,4326)[1])[1]) for x,y in zip(xm,ym)]
#buffers = [Point(x,y).buffer(np.nanmax(coord_unc(x,y))) for x,y in zip(xs,ys)]
In [21]:
# make the buffer one geometry
buffer_union = cascaded_union(buffers)
In [22]:
# create a series of true/false of those that do not intersect with the buffer should be tried to be directly imported into FoG (there WILL be mistakes) !!!
disjoint = lp.disjoint(buffer_union)
In [23]:
disjoint_ix = disjoint[disjoint == True].index
intersect_ix = disjoint[disjoint == False].index
In [24]:
# get the real geometries
fast_to_FoG = lp[lp.index.isin(disjoint_ix)]
check_closely = lp[lp.index.isin(intersect_ix)]
In [25]:
assert len(fast_to_FoG)+len(check_closely)==len(lp)
In [26]:
len(fast_to_FoG)
Out[26]:
In [50]:
template_a = pd.read_excel('C:\\users\\jlandman\\Desktop\\FoG_Subm_for_pandas.xls', sheetname='A - GENERAL INFORMATION')
template_b = pd.read_excel('C:\\users\\jlandman\\Desktop\\FoG_Subm_for_pandas.xls', sheetname='B - STATE')
template_d = pd.read_excel('C:\\users\\jlandman\\Desktop\\FoG_Subm_for_pandas.xls', sheetname='D - CHANGE')
template_a['GEO-REGION_CODE'] = ''
template_a['GEO-SUBREGION_CODE'] = ''
In [51]:
assign_dict = {(3373, 'AIALIK'): np.nan, # no equivalent
(3543, 'ALLEN'): np.nan, # no equivalent
(92, 'APPLEGATE'): 119,
(102, 'BAKER'): 280,
(105, 'BALTIMORE'): 214,
(165, 'BARNARD'): 184,
(168, 'BARRY'): 971,
(1390, 'BARTLETT'): np.nan, # no equ
(3372, 'BEAR'): np.nan, # no equ
(3376, 'BEAR LAKE'): np.nan, # no eq
(97, 'BELOIT'): 302,
(3377, 'BENCH'): np.nan,
(98,'BLACKSTONE'): 980,
(157,'BRILLIANT'): np.nan,
(162,'BRYN MAWR'): 977,
(169,'CASCADE'): 973,
(100,'CATARACT'): 33,
(180,'CHENEGA'): np.nan,
(3379,'CHERNOF'): np.nan,
(152,'CHILDS'): np.nan,
(176,'CLAREMONT NORTH'): 311,
(177,'CLAREMONT WEST'): 312,
(3544,'COLONY'): 203,
(156,'COLUMBIA (627)'): 976,
(178,'CONTACT'): 325, # is okay as only front variation present in FoG
(167,'COXE'): 978,
(3381,'DESERTED'): np.nan,
(101,'DETACHED'): np.nan,
(3382,'DINGLESTADT'): 364,
(3383,'DOUBLE'): 774, # ATTENTION: Name is different: DOUBLE in FoG, Big RIver Glacier in LP and Big River Lobe Double Glac in GLIMS
(85,'EKLUTNA'): np.nan,
(3790,'EXCELSIOR'): 981,
(86,'EXIT'): np.nan,
(91,'FALLING'): np.nan,
(3386,'FORK TLIKAKILA'): 777,
(3791,'GREWINGK'): np.nan,
(172,'HARRIMAN'): 975,
(160,'HARVARD'): np.nan,
(3390,'HOLGATE'): 982,
(166,'HOLYOKE'): 185,
(3391,'KACHEMAK'): np.nan,
(3310,'KNIK'): 204,
(4333,'KNIK NORTH'): np.nan, # part of KNIK (only one polygon in LP)
(4332,'KNIK SOUTH'): np.nan, # part of KNIK (only one polygon in LP)
(95,'LAWRENCE'): 298,
(173,'LEARNARD'): 100,
(3394,'LITTLE DINGLESTADT'): 360,
(3545,'MARCUS BAKER'): 177,
(96,'MARQUETTE'): 301,
(3546,'MATANUSKA E'): np.nan, # one polygon together with MATANUSKA W in LP
(3547,'MATANUSKA W'): np.nan, # one polygon together with MATANUSKA E in LP
(3396,'MC CARTY'): np.nan,
(158,'MEARES'): np.nan,
(3548,'NELCHINA'): np.nan,
(179,'NELLIE JUAN'): np.nan,
(3414,'NORTH FORK TLIKAKILA'): 769,
(3399,'NORTHEASTERN'): 366, # not absolutely sure if this is the one
(3793,'NORTHWESTERN'): np.nan,
(3794,'NUKA'): np.nan, # VERY DIFFICULT CASE
(103,'PENNIMAN EAST'): np.nan, # not really clear where the transition is
(104,'PENNIMAN WEST'): 279, # probably true but not 100% clear
(174,'PORTAGE'): 292,
(3335,'RED'): np.nan,
(99,'ROARING'): np.nan,
(108,'SADDLEBAG'): np.nan,
(4334,'SCOTT'): np.nan,
(170,'SERPENTINE'): 979,
(3413,'SHAMROCK'): 781,
(151,'SHERIDAN'): np.nan,
(107,'SHERMAN'): np.nan,
(155,'SHOUP'): 970,
(3401,'SKILAK'): np.nan,
(161,'SMITH'): np.nan,
(3549,'SOUTH FORK TSINA'): 179, # very strange naming, but this is it
(1391,'SPENCER'): np.nan, # FoG coordinate was wrong
(171,'SURPRISE'): 974,
(3403,'TANAINA'): 792,
(93,'TAYLOR US'): 310,
(3405,'TAZLINA'): np.nan,
(175,'TEBENKOF'): np.nan,
(3550,'TONSINA'): 195,
(1389,'TRAIL'): 307,
(3551,'TSINA'): 178,
(3407,'TURQUOISE'): 787,
(3408,'TUSTUMENA'): np.nan,
(3409,'TUXEDNI'): 928,
(1387,'UNNAMED US0623'): np.nan, # seems to be part of ID 976 in LP as B table says only 4km long
(106,'UNNAMED US624'): np.nan, # unclear which is meant; probably no equivalent
(154,'VALDEZ'): 211,
(163,'VASSAR'): 215,
(164,'WELLESLEY'): 972,
(94,'WOLVERINE'): np.nan,
(3580,'WOODWORTH'): np.nan,
(153,'WORTHINGTON'): 181,
(3412,'WORTHMANNS'): np.nan,
(159,'YALE'): np.nan,
(3797,'YALIK'):409}
In [52]:
# constants
political_unit = 'US'
geogr_loc = 'Alaska'
year = 2006
survey_date = 20069999
ref_date = 19509999
geo_code = 'ALA'
surv_plat_meth = 'sP'
ref_plat_meth = 'aM'
remarks_a = ''
remarks_b = 'Outlines derived by manual/automated digitizing based on large-scale 7.5 min topographic quadrangle maps - some corrected manually. Geometrical attributes are taken from RGI5.0 and might be from up to four years around the survey date.'
remarks_d = '1950s: Outlines and elevation derived from topographic maps, 2006: elevation derived from SPOT5 HRS (SPIRIT). '
investigator = 'Raymond LE BRIS and Frank PAUL'
spons_agenc = 'Dept. of Geography, University of Zurich, Winterthurerstrasse 190, 8057 Zurich, Switzerland'
ref = 'Le Bris, R. and Paul, F. (2015); Annals of Glaciology, 56(70), pp.184-192.'
In [53]:
lp = lp.to_crs({'init':'epsg:4326'})
In [54]:
lp.to_file('c:\\users\jlandman\\desktop\\LPGeometries.shp', driver='ESRI Shapefile')
In [55]:
rgi_ak.to_file('c:\\users\jlandman\\desktop\\RGIGeometries.shp', driver='ESRI Shapefile')
In [56]:
new_ids = range(5871, 6871, 1)
# create lists with areas in order to check whether LeBris/RGI polygons match
area_list_LP = []
area_list_RGI = []
no_probs_LP = []
no_probs_RGI = []
ct=0
fpct=0
#lp_copy = lp.copy()
for ind, row in lp.iterrows():
ct+=1
#important to reset!
match_ix = None
# this is to check agreement with RGI5.0 on the geometric way....catastrophe (161 matches when maximum intersection plus 10% area tolerance is chosen, 11 matches when RGI point loaction is checked in polygons plus 10% area tolerance)
"""
try:
#match_ix = rgi_ak.intersects(row.geometry).values.tolist().index(True)
match_ix = np.where(rgi_ak.intersects(row.geometry).values == True)[0].tolist()
except ValueError: # in case there is no intersection, there won't be any TRUE in the list
pass
if len(match_ix)==1:
rgi_match = rgi_ak[rgi_ak.index.isin(match_ix)]
elif len(match_ix) > 1: # if there is more than one intersection, take the one with the biggest area overlap
rgi_match = rgi_ak[rgi_ak.index.isin(match_ix)]
rgi_match['union_area'] = lp[lp.index==ind].geometry.intersection(rgi_match.geometry).area
rgi_match = rgi_match[rgi_match.index ==rgi_match['union_area'].argmax()]
else:
rgi_match = pd.DataFrame([])
if rgi_match.empty:
rgi_match = pd.DataFrame(data=np.nan * np.ones(shape=(1,19)), columns=rgi_ak.columns.values)
#print('RGI Area: %s, LP Area: %s' % (rgi_match.Area, row.AREA))
if not (rgi_match.Area.iloc[0]*0.9 < row.AREA < rgi_match.Area.iloc[0]*1.1):
print('Flaechenproblem:', rgi_match.RGIId.iloc[0])
fpct+=1
else:
no_probs_LP.append(row.AREA)
no_probs_RGI.append(rgi_match.Area.iloc[0])
"""
remarks_a_thistime = ''
remarks_d_thistime = ''
if row.ID not in assign_dict.values():
wgms_id = new_ids[0]
new_ids = new_ids[1:]
template_a.loc[ind, 'POLITCAL_UNIT'] = political_unit
template_a.loc[ind, 'WGMS_ID'] = wgms_id
if row.glacier_na != None:
template_a.loc[ind, 'GLACIER_NAME'] = row.glacier_na.upper()
else:
template_a.loc[ind, 'GLACIER_NAME'] = 'UNNAMED %s' % row.ID
remarks_a_thistime = 'Glacier number in GLACIER_NAME from Le Bris, R. and Paul, F. (2015): Annals of Glaciology 56(70), pp.184-192.'
template_a.loc[ind, 'GEOGRAPHICAL_LOCATION_GENERAL'] = geogr_loc
#print(rgi_match.empty)
#if pd.isnull(rgi_match.CenLat.iloc[0]) or pd.isnull(rgi_match.CenLon.iloc[0]):
repres_point_ll = row.geometry.representative_point()
#repres_point_ll = transform_coords(repres_point.x, repres_point.y,32605,4326)
#else:
# repres_point_ll = Point((rgi_match.CenLon.iloc[0], rgi_match.CenLat.iloc[0]))
template_a.loc[ind, 'LATITUDE'] = round(repres_point_ll.y, 6)
template_a.loc[ind, 'LONGITUDE'] = round(repres_point_ll.x, 6)
template_a.loc[ind, 'REMARKS'] = remarks_a + remarks_a_thistime
if row.Type == 1 or row.Type == 2: # means: lake-/tidewater-terminating
template_a.loc[ind, 'FRONTAL_CHARACTERISTICS'] = 4 # means: calving
template_a.loc[ind, 'GEO-REGION_CODE'] = geo_code
for k,r in wgms_subregions.iterrows():
if r.geometry.contains(Point(round(repres_point_ll.x, 6), round(repres_point_ll.y, 6))):
geo_sub_code = r.WGMS_Code
template_a.loc[ind, 'GEO-SUBREGION_CODE'] = geo_sub_code
"""template_b.loc[ind, 'POLITICAL_UNIT'] = political_unit
template_b.loc[ind, 'GLACIER_NAME'] = row.glacier_na.upper()
template_b.loc[ind, 'WGMS_ID'] = wgms_id
template_b.loc[ind, 'YEAR'] = year
template_b.loc[ind, 'MAXIMUM_ELEVATION_OF_GLACIER'] = rgi_match.Zmax.iloc[0]
template_b.loc[ind, 'MEDIAN_ELEVATION_OF_GLACIER'] = rgi_match.Zmed.iloc[0]
template_b.loc[ind, 'MINIMUM_ELEVATION_OF_GLACIER'] = rgi_match.Zmin.iloc[0]
template_b.loc[ind, 'LENGTH'] = rgi_match.Lmax.iloc[0]
template_b.loc[ind, 'AREA'] = rgi_match.Area.iloc[0]
template_b.loc[ind, 'SURVEY_DATE'] = rgi_match.BgnDate.iloc[0]
template_b.loc[ind, 'SURVEY_PLATFORM_METHOD'] = surv_plat_meth
template_b.loc[ind, 'INVESTIGATOR'] = investigator
template_b.loc[ind, 'SPONSORING_AGENCY'] = spons_agenc
template_b.loc[ind, 'REFERENCE'] = ref
template_b.loc[ind, 'REMARKS'] = remarks_b + ''
"""
template_d.loc[ind, 'POLITICAL_UNIT'] = political_unit
if row.glacier_na != None:
template_d.loc[ind, 'GLACIER_NAME'] = row.glacier_na.upper()
else:
template_d.loc[ind, 'GLACIER_NAME'] = 'UNNAMED %s' % row.ID
remarks_d_thistime = 'Glacier number in GLACIER_NAME from Le Bris, R. and Paul, F. (2015): Annals of Glaciology 56(70), pp.184-192. '
template_d.loc[ind, 'WGMS_ID'] = wgms_id
template_d.loc[ind, 'YEAR'] = year
template_d.loc[ind, 'LOWER_BOUND'] = 9999.
template_d.loc[ind, 'UPPER_BOUND'] = 9999.
template_d.loc[ind, 'AREA_SURVEY_YEAR'] = row.AREA
template_d.loc[ind, 'THICKNESS_CHANGE'] = round(row.dh_mean * 1000.) # unit: mm
template_d.loc[ind, 'VOLUME_CHANGE'] = round(row.dh_mean * row.AREA * 1000.) # unit: m*km2*1000 = 1000m3
template_d.loc[ind, 'REFERENCE_DATE'] = ref_date
template_d.loc[ind, 'SURVEY_DATE'] = survey_date
template_d.loc[ind, 'SURVEY_DATE_PLATFORM_METHOD'] = surv_plat_meth
template_d.loc[ind, 'REFERENCE_DATE'] = ref_date
template_d.loc[ind, 'REFERENCE_DATE_PLATFORM_METHOD'] = ref_plat_meth
template_d.loc[ind, 'INVESTIGATOR'] = investigator
template_d.loc[ind, 'SPONSORING_AGENCY'] = spons_agenc
template_d.loc[ind, 'REFERENCE'] = ref
template_d.loc[ind, 'REMARKS'] = remarks_d + remarks_d_thistime + 'Data void on this glacier is %s%%.' % "{0:.2f}".format(row.data_void)
else:
wgms_id = list(assign_dict.keys())[list(assign_dict.values()).index(row.ID)][0]
gname = list(assign_dict.keys())[list(assign_dict.values()).index(row.ID)][1]
"""
template_b.loc[ind, 'POLITICAL_UNIT'] = political_unit
template_b.loc[ind, 'GLACIER_NAME'] = gname.upper()
template_b.loc[ind, 'WGMS_ID'] = wgms_id
template_b.loc[ind, 'YEAR'] = year
template_b.loc[ind, 'MAXIMUM_ELEVATION_OF_GLACIER'] = rgi_match.Zmax.iloc[0]
template_b.loc[ind, 'MEDIAN_ELEVATION_OF_GLACIER'] = rgi_match.Zmed.iloc[0]
template_b.loc[ind, 'MINIMUM_ELEVATION_OF_GLACIER'] = rgi_match.Zmin.iloc[0]
template_b.loc[ind, 'LENGTH'] = rgi_match.Lmax.iloc[0]
template_b.loc[ind, 'AREA'] = rgi_match.Area.iloc[0]
template_b.loc[ind, 'SURVEY_DATE'] = rgi_match.BgnDate.iloc[0]
template_b.loc[ind, 'SURVEY_PLATFORM_METHOD'] = surv_plat_meth
template_b.loc[ind, 'INVESTIGATOR'] = investigator
template_b.loc[ind, 'SPONSORING_AGENCY'] = spons_agenc
template_b.loc[ind, 'REFERENCE'] = ref
template_b.loc[ind, 'REMARKS'] = remarks_b + ''
"""
template_d.loc[ind, 'POLITICAL_UNIT'] = political_unit
template_d.loc[ind, 'GLACIER_NAME'] = gname.upper()
template_d.loc[ind, 'WGMS_ID'] = wgms_id
template_d.loc[ind, 'YEAR'] = year
template_d.loc[ind, 'LOWER_BOUND'] = 9999.
template_d.loc[ind, 'UPPER_BOUND'] = 9999.
template_d.loc[ind, 'AREA_SURVEY_YEAR'] = row.AREA
template_d.loc[ind, 'THICKNESS_CHANGE'] = round(row.dh_mean * 1000.) # unit: mm
template_d.loc[ind, 'VOLUME_CHANGE'] = round(row.dh_mean * row.AREA * 1000.) # unit: m*km2*1000 = 1000m3
template_d.loc[ind, 'REFERENCE_DATE'] = ref_date
template_d.loc[ind, 'SURVEY_DATE'] = survey_date
template_d.loc[ind, 'SURVEY_DATE_PLATFORM_METHOD'] = surv_plat_meth
template_d.loc[ind, 'REFERENCE_DATE'] = ref_date
template_d.loc[ind, 'REFERENCE_DATE_PLATFORM_METHOD'] = ref_plat_meth
template_d.loc[ind, 'INVESTIGATOR'] = investigator
template_d.loc[ind, 'SPONSORING_AGENCY'] = spons_agenc
template_d.loc[ind, 'REFERENCE'] = ref
template_d.loc[ind, 'REMARKS'] = remarks_d + remarks_d_thistime + 'Data void on this glacier is %s%%.' % "{0:.2f}".format(row.data_void)
import matplotlib.pyplot as plt
area_list_LP_na = area_list_LP area_list_RGI_na = area_list_RGI for n in range(len(area_list_LP)): if not (area_list_RGI[n]0.9 < area_list_LP[n] < area_list_RGI[n]1.1): area_list_LP_na[n] = np.nan area_list_RGI_na[n] = np.nan len(area_list_RGI_na)
area_list_RGI_na.count(np.nan)
plt.scatter(area_list_LP_na, area_list_RGI_na) plt.xlim((0,12)) plt.ylim((0,12))
plt.scatter(no_probs_LP, no_probs_RGI) plt.xlim((0,100)) plt.ylim((0,100))
In [57]:
template_a.to_excel('C:\\users\\jlandman\\Desktop\\RL_FP_to_FoG_automatic_A.xls', index=False)
template_b.to_excel('C:\\users\\jlandman\\Desktop\\RL_FP_to_FoG_automatic_B.xls', index=False)
template_d.to_excel('C:\\users\\jlandman\\Desktop\\RL_FP_to_FoG_automatic_D.xls', index=False)
In [44]:
glims = salem.utils.read_shapefile('C:\\Users\\jlandman\\Desktop\\glims_db_20160429\\glims_polygons_alaska.shp')
In [45]:
glims_id_success = {}
glims_id_fail = {}
invest_date = dt.datetime.strptime(str(survey_date), '%Y%m%d') if str(survey_date)[4:]!='9999' else dt.datetime.strptime(str(survey_date)[:4]+'0630', '%Y%m%d')
ct=0
for k, row in template_a.iterrows():
gp = shpg.Point(row.LONGITUDE, row.LATITUDE)
rectangle = shpg.Polygon([(row.LONGITUDE-0.00001, row.LATITUDE-0.00001), (row.LONGITUDE+0.00001, row.LATITUDE-0.00001), (row.LONGITUDE+0.00001, row.LATITUDE+0.00001), (row.LONGITUDE-0.00001, row.LATITUDE+0.00001)])
subset = glims[glims.intersects(rectangle)]
if len(subset) == 1:
if subset.geometry.iloc[0].contains(gp):
glims_id_success[row.WGMS_ID] = subset['glac_id'].iloc[0]
else:
glims_id_fail[row.WGMS_ID] = np.nan
elif len(subset) > 1:
glims_inv = [dt.datetime.strptime(d, '%Y-%m-%dT%H:%M:%S') for d in subset.src_date.values] # GLIMS investigation date
date_diffs = [abs(invest_date-i) for i in glims_inv] # absolute date differences to get the closest investigation date
ix_date_close = date_diffs.index(min(date_diffs)) # find index of the closest date
glims_id_close = subset.iloc[ix_date_close].glac_id # get closest GLIMS ID
glims_id_success[row.WGMS_ID] = glims_id_close
else:
glims_id_fail[row.WGMS_ID] = np.nan
In [46]:
len(glims_id_success)
Out[46]:
In [47]:
len(glims_id_fail)
Out[47]:
In [48]:
len(template_a)
Out[48]:
In [49]:
keyList = glims_id_success.keys()
valueList = glims_id_success.values()
rows = zip(keyList, valueList)
with open('C:\\Users\\jlandman\\Desktop\\Links_GLIMS_FoG_AK.csv', 'w') as f:
writer = csv.writer(f, lineterminator='\n')
for row in rows:
writer.writerow((row))
GLIMS_pairs = list(set(valueList)-set(rgi_ak.GLIMSId.values)) len(GLIMS_pairs)
gl_rgi = [] gl_lp = []
for gl_id in GLIMS_pairs: gl_rgi.append(rgi_ak[rgi_ak.GLIMSId == gl_id]) gl_lp.append(lp[lp.])
lp.ID.values
In [ ]:
In [ ]: