Get basin from near outlet coordinates¶

Dataset¶

Rivers and Watersheds from River Environment Classification (REC2) New Zealand

Source: https://data-niwa.opendata.arcgis.com/datasets/river-environment-classification-web-map-rec2-v5

Description: This is a Feature Layer Representation of the River Environment Classification (REC2) Version 5, June 2019 - a database of catchment spatial attributes - summarised for every segment in NZ network of rivers. [Rivers as lines and catchments polygons.]

Note: Download both shapefile and csv table.

In [1]:

from shapely.geometry import Polygon
import branca
import fiona
import folium
import geopandas as gpd
import numpy as np
import pandas as pd

def get_shapefile_rows(filename, irows):
    ''' https://gis.stackexchange.com/questions/220023/only-read-specific-rows-of-a-shapefile-with-geopandas-fiona'''
    irows = sorted(irows) # if the elements of the list    are not sorted
    features = []
    imin = min(irows)
    imax = max(irows)+1
    print('getting {} rows from {} to {}'.format(imax-imin+1, imin,imax))
    with fiona.open(filename) as source:
        for i, feature in enumerate(source[imin:imax]):
            irow = i + imin
            if irow in irows:
                features += [feature]
    return gpd.GeoDataFrame.from_features(features)


def get_upstream_ids (tbr, rid, maxiter=50):
    '''
        tbr: pandas table with reaches attributes
        rid: HydroID of outlet reach
    '''
    reaches = [rid]
    i = 0
    newids = reaches
    while i < maxiter:
        newids2 = []
        for rid in newids:
            # get upstream ids
            newids2 += list(tbr[tbr['NextDownID'] == rid]['HydroID'])
        print('iter {}: found {} upstream reaches.'.format(i,len(newids2)))
        reaches += newids2
        if len(newids) == 0:
            i = maxiter + 1
        else:
            i = i + 1
            newids = newids2
        reaches.sort()
    return reaches

from shapely.geometry import Polygon import branca import fiona import folium import geopandas as gpd import numpy as np import pandas as pd def get_shapefile_rows(filename, irows): ''' https://gis.stackexchange.com/questions/220023/only-read-specific-rows-of-a-shapefile-with-geopandas-fiona''' irows = sorted(irows) # if the elements of the list are not sorted features = [] imin = min(irows) imax = max(irows)+1 print('getting {} rows from {} to {}'.format(imax-imin+1, imin,imax)) with fiona.open(filename) as source: for i, feature in enumerate(source[imin:imax]): irow = i + imin if irow in irows: features += [feature] return gpd.GeoDataFrame.from_features(features) def get_upstream_ids (tbr, rid, maxiter=50): ''' tbr: pandas table with reaches attributes rid: HydroID of outlet reach ''' reaches = [rid] i = 0 newids = reaches while i < maxiter: newids2 = [] for rid in newids: # get upstream ids newids2 += list(tbr[tbr['NextDownID'] == rid]['HydroID']) print('iter {}: found {} upstream reaches.'.format(i,len(newids2))) reaches += newids2 if len(newids) == 0: i = maxiter + 1 else: i = i + 1 newids = newids2 reaches.sort() return reaches

In [2]:

crs_wgs84 = 'epsg:4326'
crs_nz = 'epsg:2193'

shpfile_rivers = 'data/River_Lines.shp'
shpfile_wsheds = 'data/Watersheds.shp'

## also need this to get full basin without reading entire shapefile
csvfile_rivers = 'data/River_Lines.csv'
csvfile_wsheds = 'data/Watersheds.csv'


## output files
basin_name = 'Reporoa'

fileout_html = 'basin_{}.html'.format(basin_name)
fileout_rivers_json = '{}_rivers_wgs84.json'.format(basin_name)
fileout_wsheds_json = '{}_watersheds_wgs84.json'.format(basin_name)
fileout_outershed_json = '{}_outershed_wgs84.json'.format(basin_name)

crs_wgs84 = 'epsg:4326' crs_nz = 'epsg:2193' shpfile_rivers = 'data/River_Lines.shp' shpfile_wsheds = 'data/Watersheds.shp' ## also need this to get full basin without reading entire shapefile csvfile_rivers = 'data/River_Lines.csv' csvfile_wsheds = 'data/Watersheds.csv' ## output files basin_name = 'Reporoa' fileout_html = 'basin_{}.html'.format(basin_name) fileout_rivers_json = '{}_rivers_wgs84.json'.format(basin_name) fileout_wsheds_json = '{}_watersheds_wgs84.json'.format(basin_name) fileout_outershed_json = '{}_outershed_wgs84.json'.format(basin_name)

Search for max drainage area within bbox ("outlet")¶

In [3]:

# Reporoa
y0,x0 = -38.436619, 176.339975
bufx=0.1
bufy=0.1
x1,x2 = x0-bufx/2, x0+bufx/2.
y1,y2 = y0-bufy/2, y0+bufy/2.
print(x1,x2,y1,y2)
#x1,x2,y1,y2 = 175,177,-39,-37
bbox = gpd.GeoSeries(Polygon(((x1, y1), (x2,y1), (x2,y2), (x1,y2))))
bbox.crs = crs_wgs84
print(bbox)

bbox_nz = bbox.to_crs(crs_nz)
print(bbox_nz)

# Reporoa y0,x0 = -38.436619, 176.339975 bufx=0.1 bufy=0.1 x1,x2 = x0-bufx/2, x0+bufx/2. y1,y2 = y0-bufy/2, y0+bufy/2. print(x1,x2,y1,y2) #x1,x2,y1,y2 = 175,177,-39,-37 bbox = gpd.GeoSeries(Polygon(((x1, y1), (x2,y1), (x2,y2), (x1,y2)))) bbox.crs = crs_wgs84 print(bbox) bbox_nz = bbox.to_crs(crs_nz) print(bbox_nz)

176.289975 176.38997500000002 -38.486619 -38.386619
0    POLYGON ((176.28998 -38.48662, 176.38998 -38.4...
dtype: geometry
0    POLYGON ((1886967.769 5735060.81, 1895692.538 ...
dtype: geometry

In [4]:

%%time
# 4s
rivers_nz = gpd.read_file(shpfile_rivers, bbox=bbox_nz)
print(len(rivers_nz))

%%time # 4s rivers_nz = gpd.read_file(shpfile_rivers, bbox=bbox_nz) print(len(rivers_nz))

257
CPU times: user 166 ms, sys: 235 ms, total: 401 ms
Wall time: 707 ms

In [5]:

outlet = rivers_nz[rivers_nz['CUM_AREA'] == rivers_nz['CUM_AREA'].max()]
outlet

outlet = rivers_nz[rivers_nz['CUM_AREA'] == rivers_nz['CUM_AREA'].max()] outlet

Out[5]:

	OBJECTID_1	HydroID	NextDownID	CATAREA	CUM_AREA	nzsegment	Enabled	LENGTHDOWN	Headwater	Hydseq	...	headw_dist	segslpmean	LID	reachtype	FROM_NODE	TO_NODE	Shape_Leng	FLOWDIR	GlobalID	geometry
242	119246	119223	118992	295668.106341	4.162465e+09	3119272	1	290959.754971	0	184705	...	163268	0.008472	0	0	124044	123807	1113.01299	1	c85adc9a-c5f3-4c12-bca8-fa7c88bf4b81	LINESTRING (1887546.652 5734893.334, 1887516.5...

1 rows × 30 columns

In [6]:

rivers_nz.crs = crs_nz
rivers = rivers_nz.to_crs(crs_wgs84)
rivers['geoid'] = rivers.index.astype(str)

fig = folium.Figure(width=800, height=400)
m = folium.Map(location=[x0,y0]).add_to(fig)

folium.Marker([y0,x0], tooltip=basin_name, popup=basin_name).add_to(m)

folium.Choropleth(bbox,
                 fill_color='transparent').add_to(m)


layer_name = 'HydroID'

folium.GeoJson(rivers,
               name=layer_name,
               style_function=lambda feature: dict(color='blue', weight=1),
               tooltip=folium.GeoJsonTooltip(fields=[layer_name], aliases=[layer_name], labels=True, sticky=False)
              ).add_to(m)



outlet = rivers[rivers['CUM_AREA'] == rivers['CUM_AREA'].max()]
print(outlet[layer_name])
folium.GeoJson(outlet,
               name=layer_name,
               style_function=lambda feature: dict(color='yellow', weight=2),
               tooltip=folium.GeoJsonTooltip(fields=[layer_name], aliases=[layer_name], labels=True, sticky=False)
              ).add_to(m)


outlet = rivers[rivers['HydroID'] == 118993]
folium.GeoJson(outlet,
               name=layer_name,
               style_function=lambda feature: dict(color='red', weight=2),
               tooltip=folium.GeoJsonTooltip(fields=[layer_name], aliases=[layer_name], labels=True, sticky=False)
              ).add_to(m)

m.fit_bounds(m.get_bounds(),max_zoom=14)
m

rivers_nz.crs = crs_nz rivers = rivers_nz.to_crs(crs_wgs84) rivers['geoid'] = rivers.index.astype(str) fig = folium.Figure(width=800, height=400) m = folium.Map(location=[x0,y0]).add_to(fig) folium.Marker([y0,x0], tooltip=basin_name, popup=basin_name).add_to(m) folium.Choropleth(bbox, fill_color='transparent').add_to(m) layer_name = 'HydroID' folium.GeoJson(rivers, name=layer_name, style_function=lambda feature: dict(color='blue', weight=1), tooltip=folium.GeoJsonTooltip(fields=[layer_name], aliases=[layer_name], labels=True, sticky=False) ).add_to(m) outlet = rivers[rivers['CUM_AREA'] == rivers['CUM_AREA'].max()] print(outlet[layer_name]) folium.GeoJson(outlet, name=layer_name, style_function=lambda feature: dict(color='yellow', weight=2), tooltip=folium.GeoJsonTooltip(fields=[layer_name], aliases=[layer_name], labels=True, sticky=False) ).add_to(m) outlet = rivers[rivers['HydroID'] == 118993] folium.GeoJson(outlet, name=layer_name, style_function=lambda feature: dict(color='red', weight=2), tooltip=folium.GeoJsonTooltip(fields=[layer_name], aliases=[layer_name], labels=True, sticky=False) ).add_to(m) m.fit_bounds(m.get_bounds(),max_zoom=14) m

242    119223
Name: HydroID, dtype: int64

Out[6]:

Get upstream basin of chosen outlet¶

In [7]:

outlet  # 118993

outlet # 118993

Out[7]:

	OBJECTID_1	HydroID	NextDownID	CATAREA	CUM_AREA	nzsegment	Enabled	LENGTHDOWN	Headwater	Hydseq	...	segslpmean	LID	reachtype	FROM_NODE	TO_NODE	Shape_Leng	FLOWDIR	GlobalID	geometry	geoid
237	119016	118993	118992	562488.727169	341314336.0	3119079	1	290959.754971	0	116915	...	0.014783	0	0	123607	123807	880.107246	1	7cc2003a-48d0-44fc-a35d-23b62ad6ff01	LINESTRING (176.29474 -38.4834, 176.29242 -38....	237

1 rows × 31 columns

In [8]:

%%time

# Read all watersheds and rivers but from csv file otherwise too slow
# 6.51s

# needs this to get upstream full basin
tbw = pd.read_csv(csvfile_wsheds)
print(tbw.shape)
tbw.head(2)

tbr = pd.read_csv(csvfile_rivers)
print(tbr.shape)
tbr.head(2)

%%time # Read all watersheds and rivers but from csv file otherwise too slow # 6.51s # needs this to get upstream full basin tbw = pd.read_csv(csvfile_wsheds) print(tbw.shape) tbw.head(2) tbr = pd.read_csv(csvfile_rivers) print(tbr.shape) tbr.head(2)

(593517, 8)
(593517, 31)
CPU times: user 1.92 s, sys: 792 ms, total: 2.71 s
Wall time: 3.02 s

Out[8]:

	OBJECTID_1	OBJECTID	HydroID	NextDownID	CATAREA	CUM_AREA	nzsegment	Enabled	LENGTHDOWN	Headwater	...	headw_dist	segslpmean	LID	reachtype	FROM_NODE	TO_NODE	Shape_Leng	FLOWDIR	Shape__Length	GlobalID
0	1	1	1	9	218553.786394	218553.80	1000005	1	1801.260253	1	...	0	8.388180	0	0	1	2	213.567866	1	213.567866	6304120b-27ba-41db-8dca-1197d7e05ea1
1	2	2	2	9	455995.848576	455995.81	1000003	1	1801.260253	1	...	0	7.678527	0	0	3	2	581.789877	1	581.789877	8d81048e-c368-4dc6-86b1-672ebb2c9f0c

2 rows × 31 columns

In [9]:

%%time
reaches = get_upstream_ids(tbr, rid=int(outlet['HydroID']), maxiter=100)
print(len(reaches))

%%time reaches = get_upstream_ids(tbr, rid=int(outlet['HydroID']), maxiter=100) print(len(reaches))

iter 0: found 2 upstream reaches.
iter 1: found 2 upstream reaches.
iter 2: found 2 upstream reaches.
iter 3: found 2 upstream reaches.
iter 4: found 2 upstream reaches.
iter 5: found 2 upstream reaches.
iter 6: found 4 upstream reaches.
iter 7: found 4 upstream reaches.
iter 8: found 4 upstream reaches.
iter 9: found 6 upstream reaches.
iter 10: found 6 upstream reaches.
iter 11: found 6 upstream reaches.
iter 12: found 8 upstream reaches.
iter 13: found 10 upstream reaches.
iter 14: found 10 upstream reaches.
iter 15: found 6 upstream reaches.
iter 16: found 12 upstream reaches.
iter 17: found 12 upstream reaches.
iter 18: found 14 upstream reaches.
iter 19: found 14 upstream reaches.
iter 20: found 18 upstream reaches.
iter 21: found 16 upstream reaches.
iter 22: found 22 upstream reaches.
iter 23: found 20 upstream reaches.
iter 24: found 18 upstream reaches.
iter 25: found 16 upstream reaches.
iter 26: found 20 upstream reaches.
iter 27: found 18 upstream reaches.
iter 28: found 14 upstream reaches.
iter 29: found 14 upstream reaches.
iter 30: found 12 upstream reaches.
iter 31: found 16 upstream reaches.
iter 32: found 12 upstream reaches.
iter 33: found 16 upstream reaches.
iter 34: found 20 upstream reaches.
iter 35: found 26 upstream reaches.
iter 36: found 20 upstream reaches.

<timed exec>:1: FutureWarning: Calling int on a single element Series is deprecated and will raise a TypeError in the future. Use int(ser.iloc[0]) instead

iter 37: found 18 upstream reaches.
iter 38: found 20 upstream reaches.
iter 39: found 22 upstream reaches.
iter 40: found 24 upstream reaches.
iter 41: found 24 upstream reaches.
iter 42: found 18 upstream reaches.
iter 43: found 17 upstream reaches.
iter 44: found 16 upstream reaches.
iter 45: found 14 upstream reaches.
iter 46: found 14 upstream reaches.
iter 47: found 10 upstream reaches.
iter 48: found 12 upstream reaches.
iter 49: found 8 upstream reaches.
iter 50: found 12 upstream reaches.
iter 51: found 10 upstream reaches.
iter 52: found 6 upstream reaches.
iter 53: found 8 upstream reaches.
iter 54: found 8 upstream reaches.
iter 55: found 7 upstream reaches.
iter 56: found 2 upstream reaches.
iter 57: found 2 upstream reaches.
iter 58: found 2 upstream reaches.
iter 59: found 0 upstream reaches.
iter 60: found 0 upstream reaches.
701
CPU times: user 333 ms, sys: 631 μs, total: 334 ms
Wall time: 331 ms

In [10]:

irows1 = [tbr[tbr['HydroID'] == rid].index[0] for rid in reaches]
irows2 = [tbw[tbw['HydroID'] == rid].index[0] for rid in reaches]

irows1 = [tbr[tbr['HydroID'] == rid].index[0] for rid in reaches] irows2 = [tbw[tbw['HydroID'] == rid].index[0] for rid in reaches]

In [11]:

%%time
dfr = get_shapefile_rows(shpfile_rivers, irows1)
print(dfr.shape)
dfr.head(2)

%%time dfr = get_shapefile_rows(shpfile_rivers, irows1) print(dfr.shape) dfr.head(2)

getting 17283 rows from 101711 to 118993
(701, 30)
CPU times: user 783 ms, sys: 20.4 ms, total: 803 ms
Wall time: 800 ms

Out[11]:

	geometry	OBJECTID_1	HydroID	NextDownID	CATAREA	CUM_AREA	nzsegment	Enabled	LENGTHDOWN	Headwater	...	nzreach_re	headw_dist	segslpmean	LID	reachtype	FROM_NODE	TO_NODE	Shape_Leng	FLOWDIR	GlobalID
0	LINESTRING (1892548.747 5757627.634, 1892593.7...	101733	101712	101998	982613.943736	982613.88	3101348	1	324554.756283	1	...	3030789	0	2.251327	0	0	106499	106500	1167.502712	1	592d7dc7-0928-4d7c-b3e0-389c1ecb1cbe
1	LINESTRING (1892699.926 5756787.134, 1892805.0...	101757	101736	101998	480487.744369	480487.81	3102006	1	324554.756283	1	...	3031057	0	3.943834	0	0	106524	106500	412.618596	1	5803fe51-b3f9-4ef4-9fd8-5154a39b26d7

2 rows × 30 columns

In [12]:

%%time
dfw = get_shapefile_rows(shpfile_wsheds, irows2)
print(dfw.shape)
dfw.head(2)

%%time dfw = get_shapefile_rows(shpfile_wsheds, irows2) print(dfw.shape) dfw.head(2)

getting 17731 rows from 101323 to 119053
(701, 6)
CPU times: user 841 ms, sys: 160 ms, total: 1 s
Wall time: 1.01 s

Out[12]:

	geometry	HydroID	nzsegment	nzreach_re	Area	GlobalID
0	POLYGON ((1892788.395 5758063.294, 1892668.297...	101712	3101348	3030789	982613.943736	51677ddf-9b0d-4098-8907-5869ab52d3f6
1	POLYGON ((1892430.136 5756441.512, 1892430.099...	101736	3102006	3031057	480487.744369	c16a4b88-ce6a-493e-8843-413e86858e05

In [13]:

dfw.crs = crs_nz
dfr.crs = crs_nz

dfw = dfw.to_crs(crs_wgs84)
dfr = dfr.to_crs(crs_wgs84)

# dfw['geoid'] = dfw.index.astype(str)
# dfr['geoid'] = dfr.index.astype(str)

dfw.crs = crs_nz dfr.crs = crs_nz dfw = dfw.to_crs(crs_wgs84) dfr = dfr.to_crs(crs_wgs84) # dfw['geoid'] = dfw.index.astype(str) # dfr['geoid'] = dfr.index.astype(str)

In [14]:

# check if we are reading correctly
r1 = np.sort(reaches)
r2 = np.sort(dfr['HydroID'].values)
r3 = np.sort(dfw['HydroID'].values)
print(np.all(r2 == r1))
print(np.all(r3 == r1))

# check if we are reading correctly r1 = np.sort(reaches) r2 = np.sort(dfr['HydroID'].values) r3 = np.sort(dfw['HydroID'].values) print(np.all(r2 == r1)) print(np.all(r3 == r1))

True
True

In [15]:

# copy some arguments from rivers to watershed for plotting
dfw.set_index('HydroID', inplace=True)
dfr.set_index('HydroID', inplace=True)

dfw.loc[reaches, 'CUM_AREA'] = dfr.loc[reaches, 'CUM_AREA']
dfw.loc[reaches, 'headw_dist'] = dfr.loc[reaches, 'headw_dist']

dfw.reset_index(inplace=True)
dfr.reset_index(inplace=True)

# copy some arguments from rivers to watershed for plotting dfw.set_index('HydroID', inplace=True) dfr.set_index('HydroID', inplace=True) dfw.loc[reaches, 'CUM_AREA'] = dfr.loc[reaches, 'CUM_AREA'] dfw.loc[reaches, 'headw_dist'] = dfr.loc[reaches, 'headw_dist'] dfw.reset_index(inplace=True) dfr.reset_index(inplace=True)

In [16]:

%%time
dfouter = gpd.GeoSeries(dfw.unary_union)
dfouter.crs = crs_wgs84
print(type(dfw), type(dfouter))

%%time dfouter = gpd.GeoSeries(dfw.unary_union) dfouter.crs = crs_wgs84 print(type(dfw), type(dfouter))

<timed exec>:1: DeprecationWarning: The 'unary_union' attribute is deprecated, use the 'union_all()' method instead.

<class 'geopandas.geodataframe.GeoDataFrame'> <class 'geopandas.geoseries.GeoSeries'>
CPU times: user 234 ms, sys: 756 μs, total: 235 ms
Wall time: 231 ms

In [17]:

dfw['CUM_AREA_KM2'] = dfw['CUM_AREA']/1e6
print(dfw['CUM_AREA_KM2'].min(), dfw['CUM_AREA_KM2'].max())

dfw['CUM_AREA_KM2'] = dfw['CUM_AREA']/1e6 print(dfw['CUM_AREA_KM2'].min(), dfw['CUM_AREA_KM2'].max())

0.076633 341.314336

In [18]:

colorscale = branca.colormap.linear.YlOrRd_06.scale(0,400)
print(colorscale(0.5))
colorscale

colorscale = branca.colormap.linear.YlOrRd_06.scale(0,400) print(colorscale(0.5)) colorscale

#ffffb2ff

Out[18]:

In [19]:

layer_name = 'CUM_AREA_KM2'
dfw.__geo_interface__['features'][0]['properties'][layer_name]

layer_name = 'CUM_AREA_KM2' dfw.__geo_interface__['features'][0]['properties'][layer_name]

Out[19]:

0.98261388

Plot folium¶

In [20]:

fig = folium.Figure(width=1200, height=600)
m = folium.Map(location=[x0,y0]).add_to(fig)
folium.Marker([y0,x0], tooltip=basin_name, popup=basin_name).add_to(m)

folium.Choropleth(bbox, fill_color='transparent', name='bbox').add_to(m)

folium.Choropleth(dfouter, fill_color='transparent', name='full basin').add_to(m)

folium.Choropleth(dfr, key_on='feature.id', line_color='blue', name='river').add_to(m)

folium.GeoJson(dfw,
               name=layer_name,
               style_function=lambda feature: {
                   'fillColor': colorscale(feature['properties'][layer_name]),
                   'color': 'transparent',
                   'fillOpacity': 0.4,
               },
              tooltip=folium.GeoJsonTooltip(fields=[layer_name],
                                            aliases=[layer_name],
                                            labels=True,
                                            sticky=False)
              ).add_to(m)

colorscale.caption = layer_name
colorscale.add_to(m)

folium.LayerControl().add_to(m)
m.fit_bounds(m.get_bounds(),max_zoom=12)
print(fileout_html)
m.save(fileout_html)
m

fig = folium.Figure(width=1200, height=600) m = folium.Map(location=[x0,y0]).add_to(fig) folium.Marker([y0,x0], tooltip=basin_name, popup=basin_name).add_to(m) folium.Choropleth(bbox, fill_color='transparent', name='bbox').add_to(m) folium.Choropleth(dfouter, fill_color='transparent', name='full basin').add_to(m) folium.Choropleth(dfr, key_on='feature.id', line_color='blue', name='river').add_to(m) folium.GeoJson(dfw, name=layer_name, style_function=lambda feature: { 'fillColor': colorscale(feature['properties'][layer_name]), 'color': 'transparent', 'fillOpacity': 0.4, }, tooltip=folium.GeoJsonTooltip(fields=[layer_name], aliases=[layer_name], labels=True, sticky=False) ).add_to(m) colorscale.caption = layer_name colorscale.add_to(m) folium.LayerControl().add_to(m) m.fit_bounds(m.get_bounds(),max_zoom=12) print(fileout_html) m.save(fileout_html) m

basin_Reporoa.html

Out[20]:

Save output geojson files¶

In [21]:

print('Saving {}'.format(fileout_rivers_json))
dfr.to_file(fileout_rivers_json, driver="GeoJSON")

print('Saving {}'.format(fileout_wsheds_json))
dfw.to_file(fileout_wsheds_json, driver="GeoJSON")

print('Saving {}'.format(fileout_outershed_json))
dfouter.to_file(fileout_outershed_json, driver="GeoJSON")

print('Saving {}'.format(fileout_rivers_json)) dfr.to_file(fileout_rivers_json, driver="GeoJSON") print('Saving {}'.format(fileout_wsheds_json)) dfw.to_file(fileout_wsheds_json, driver="GeoJSON") print('Saving {}'.format(fileout_outershed_json)) dfouter.to_file(fileout_outershed_json, driver="GeoJSON")

Saving Reporoa_rivers_wgs84.json
Saving Reporoa_watersheds_wgs84.json
Saving Reporoa_outershed_wgs84.json

In [ ]: