{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"import geopandas as gpd\n\npd.options.display.max_colwidth=None\npd.options.display.max_rows=200","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"PREFIX = '../input/cdp-unlocking-climate-solutions/Cities/Cities Disclosing/'\n\ncdp_cities_disclosing_df = pd.concat([pd.read_csv(PREFIX+'2018_Cities_Disclosing_to_CDP.csv'),\\\n                                           pd.read_csv(PREFIX+'2019_Cities_Disclosing_to_CDP.csv'),\\\n                                           pd.read_csv(PREFIX+'2020_Cities_Disclosing_to_CDP.csv')])\ncdp_cities_disclosing_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# # geometries often come in the WKT format, as in case of City Location:\n\n## TODO: That should be fixed to handle missing/Null values\n\n# import geopandas as gpd\n# import shapely.wkt\n\n# # for null geometries\n# from shapely.geometry import GeometryCollection\n\n# geometry = cdp_cities_disclosing_df['City Location'].map(shapely.wkt.loads)\n# cdp_cities_disclosing_df = cdp_cities_disclosing_df.drop('City Location', axis=1)\n# cdp_cities_disclosing_df = gpd.GeoDataFrame(cdp_cities_disclosing_df, crs=\"EPSG:4326\", geometry=geometry)\n# cdp_cities_disclosing_df = cdp_cities_disclosing_df.rename({'geometry': 'City Location'}, axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# ## convert 'City Location' to GeoSeries\n# cdp_cities_disclosing_df['City Location x'] = cdp_cities_disclosing_df['City Location'].str.replace('POINT \\(', '')\n# cdp_cities_disclosing_df['City Location x'] = cdp_cities_disclosing_df['City Location x'].str.replace(')', '')\n# cdp_cities_disclosing_df['City Location LONG'] = cdp_cities_disclosing_df['City Location x'].apply(lambda x: x.split(' ')[0] if (pd.isnull(x)==False) else None)\n# cdp_cities_disclosing_df['City Location LAT'] = cdp_cities_disclosing_df['City Location x'].apply(lambda x: x.split(' ')[1] if (pd.isnull(x)==False) else None)\n# cdp_cities_disclosing_df['City Location'] = gpd.points_from_xy(cdp_cities_disclosing_df['City Location LONG'], cdp_cities_disclosing_df['City Location LAT'])\n# cdp_cities_disclosing_df.loc[cdp_cities_disclosing_df['City Location LONG'].isnull() |  cdp_cities_disclosing_df['City Location LAT'].isnull(), ]\n# cdp_cities_disclosing_df.loc[cdp_cities_disclosing_df['City Location x'].isna(), 'City Location'] = np.NaN\n# cdp_cities_disclosing_df = cdp_cities_disclosing_df.drop(['City Location x', 'City Location LONG', 'City Location LAT'],axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_unique_cities = cdp_cities_disclosing_df[['Account Number', 'City', 'Country', 'City Location']].drop_duplicates().reset_index(drop=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_unique_cities.count()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(f'There are {1032-688} cities without coordinates')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Let's geocode City Name"},{"metadata":{"trusted":true},"cell_type":"code","source":"!pip install geopy","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Geocoding\n\nimport geopy\n\n# create an address if City is not empty (otherwise it would be misleading)\ndf_unique_cities['addr'] = df_unique_cities[['City', 'Country']].fillna('').apply(lambda x: str(x['City']) + ',' + str(x['Country']), axis=1)\ndf_unique_cities.loc[df_unique_cities['City'].isna(), 'addr'] = np.nan\n\nfrom geopy.geocoders import Nominatim\ngeolocator = Nominatim(user_agent=\"Kaggle_CDP_challenge\")\n\nimport geopandas as gpd\n\n# let's geocode\nfrom geopy.geocoders import Nominatim\nfrom geopy.extra.rate_limiter import RateLimiter\ngeocode = RateLimiter(geolocator.geocode, min_delay_seconds=1)\n\ndf_unique_cities['geolocated_string'] = None\n\nfrom tqdm import tqdm\ntqdm.pandas()\n\n# GEOCODE:\ndf_unique_cities.loc[df_unique_cities['addr'].notnull(), 'geoloc_string'] = df_unique_cities.loc[df_unique_cities['addr'].notnull(), 'addr'].apply(geocode)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_unique_cities.iloc[0]['geoloc_string'].raw","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# in GeoCode - LAT, LON are in the reverse order:(location.latitude, location.longitude). As in Google Maps.\n# In GeoPandas (x,y) it suppose to be (LONGITUDE aka x, LATITUDE aka y).\nfrom shapely.geometry import Point\ndf_unique_cities['geoloc_point'] = df_unique_cities['geoloc_string'].apply(lambda loc: Point(tuple(loc.point)[1], tuple(loc.point)[0]) if (pd.isnull(loc)==False) else None)\ndf_unique_cities['geoloc_point'] = gpd.GeoSeries(df_unique_cities['geoloc_point'],  crs=\"EPSG:4326\")\n\n\ndf_unique_cities['geoloc_addr'] = df_unique_cities['geoloc_string'].apply(lambda loc: loc.address if (pd.isnull(loc)==False) else None)\n#df_unique_cities['geoloc_addr'] = df_unique_cities['geoloc_string'].apply(lambda loc: loc.address if loc else None)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_unique_cities[['City', 'Country', 'City Location', 'geoloc_point', 'geoloc_addr']].count()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print('New addresses geolocated:', df_unique_cities.loc[(df_unique_cities['City Location'].isnull()) & (df_unique_cities['geoloc_addr'].notnull())].shape[0])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_unique_cities.dtypes","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# We use shapely.wkt sub-module to parse geoinfo (City Location) in wkt format:\nfrom shapely import wkt\n\n# for WKT to convert correctly, in columns must be valid strings. Execute only once!\ndf_unique_cities.loc[~df_unique_cities['City Location'].isnull(), 'City Location'] = df_unique_cities.loc[~df_unique_cities['City Location'].isnull(), 'City Location'].apply(wkt.loads)\n\n# cast Pandas DF into Geo Pandas DF\ngdf = gpd.GeoDataFrame(df_unique_cities, geometry='City Location', crs=\"EPSG:4326\")\n\npoints_df = gdf.reset_index(drop=True).copy()     #.loc[gdf['Country']=='Poland'].reset_index(drop=True)\npoints_df.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"points_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Checking distances"},{"metadata":{"trusted":true},"cell_type":"code","source":"#points_df['distance'] = \ndf_tmp = points_df.loc[(points_df['City Location'].isna()==False) & (points_df['geoloc_point'].isna()==False)].reset_index(drop=True)\n# I need to change projection. If epsg:4326 then distance will be in degrees. So epsg=3310 should be used according to this:\n# https://gis.stackexchange.com/a/293342\ndf_tmp['City Location'] = df_tmp['City Location'].to_crs(epsg=3310)\ndf_tmp['geoloc_point'] = df_tmp['geoloc_point'].to_crs(epsg=3310)\n# distance in km\ndf_tmp['distance'] = (df_tmp['City Location'].distance(df_tmp['geoloc_point'])/1000).round(2)\n# back to WGS for google maps.\ndf_tmp['City Location'] = df_tmp['City Location'].to_crs(epsg=4326)\ndf_tmp['geoloc_point'] = df_tmp['geoloc_point'].to_crs(epsg=4326)\n\n# over 50 km difference\ndf_tmp[df_tmp['distance']>50]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# over 50 km difference\nprint(df_tmp[df_tmp['distance']>50].shape[0], ' out of ', df_tmp.shape[0], 'were geocoded with over 50 km away. ~', round(df_tmp[df_tmp['distance']>50].shape[0]/df_tmp.shape[0]*100, 2), '%' )\nprint('Mistakes are on both sides. In case of wrong geocoding, I see that comparing countries could help a lot.')\nprint('If country is ok, than often coordinates of a City Location were given in the wrong order or without - sign.')\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Load cities polygons"},{"metadata":{"trusted":true},"cell_type":"code","source":"cities_poly_df = gpd.read_file('../input/1-cities-landscan/1_CITIES_landscan/ne_10m_urban_areas_landscan/ne_10m_urban_areas_landscan.shp')\ncities_poly_df.crs","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import matplotlib.pyplot as plt\n\n# Plot the WGS84\ncities_poly_df.loc[cities_poly_df['name_conve']=='Warsaw'].plot(facecolor='gray');\n# Add title\nplt.title(\"WGS84 (EPSG:4326) projection\");\n# Remove empty white space around the plot\nplt.tight_layout()\n\n# Plot the one with ETRS-LAEA projection (more realistic for this region)\ncities_poly_df.loc[cities_poly_df['name_conve']=='Warsaw'].to_crs(epsg=3035).plot(facecolor='blue');\n# Add title\nplt.title(\"ETRS Lambert Azimuthal Equal Area projection\");\n# Remove empty white space around the plot\nplt.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"poly_df = cities_poly_df.reset_index(drop=True).copy()     #.loc[cities_poly_df['name_conve'].isin(['Warsaw', 'Moscow'])].reset_index(drop=True)\npoly_df.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots()\npoly_df.plot(ax=ax, facecolor='red');\npoints_df.plot(ax=ax, color='blue', markersize=5);\nplt.tight_layout();\n\n\n\nassert points_df.crs == poly_df.crs, 'Coordinate Systems do not match!!'\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Spatial join with citie's polygons"},{"metadata":{"trusted":true},"cell_type":"code","source":"## Assigning polygons:\n\nprint(f'Initially we have {points_df.shape} unique cities')\n# this spatial join works well here, but two problems: some polygons in our dataset are wrong like (Yokohama inside Tokyo). \nresult = gpd.sjoin(points_df, poly_df, how=\"left\", op=\"within\")\nprint(f'after left joining: {result.shape}')\n# We lost the geometry so now, we just have to add it back again.\nresult = result.merge(poly_df[['name_conve', 'geometry']], left_on='name_conve', right_on='name_conve', how='left')\nprint(f'and then: {result.shape}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# But problem: city Tokyo was found in Tokio and in Yokohama polygon \nresult.loc[result['Account Number']==31111].iloc[:,:-1]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# The problem is, that more cities have their polgyons overlapping over each other like Tokyo and Yokohama.\n# I'll assign the polygon with the smaller area in that cases: using 'min_areakm' column (so Yokohama in this case)\n# when duplicated rows (by keys) that should keep only the first, smaller polygon (because sorted ascending)\nresult = result.sort_values(by=['min_areakm']).drop_duplicates(subset=['Account Number', 'City', 'Country', 'City Location'])\n\nprint('Shape: ', result.shape)\n\n# check if Yokohama.\nresult.loc[result['Account Number']==31111].iloc[:,:-3]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## how many points don't have the pologyon assigned\nprint(result.loc[result['City Location'].notnull() & (result['geometry'].isnull())].shape)\nresult.loc[result['City Location'].notnull() & (result['geometry'].isnull())]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# so we have polygons only for \nprint(f'we have polygons for:', result.loc[(result['geometry'].notnull())].shape[0])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# test\nfig, ax = plt.subplots()\nresult.loc[result['name_conve']=='Warsaw', 'geometry'].plot(ax=ax, facecolor='gray');\nresult.loc[result['name_conve']=='Warsaw', 'City Location'].plot(ax=ax, facecolor='red');\nplt.tight_layout();","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# test\nfor city in ['Santiago', 'City of Sydney', 'Kasama', 'San Francisco']:\n    try:\n        fig, ax = plt.subplots()\n        result.loc[result['City']==city, 'geometry'].plot(ax=ax, facecolor='gray');\n        result.loc[result['City']==city, 'City Location'].plot(ax=ax, facecolor='red');\n        plt.tight_layout();\n    except:\n        pass;","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Save results"},{"metadata":{"trusted":true},"cell_type":"code","source":"# create an output directory\n!mkdir /kaggle/working/CDP","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## Saving results\ndf_unique_cities.to_pickle('/kaggle/working/CDP/CDP_cities_with_polygons_v1.pkl')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}