Charles Karpati | Geo-Data Intake and Operations

In this notebook, the basics of working with geographic data are introduced.

• Reading in data (points/ geoms)
• • Convert lat/lng columns to point coordinates
• • Geocoding address to coordinates
• • Changing coordinate reference systems
• • Connecting to PostGisDB's
• Basic Operations
• Saving shape data
• Get Polygon Centroids
• Working with Points and Polygons
• • Map Points and Polygons
• • Get Points in Polygons
• • Create Choropleths
• • Create Heatmaps (KDE?)

Objectives

By the end of this tutorial users should have an understanding of:

• How to read in and process geo-data asa geo-dataframe.
• The Coordinate Reference System and Coordinate Encoding
• Basic geo-visualization strategies

An expansice discursive on programming and cartography can be found here

yourGeoDataframe.dtype

yourGeoDataframe.CRS

As mentioned earlier:

When you use a pandas function to 'read-in' a dataset, the returned value is of a datatype called a 'Dataframe'.

We need a 'Geo-Dataframe', however, to effectively work with spatial data.

While Pandas does not support Geo-Dataframes; Geo-pandas does.

Geopandas has everything you love about pandas, but with added support for geo-spatial data.

Principle benefits of using Geopandas over Pandas when working with spatial data:

• The geopandas plot function will now render a map by default using your 'spatial-geometries' column.
• Libraries exist spatial-operations and interactive map usage.

There are many ways to have our spatial-data be read-in using geo-pandas into a geo-dataframe.

Namely, it means reading in Geo-Spatial-data from a:

1. (.geojson or .shp) file directly using Geo-pandas
2. (.csv, .json) file using Pandas and convert it to Geo-Pandas

• using a prepared 'geometry' column
• by transformting latitude and longitude columns into a 'geometry' column.
• acquiring coordinates from an address
• mapping your non-spatial-data to data-with-space

3. Connecting to a DB

We will review each one below

# BNIA ArcGIS Homepage: https://data-bniajfi.opendata.arcgis.com/
 csa_gdf = intaker.Intake.getData("https://services1.arcgis.com/mVFRs7NF4iFitgbY/ArcGIS/rest/services/Hhchpov/FeatureServer/0/query?where=1%3D1&outFields=*&returnGeometry=true&f=pgeojson")

url = "example.csv"
 geom = 'geometry'
 # An example of loading in an internal BNIA file
 crs = {'init' :'epsg:2248'} 
  
 # Read in the dataframe
 csa_gdf = intaker.Intake.getData(url)

csa_gdf[geom] = csa_gdf[geom].apply(lambda x: loads( str(x) ))

csa_gdf = GeoDataFrame(csa_gdf, crs=crs, geometry=geom)

 
 # If your data has coordinates in two columns run this cell
 # It will create a geometry column from the two.
 # A public dataset is not provided for this example and will not run.
  
 # Load DF HERE. Accidently deleted the link. Need to refind. 
 # Just rely on example 2 for now. 
 """
 exe_df['x'] = pd.to_numeric(exe_df['x'], errors='coerce')
 exe_df['y'] = pd.to_numeric(exe_df['y'], errors='coerce')
 # exe_df = exe_df.replace(np.nan, 0, regex=True)
  
 # An example of loading in an internal BNIA file
 geometry=[Point(xy) for xy in zip(exe_df.x, exe_df.y)]
 exe_gdf = gpd.GeoDataFrame( exe_df.drop(['x', 'y'], axis=1), crs=crs, geometry=geometry)
 """

# Table: Geoloom, 
 # Columns:  
 # In this example, we are going to read in a shapefile
 geoloom_gdf = gpd.read_file("https://services1.arcgis.com/mVFRs7NF4iFitgbY/ArcGIS/rest/services/Geoloom_Crowd/FeatureServer/0/query?where=1%3D1&outFields=*&returnGeometry=true&f=pgeojson");
 # then create columns for its x and y coords
 geoloom_gdf['POINT_X'] = geoloom_gdf['geometry'].centroid.x
 geoloom_gdf['POINT_Y'] = geoloom_gdf['geometry'].centroid.y
 # Now lets just drop the geometry column and save it to have our example dataset. 
 geoloom_gdf = geoloom_gdf.dropna(subset=['geometry'])
 geoloom_gdf.to_csv('example.csv')

# We already know the x and y columns because we just saved them as such.
 geoloom_df['POINT_X'] = pd.to_numeric(geoloom_df['POINT_X'], errors='coerce')
 geoloom_df['POINT_Y'] = pd.to_numeric(geoloom_df['POINT_Y'], errors='coerce')
 # df = df.replace(np.nan, 0, regex=True)
  
 # And filter out for points only in Baltimore City. 
 geoloom_df = geoloom_df[ geoloom_df['POINT_Y'] > 39.3  ]
 geoloom_df = geoloom_df[ geoloom_df['POINT_Y'] < 39.5  ]

crs = {'init' :'epsg:2248'} 
 geometry=[Point(xy) for xy in zip(geoloom_df['POINT_X'], geoloom_df['POINT_Y'])]
 geoloom_gdf = gpd.GeoDataFrame( geoloom_df.drop(['POINT_X', 'POINT_Y'], axis=1), crs=crs, geometry=geometry)
 # 39.2904° N, 76.6122°

ax = csa_gdf.plot(column='hhchpov18', edgecolor='black')
  
 # now plot our points over it.
 geoloom_gdf.plot(ax=ax, color='red')
  
 plt.show()

# Change these values in the cell below using different geographic reference codes will change those parameters
 tract = '*'
 county = '510' # '059' # 153 '510'
 state = '24' #51
  
 # Specify the download parameters the function will receieve here
 tableId = 'B19049' # 'B19001'
 year = '17'
 saveAcs = True 

IPython.core.display.HTML("")
 # state, county, tract, tableId, year, saveOriginal, save 
 df = retrieve_acs_data(state, county, tract, tableId, year)
 df.head(1)
 df.to_csv('tracts_data.csv')

crosswalk.tail(1)

 merged_df_geom = mergeDatasets(left_ds=df, right_ds=crosswalk, crosswalk_ds=False,
                   left_col='tract', right_col='TRACTCE10',
                   crosswalk_left_col = False, crosswalk_right_col = False,
                   merge_how='outer', # left right or columnname to retrieve
                   interactive=False)
 merged_df_geom.head(1)

Hhchpov = gpd.read_file("https://services1.arcgis.com/mVFRs7NF4iFitgbY/ArcGIS/rest/services/Hhchpov/FeatureServer/0/query?where=1%3D1&outFields=*&returnGeometry=true&f=pgeojson")
 Hhchpov = Hhchpov[['CSA2010', 'hhchpov15',	'hhchpov16',	'hhchpov17',	'hhchpov18', 'geometry']]
 Hhchpov.to_file("Hhchpov.geojson", driver='GeoJSON')
 Hhchpov.to_csv('Hhchpov.csv')
 gpd.read_file("Hhchpov.geojson").head(1)

# df.merge(crosswalk, left_on='tract', right_on='TRACTCE10')

merged_df = mergeDatasets(left_ds=merged_df_geom, right_ds=Hhchpov, crosswalk_ds=False,
                   left_col='CSA2010', right_col='CSA2010',
                   crosswalk_left_col = False, crosswalk_right_col = False,
                   merge_how='outer', # left right or columnname to retrieve
                   interactive=False)

# The attributes are what we will use.
 in_crs = 2248 # The CRS we recieve our data 
 out_crs = 4326 # The CRS we would like to have our data represented as
 geom = 'geometry' # The column where our spatial information lives.
 
 # To create this dataset I had to commit a full outer join. 
 # In this way geometries will be included even if there merge does not have a direct match. 
 # What this will do is that it means at least one (near) empty record for each community will exist that includes (at minimum) the geographic information and name of a Community.
 # That way if no point level information existed in the community, that during the merge the geoboundaries are still carried over.
 
 # Primary Table
 # Description: I created a public dataset from a google xlsx sheet 'Bank Addresses and Census Tract'.
 # Table: FDIC Baltimore Banks
 # Columns: Bank Name, Address(es), Census Tract
 left_ds = 'tracts_data.csv'
 left_col = 'tract'
 
 # Crosswalk Table
 # Table: Crosswalk Census Communities
 # 'TRACT2010', 'GEOID2010', 'CSA2010'
 crosswalk_ds = 'CSA-to-Tract-2010.csv'
 use_crosswalk = True
 crosswalk_left_col = 'TRACTCE10'
 crosswalk_right_col = 'CSA2010'
 
 # Secondary Table
 # Table: Baltimore Boundaries => HHCHPOV
 # 'TRACTCE10', 'GEOID10', 'CSA', 'NAME10', 'Tract', 'geometry'
 right_ds = 'Hhchpov.geojson'
 right_col ='CSA2010'
 
 interactive = True
 merge_how = 'outer'
 
 # reutrns a pandas dataframe
 mergedf = merge.mergeDatasets( left_ds=left_ds, left_col=left_col, 
               crosswalk_ds=crosswalk_ds,
               crosswalk_left_col = crosswalk_left_col, crosswalk_right_col = crosswalk_right_col,
               right_ds=right_ds, right_col=right_col, 
               merge_how=merge_how, interactive = interactive )

# mergedf[geom] = mergedf[geom].apply(lambda x: loads( str(x) ) ) 
 
 # Process the dataframe as a geodataframe with a known CRS and geom column 
 mergedGdf = GeoDataFrame(mergedf, crs=in_crs, geometry=geom) 

geoloom = geoloom.dropna(subset=['geometry'])
 geoloom = geoloom.drop(columns=['geometry','GlobalID', 'POINT_X',	'POINT_Y'])
 geoloom.head(1)

    'Location' : pd.Series([
     '100 N. Holliday St, Baltimore, MD 21202',
     '200 E Pratt St, Baltimore, MD',
     '2401 Liberty Heights Ave, Baltimore, MD',
     '201 E Pratt St, Baltimore, MD',
     '3501 Boston St, Baltimore, MD',
     '857 E Fort Ave, Baltimore, MD',
     '2413 Frederick Ave, Baltimore, MD'
   ]),
     'Address' : pd.Series([ 
     'Baltimore City Council',
     'The Gallery at Harborplace',
     'Mondawmin Mall',
     'Harborplace',
     'The Shops at Canton Crossing',
     'Southside Marketplace',
     'Westside Shopping Center'
   ])
 })
 
 address_df.head()

addrCol = 'Location'

 # In this example we retrieve and map a dataset with no lat/lng but containing an address
 
 # In this example our data is stored in the 'STREET' attribute
 geometry = []
 geolocator = Nominatim(user_agent="my-application")
 
 for index, row in address_df.iterrows():
   # We will try and return an address for each Street Name
   try: 
       # retrieve the geocoded information of our street address
       geol = geolocator.geocode(row[addrCol], timeout=None)
 
       # create a mappable coordinate point from the response object's lat/lang values.
       pnt = Point(geol.longitude, geol.latitude)
       
       # Append this value to the list of geometries
       geometry.append(pnt)
       
   except: 
       # If no street name was found decide what to do here.
       # df.loc[index]['geom'] = Point(0,0) # Alternate method
       geometry.append(Point(0,0))
       
 # Finally, we stuff the geometry data we created back into the dataframe
 address_df['geometry'] = geometry

gdf = gdf[ gdf.centroid.y > 39.3  ]
 gdf = gdf[ gdf.centroid.y < 39.5  ]

ax = csa_gdf.plot(column='hhchpov18', edgecolor='black')
 
 # now plot our points over it.
 geoloom_gdf.plot(ax=ax, color='red')

'''
 conn = psycopg2.connect(host='', dbname='', user='', password='', port='')
 
 # DB Import Method One
 sql1 = 'SELECT the_geom, gid, geogcode, ooi, address, addrtyp, city, block, lot, desclu, existing FROM housing.mdprop_2017v2 limit 100;'
 pointData = gpd.io.sql.read_postgis(sql1, conn, geom_col='the_geom', crs=2248)
 pointData = pointData.to_crs(epsg=4326)
 
 # DB Import Method Two
 sql2 = 'SELECT ST_Transform(the_geom,4326) as the_geom, ooi, desclu, address FROM housing.mdprop_2017v2;'
 pointData = gpd.GeoDataFrame.from_postgis(sql2, conn, geom_col='the_geom', crs=4326)
 pointData.head()
 pointData.plot()
 '''

Create Geospatial Functions

Operations:

• Reading in data (points/ geoms) -- Convert lat/lng columns to point coordinates -- Geocoding address to coordinates -- Changing coordinate reference systems -- Connecting to PostGisDB's
• Basic Operations
• Saving shape data
• Get Polygon Centroids
• Working with Points and Polygons -- Map Points and Polygons -- Get Points in Polygons

Input(s):

• Dataset (points/ bounds) url
• Points/ bounds geometry column(s)
• Points/ bounds crs's
• Points/ bounds mapping color(s)
• New filename

Output: File

This function will handle common geo spatial exploratory methods. It covers everything discussed in the basic operations and more!

#
 # Work With Geometry Data
 # Description: geomSummary, getPointsInPolygons, getPolygonOnPoints, mapPointsInPolygons, getCentroids
 #
 def workWithGeometryData(method=False, df=False, polys=False, ptsCoordCol=False, polygonsCoordCol=False, polyColorCol=False, polygonsLabel='polyOnPoint', pntsClr='red', polysClr='white', interactive=False):
   def geomSummary(df): return type(df), df.crs, df.columns;
   def getCentroid(df, col): return df[col].representative_point() # df['geometry'].centroid
 
   # To 'import' a script you wrote, map its filepath into the sys
   def getPolygonOnPoints(pts, polygons, ptsCoordCol, polygonsCoordCol, polygonsLabel, interactive):
       count = 0
       # We're going to keep a list of how many points we find.
       boundaries = []
 
       # Loop over polygons with index i.
       for i, pt in pts.iterrows():
           # print('Searching for point within Geom:', pt )
           # Only one Label is accepted.
           poly_on_this_point = []
           # Now loop over all polygons with index j.
           for j, poly in polygons.iterrows():
               if poly[polygonsCoordCol].contains(pt[ptsCoordCol]):
                   # Then it's a hit! Add it to the list
                   poly_on_this_point.append(poly[polygonsLabel])
                   count = count + 1
                   # pts = pts.drop([j])
 
           # We could do all sorts, like grab a property of the
           # points, but let's just append the number of them.
           boundaries.append(poly_on_this_point)
           clear_output(wait=True)
 
       # Add the number of points for each poly to the dataframe.
       pts = pts.assign(CSA2010 = boundaries)
       if (interactive):
         print( 'Total Points: ', (pts.size / len(pts.columns) ) )
         print( 'Total Points in Polygons: ', count )
         print( 'Prcnt Points in Polygons: ', count / (pts.size / len(pts.columns) ) )
       return pts
 
   # To 'import' a script you wrote, map its filepath into the sys
   def getPointsInPolygons(pts, polygons, ptsCoordCol, polygonsCoordCol, interactive):
     count = 0
     total = pts.size / len(pts.columns)
     # We're going to keep a list of how many points we find.
     pts_in_polys = []
 
     # Loop over polygons with index i.
     for i, poly in polygons.iterrows():
         # print('Searching for point within Geom:', poly )
         # Keep a list of points in this poly
         pts_in_this_poly = 0
 
         # Now loop over all points with index j.
         for j, pt in pts.iterrows():
             if poly[polygonsCoordCol].contains(pt[ptsCoordCol]):
                 # Then it's a hit! Add it to the list,
                 pts_in_this_poly += 1
                 # and drop it so we have less hunting. # pts = pts.drop([j])
 
         # We could do all sorts, like grab a property of the
         # points, but let's just append the number of them.
         pts_in_polys.append(pts_in_this_poly)
         if (interactive): print('Found this many points within the Geom:', pts_in_this_poly )
         count += pts_in_this_poly
         clear_output(wait=True)
 
     # Add the number of points for each poly to the dataframe.
     polygons['pointsinpolygon'] = gpd.GeoSeries(pts_in_polys)
     if (interactive):
       print( 'Total Points: ', total )
       print( 'Total Points in Polygons: ', count )
       print( 'Prcnt Points in Polygons: ', count / total )
     return polygons
 
   def mapPointsandPolygons(pnts, polys, pntsCl, polysClr, polyColorCol):
     print('mapPointsandPolygons');
     # We restrict to South America.
     ax = 1
     if polyColorCol:
       ax = polys.plot( column=polyColorCol, legend=True)
     else:
       ax = polys.plot( color=polysClr, edgecolor='black')
 
     # We can now plot our ``GeoDataFrame``.
     pnts.plot(ax=ax, color=pntsClr)
 
     return plt.show()
 
   if method=='summary': return geomSummary(df);
   if method=='ponp': return getPolygonOnPoints(df, polys, ptsCoordCol, polygonsCoordCol, polygonsLabel, interactive);
   if method=='pinp': return getPointsInPolygons(df, polys, ptsCoordCol, polygonsCoordCol, interactive);
   if method=='pandp': return mapPointsandPolygons(df, polys, pntsClr, polysClr, polyColorCol);
   if method=='centroid': return getCentroid(df, col);

                plot_points=False, pt_radius=15, \
                 draw_heatmap=True, heat_map_weights_col=None, \
                 heat_map_weights_normalize=True, heat_map_radius=15):
     """Creates a map given a dataframe of points. Can also produce a heatmap overlay
     Arg:
         df: dataframe containing points to maps
         lat_col: Column containing latitude (string) 
     """
 
     ## center map in the middle of points center in
     middle_lat = df[lat_col].median()
     middle_lon = df[lon_col].median()
 
     curr_map = folium.Map(location=[middle_lat, middle_lon],
                           zoom_start=zoom_start)
 
     # add points to map
     if plot_points:
         for _, row in df.iterrows():
             folium.CircleMarker([row[lat_col], row[lon_col]],
                                 radius=pt_radius,
                                 popup=row['name'],
                                 fill_color="#3db7e4", # divvy color
                                ).add_to(curr_map)
 
     # add heatmap
     if draw_heatmap:
         # convert to (n, 2) or (n, 3) matrix format
         if heat_map_weights_col is None:
             stations = zip(df[lat_col], df[lon_col])
         else:
             # if we have to normalize
             if heat_map_weights_normalize:
                 df[heat_map_weights_col] = \
                     df[heat_map_weights_col] / df[heat_map_weights_col].sum()
 
             stations = zip(df[lat_col], df[lon_col], df[heat_map_weights_col])
 
         curr_map.add_child(plugins.HeatMap(stations, radius=heat_map_radius))
 
     return curr_map

 # reverseGeoCode, readFile, getGeoParams, main
 def readInGeometryData(url=False, porg=False, geom=False, lat=False, lng=False, revgeocode=False,  save=False, in_crs=4326, out_crs=False):
 
   def reverseGeoCode(df, lat ):
     # STREET	CITY	STATE ZIP NAME
     # , format_string="%s, BALTIMORE MD"
     geometry = []
     geolocator = Nominatim(user_agent="my-application")
     for index, row in df.iterrows():
       try:
           geol = geolocator.geocode(row[lat], timeout=None)
           pnt = Point(geol.longitude, geol.latitude)
           geometry.append(pnt)
       except:
           geometry.append(Point(-76, 39) )
           print(row[lat]);
     return geometry
 
   def readFile(url, geom, lat, lng, revgeocode, in_crs, out_crs):
     df = False
     gdf = False
     ext = isinstance(url, pd.DataFrame)
     if ext: ext='csv'
     else: ext = url[-3:]
 
     #XLS
     # b16 = pd.read_excel('Jones.BirthsbyCensus2016.XLS', sheetname='Births')
 
     # The file extension is used to determine the appropriate import method.
     if ext in ['son', 'kml', 'shp', 'pgeojson']: gdf = gpd.read_file(url)
     if ext == 'csv':
       df = url if isinstance(url, pd.DataFrame) else pd.read_csv(url)
       # Read using Geom, Lat, Lat/Lng, revGeoCode
       if revgeocode=='y': df['geometry'] = reverseGeoCode(df, lat)
       elif geom: df['geometry'] = df[geom].apply(lambda x: loads( str(x) ))
       elif lat==lng: df['geometry'] = df[lat].apply(lambda x: loads( str(x) ))
       elif lat!=lng: df['geometry'] = gpd.points_from_xy(df[lng], df[lat]);
 
       gdf = GeoDataFrame(df, crs=in_crs, geometry='geometry') #crs=2248
       if not out_crs == in_crs: gdf = gdf.to_crs(epsg=out_crs)
     return gdf
 
   def getGeoParams(url, porg, geom, lat, lng, revgeocode, save, in_crs, out_crs):
     addr=False
 
     if not url: url = input("Please enter the location of your dataset: " )
     # if url[-3:] == 'csv' :
     #   df = pd.read_csv(url,index_col=0,nrows=1)
     #   print(df.columns)
 
     # Geometries inside
     if geom and not (lat and lng): porg = 'g'
     # Point data inside
     elif not geom and lat or lng:
       porg = 'p';
       if not lat: lat = lng
       if not lng: lng = lat
 
     # If the P/G could not be infered...
     if not (porg in ['p', 'g']):
       if not revgeocode in ['y', 'n']: revgeocode = input("Do your records need reverse geocoding: (Enter: y/n') " )
       if revgeocode == 'y': porg = 'p'; lng = lat = input("Please enter the column name where the address is stored: " );
       elif revgeocode == 'n': porg = input("""Do the records in this dataset use (P)oints or (g)eometric polygons?: (Enter: 'p' or 'g') """ );
       else: return getGeoParams(url, porg, geom, lat, lng, revgeocode, save, in_crs, out_crs);
 
       if porg=='p':
           if not lat: lat = input("Please enter the column name where the latitude coordinate is stored: " );
           if not lng: lng = input("Please enter the column name where the longitude cooridnate is stored: (Could be same as the lat) " );
       elif porg=='g':
         if not geom: geom = input("Please enter column name where the geometry data is stored: (*optional, skip if unkown)" );
       else: return getGeoParams(url, porg, geom, lat, lng, revgeocode, save, in_crs, out_crs)
 
     if not out_crs: out_crs=in_crs
 
     return url, porg, geom, lat, lng, revgeocode, save, in_crs, out_crs
 
   # This function uses all the other functions
   def main(url, porg, geom, lat, lng, revgeocode, save, in_crs, out_crs):
 
     # Check for missing values. retrieve them
     if (isinstance(url, pd.DataFrame)): print('Converting DF to GDF')
     elif (not (url and porg) ) or (
         not (porg == 'p' or porg == 'g') ) or (
         porg == 'g' and not geom) or (
         porg == 'p' and (not (lat and lng) ) ):
       return readInGeometryData( *getGeoParams(url, porg, geom, lat, lng, revgeocode, save, in_crs, out_crs) );
 
     # print(f"RECIEVED url: {url}, \r\n porg: {porg}, \r\n geom: {geom}, \r\n lat: {lat}, \r\n lng: {lng}, \r\n revgeocode: {revgeocode}, \r\n in_crs: {in_crs}, \r\n out_crs: {out_crs}")
 
     # Quit if the Columns dont exist -> CSV Only
     # status = checkColumns(url, geom, lat, lng)
     # if status == False: print('A specified column does not exist'); return False;
 
     # Perform operation
     gdf = readFile(url, geom, lat, lng, revgeocode, in_crs, out_crs)
 
     # Tidy up
 
     # Save
     # if save: saveGeoData(gdf, url, fileName, driver='esri')
 
     return gdf
 
   return main(url, porg, geom, lat, lng, revgeocode, save, in_crs, out_crs)

# draw_heatmap, cluster_points, plot_points,
 def map_points(data, lat_col='POINT_Y', lon_col='POINT_X', zoom_start=11, plot_points=True, cluster_points=False,
                pt_radius=15, draw_heatmap=False, heat_map_weights_col=None, heat_map_weights_normalize=True,
                heat_map_radius=15, popup=False):
     """Creates a map given a dataframe of points. Can also produce a heatmap overlay
 
     Arg:
         df: dataframe containing points to maps
         lat_col: Column containing latitude (string)
         lon_col: Column containing longitude (string)
         zoom_start: Integer representing the initial zoom of the map
         plot_points: Add points to map (boolean)
         pt_radius: Size of each point
         draw_heatmap: Add heatmap to map (boolean)
         heat_map_weights_col: Column containing heatmap weights
         heat_map_weights_normalize: Normalize heatmap weights (boolean)
         heat_map_radius: Size of heatmap point
 
     Returns:
         folium map object
     """
     df = data.copy()
     ## center map in the middle of points center in
     middle_lat = df[lat_col].median()
     middle_lon = df[lon_col].median()
 
     curr_map = folium.Map(location=[middle_lat, middle_lon], zoom_start=zoom_start)
 
     # add points to map
     if plot_points:
       for _, row in df.iterrows():
         # print([row[lat_col], row[lon_col]], row[popup])
         folium.CircleMarker([row[lat_col], row[lon_col]],
           radius=pt_radius,
           popup=row[popup],
           fill_color="#3db7e4", # divvy color
         ).add_to(curr_map)
     if cluster_points:
       marker_cluster = MarkerCluster().add_to(curr_map)
       for index, row in df.iterrows():
         folium.Marker( location=[row[lat_col],row[lon_col]], popup=row[popup], icon=None ).add_to(marker_cluster)
 
     # add heatmap
     if draw_heatmap:
         # convert to (n, 2) or (n, 3) matrix format
         if heat_map_weights_col is None:
             stations = zip(df[lat_col], df[lon_col])
         else:
             # if we have to normalize
             if heat_map_weights_normalize:
                 df[heat_map_weights_col] = \
                     df[heat_map_weights_col] / df[heat_map_weights_col].sum()
 
             stations = zip(df[lat_col], df[lon_col], df[heat_map_weights_col])
 
         curr_map.add_child(plugins.HeatMap(stations, radius=heat_map_radius))
 
     return curr_map

Processing Geometry is tedius enough to merit its own handler

               pt_radius=15, draw_heatmap=False, heat_map_weights_col=None, heat_map_weights_normalize=True,
                heat_map_radius=15, popup='ProjNm')

As you can see we have a lot of points. Lets see if there is any better way to visualize this.

Example: Using the advanced Functions

Playing with Points: Geoloom

Points In Polygons

The red dots from when we mapped the geoloom points above were a bit too noisy.

Lets create a choropleth instead!

We can do this by aggregating by CSA.

To do this, start of by finding which points are inside of which polygons!

Since the geoloom data does not have a CSA dataset, we will need merge it to one that does!

Lets use the childhood poverty link from example one and load it up because it contains the geometry data and the csa labels.

# csa_gdf = Intake.getData('https://raw.githubusercontent.com/bniajfi/bniajfi/main/CSA-to-Tract-2010.csv')

# BNIA ArcGIS Homepage: https://data-bniajfi.opendata.arcgis.com/
 csa_gdf_url = "https://services1.arcgis.com/mVFRs7NF4iFitgbY/ArcGIS/rest/services/Hhchpov/FeatureServer/0/query?where=1%3D1&outFields=*&returnGeometry=true&f=pgeojson"
 csa_gdf = readInGeometryData(url=csa_gdf_url, porg=False, geom='geometry', lat=False, lng=False, revgeocode=False,  save=False, in_crs=2248, out_crs=False)

And now lets pull in our geoloom data. But to be sure, drop the empty geometry columns or the function directly below will now work.

geoloom_gdf = readInGeometryData(url=geoloom_gdf_url, porg=False, geom='geometry', lat=False, lng=False, revgeocode=False,  save=False, in_crs=4326, out_crs=False)
 geoloom_gdf = geoloom_gdf.dropna(subset=['geometry'])
 # geoloom_gdf = geoloom_gdf.drop(columns=['POINT_X','POINT_Y'])
 geoloom_gdf.head(1)

And now use a point in polygon method 'ponp' to get the CSA2010 column from our CSA dataset added as a column to each geoloom record.

You'll see you have a 'pointsinpolygons' column now.

Polygons in Points

Alternately, you can run the ponp function and have returned the geoloom dataset

We can count the totals per CSA using value_counts

Alternately, we could map the centroid of boundaries within another boundary to find boundaries within boundaries

geoloom_w_csas['POINT_X'] = geoloom_w_csas.centroid.x
 
 # We already know the x and y columns because we just saved them as such.
 geoloom_w_csas['POINT_X'] = pd.to_numeric(geoloom_w_csas['POINT_X'], errors='coerce')
 geoloom_w_csas['POINT_Y'] = pd.to_numeric(geoloom_w_csas['POINT_Y'], errors='coerce')
 # df = df.replace(np.nan, 0, regex=True)
 
 # And filter out for points only in Baltimore City. 
 geoloom_w_csas = geoloom_w_csas[ geoloom_w_csas['POINT_Y'] > 39.3  ]
 geoloom_w_csas = geoloom_w_csas[ geoloom_w_csas['POINT_Y'] < 39.5  ]

Folium

But if that doesn't do it for you, we can also create heat maps and marker clusters

               pt_radius=1, draw_heatmap=True, heat_map_weights_col=None, heat_map_weights_normalize=True,
                heat_map_radius=15, popup='ProjNm')

And Time Sliders

Choropleth Timeslider

https://github.com/python-visualization/folium/blob/master/examples/TimeSliderChoropleth.ipynb

import numpy as np
 import pandas as pd
 from branca.colormap import linear
 # conditionally loaded ->  from dataplay import geoms
 
 u = intaker.Intake
 rdf = u.getData('https://services1.arcgis.com/mVFRs7NF4iFitgbY/ArcGIS/rest/services/Biz1_/FeatureServer/0/query?where=1%3D1&outFields=*&returnGeometry=true&f=pgeojson')
 # rdf.set_index('CSA2010', drop=True, inplace=True)
 rdf.drop(labels=['OBJECTID_1', 'Shape__Area', 'Shape__Length'], axis=1, inplace=True)
 
 ndf = rdf.filter(regex='biz1|CSA2010', axis=1)
 
 # Calculate number of years available
 n_periods = len(ndf.columns) - 1
 # Get starting year.
 startAt = "20"+ndf.columns[1][-2:]
 
 # Create a 'YEAR' index with the assumption that all following years exist
 datetime_index = pd.date_range(startAt, periods=n_periods, freq="Y")
 dt_index_epochs = datetime_index.astype(int) // 10 ** 9
 dt_index = dt_index_epochs.astype("U10")

# For the Index of each CSA
 for idx, csa in rdf.iterrows():
     df = pd.DataFrame( { "color": csa.values[1:-1] }, index=dt_index, )
     styledata[idx] = df
 
 max_color, min_color = 0, 0
 for country, data in styledata.items():
     max_color = max(max_color, data["color"].max())
     min_color = min(max_color, data["color"].min())
 
 cmap = linear.PuRd_09.scale(min_color, max_color)
 def norm(x): return (x - x.min()) / (x.max() - x.min())
 for country, data in styledata.items():
     data["color"] = data["color"].apply(cmap)
     data["opacity"] = 1
 
 styledict = { str(country): data.to_dict(orient="index") for country, data in styledata.items() }
 
 # { CSA : { timestamp: {color: value, opacity:value } }, 
 #    CSA : { timestamp: {color: value, opacity:value } }, 
 #    ... 
 # }

from folium.plugins import TimeSliderChoropleth
 
 m = folium.Map([39.28759453969165, -76.61278931706487], width='50%', height='50%', zoom_start=12)
 g = TimeSliderChoropleth( rdf.to_json(), styledict=styledict, ).add_to(m)
 m