<jupyter_start><jupyter_text>Experience station data and append additional information<jupyter_code>import numpy as np
import pandas as pd
import geopandas as gpd

# import earthpy.clip as cl  # clipping method causes plt.plot() failure
from shapely.geometry import Point, Polygon

import matplotlib.pyplot as plt

import unidecode
import fiona
# configure display options
pd.set_option('display.max_rows', 10)
pd.set_option('display.max_columns', 10)
pd.set_option('display.width', 20)

pd.options.display.float_format = '{:,.5f}'.format

%matplotlib inline
# read excel file
filename = f'../../dataset/bentre/Toa_do_vitridoman_BenTre.xlsx'
excel = pd.ExcelFile(filename)

filename
# read a sheet
df = pd.read_excel(excel, sheet_name=0, dtype=str)

df.info()
df.head(10)
# delete the first column STT
df.drop(columns=['STT'], axis=1, inplace=True)

df.head(10)
# reset column names
df.columns = ['name', 'x', 'y']

df.info()
# drop all rows that have less than 3 columns filled out: missing station name or x or y
df.dropna(thresh=3, inplace=True)

df.head(10)
df.tail(10)
# insert 'code' column
df.insert(0, 'code', df['name'].map(lambda name: unidecode.unidecode(name).upper().strip().replace(' ', '')))
        
df.head()
# convert x and y to floats
df[['x', 'y']] =df[['x', 'y']].astype(float)

df.info()
df.head()
# convert Pandas dataframe to GeoPandas dataframe with correct projection using 'proj4' related codes
df = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df['x'], df['y']), 
                      crs={
                          'proj': 'utm', 
                          'zone': 48, 
                          'ellps': 'WGS84', 
                          'datum': 'WGS84', 
                          'units': 'm', 
                          'no_defs': True})

df.crs
df.head()
# convert station to new projection
df = df.to_crs({'init': 'epsg:4326'})

df.head()
# the boundary of all stations
df.total_bounds
# print station coordinates to manually figure out its height by Google Earch app
[[c, p.y, p.x] for c, p in df[['code', 'geometry']].values]  # station code, longitude, latitude
# manually add 'ft' and 'cm' columns presenting the height from sea level
heights = np.array([7, 13, 22, 17, 31, 19, 9, 2, 12, 24, 23, 0, 8, 0, 1, 0, 0, 6, 5, 10]).astype(float)  # feet

df.insert(4, 'ft', heights)  # feet
df.insert(5, 'cm', heights * 30.48)  # centimeters
        
df.info()
df.head(10)
df.tail(10)
# support output format
fiona.supported_drivers
# save to csv file
df.to_file('../../dataset/csv/station.csv', driver='CSV', index = False)
# save to shapefile
df.to_file('../../dataset/station', encoding='utf-8')<jupyter_output><empty_output><jupyter_text>Visualization<jupyter_code># reload stations from shapefile
df = gpd.read_file('../../dataset/station/station.shp', encoding='utf-8')

df.crs
df.head(10)
df.tail(10)
df.info()
# load Vietnam boundary
vn = gpd.read_file('../../map/VietnamBound4795/VietnamBound47.shp', encoding='utf-8')

vn.crs
vn.head()
# visualize stations on the entire of Vietnam
fig, ax = plt.subplots(figsize=(20, 20))

vn.plot(ax=ax, alpha=0.2, color='grey')
df.plot(ax=ax, markersize=4, color='red', marker='^')
# create bentre dataframe for the boundary of stations
xmin, ymin, xmax, ymax = df.total_bounds

# add margins
margin = 0.3

xmin -= margin
xmax += margin
ymin -= margin
ymax += margin

p1 = Point(xmin, ymin)
p2 = Point(xmin, ymax)
p3 = Point(xmax, ymax)
p4 = Point(xmax, ymin)

np1 = (p1.coords.xy[0][0], p1.coords.xy[1][0])
np2 = (p2.coords.xy[0][0], p2.coords.xy[1][0])
np3 = (p3.coords.xy[0][0], p3.coords.xy[1][0])
np4 = (p4.coords.xy[0][0], p4.coords.xy[1][0])

bentre = gpd.GeoDataFrame(gpd.GeoSeries(Polygon([np1, np2, np3, np4])), columns=['geometry'])

bentre.head()
# join vn and df for bentre
bentre = gpd.overlay(bentre, vn, how='intersection')

bentre.head()
# load Vietnam waterways
waterways = gpd.read_file('../../map/vietnam-latest-free/gis_osm_waterways_free_1/gis_osm_waterways_free_1.shp', encoding='utf-8')

waterways.crs
waterways.info()
waterways.head()
# types of waterways
waterways['fclass'].value_counts()
# visualize stations
fig, ax = plt.subplots(figsize=(20, 20))
ax.set(xlim=[xmin, xmax], ylim=[ymin, ymax])

vn.plot(ax=ax, alpha=0.1, color='grey')

# natural
waterways[waterways['fclass'] == 'river'].plot(ax=ax, alpha=0.4, color='blue', lw=4)
waterways[waterways['fclass'] == 'stream'].plot(ax=ax, alpha=0.4, color='blue', lw=1)

# man made
waterways[waterways['fclass'] == 'canal'].plot(ax=ax, alpha=0.4, color='green', lw=3)
waterways[waterways['fclass'] == 'drain'].plot(ax=ax, alpha=0.4, color='green', lw=1)

df.plot(ax=ax, markersize=200, color='red', marker='^')
# clip waterway into bentre
# does not work
# waterways = cl.clip_shp(waterways, bentre)
# waterways.info()
# waterways.head()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1722567287784679, 'fraction_numerical': 0.03146997929606625, 'mean_word_length': 4.251086956521739, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '78442', 'license': 'Unlicense', 'n_tokens_mistral': 1927, 'n_tokens_neox': 1820, 'n_words': 435, 'path': 'salinization/dev/05. station-pandas/ben-tre-station.ipynb', 'repo_name': 'elbertsoftware/SpringboardAIC', 'size': '4830'}