<jupyter_start><jupyter_text>To install geopandas: conda install geopandasAlso needed: conda install descartesIn new environment:conda create gis python==3.8.3conda activate gisconda install -y -c conda-forge geopandas jupyterLabconda remove --force pyprojIn anaconda prompt: conda activate gisconda listWHAT I FINALLY DID:conda create --name thesis geopandasconda activate thesisconda listconda install -c conda-forge foliumconda install rasterioconda install -c conda-forge earthengine-api<jupyter_code>import os 
os.chdir('C:/Users/Lenovo/Desktop/Master_BAOR/MSc Thesis/GitHub/DiCRA_clone')
#conda install -y -c conda-forge geopandas jupyterlab
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point
import pyproj
import pprint
import folium
from IPython.display import display
import rasterio
#from pyproj import CRS
#crs=CRS('EPSG:4326')
#from pyproj import CRS
#import pyproj 

#proj_4326 = pyproj.Proj("epsg:4326")
#proj_3857 = pyproj.Proj("epsg:3857")
#pyproj.Proj('+init=espg: 3857')

boundaries_df = gpd.read_file('dicra/src/data_preprocessing/tsdm/District_Boundary.shp')
print(boundaries_df.head())

NDV_index = gpd.read_file('C:/Users/Lenovo/Desktop/Master_BAOR/MSc Thesis/dowloaded from DiCRA/NDVI/NDVI_2021/VECTOR/DISTRICT/28-08-2021.geojson')
NDV_index.head()

soil_df = gpd.read_file('C:/Users/Lenovo/Desktop/Master_BAOR/MSc Thesis/dowloaded from DiCRA/SOILM/SOILM_2022/VECTOR/DISTRICT/23-01-2022.geojson')
soil_df.head()

market_yards = pd.read_csv('C:/Users/Lenovo/Desktop/Master_BAOR/MSc Thesis/dowloaded from DiCRA/CP/CP_2022/market_yard_data.csv')
#print(market_yards.head())

relative_wealth = pd.read_csv('C:/Users/Lenovo/Desktop/Master_BAOR/MSc Thesis/GitHub/DiCRA_clone/dicra/analytics/notebooks/relative_wealth_index/data/ind_pak_relative_wealth_index.csv')
#print(relative_wealth.head())
# Viewing a geometry (change the number to see some other district)
boundaries_df.loc[0, 'geometry']
# Gives all latitude, longitude pairs that define the polygon as above
#print(boundaries_df.loc[0, 'geometry'])
boundaries_df.plot()
plt.show()
# The colors are based on the size of the area for each district
boundaries_df.plot(column = 'Area', legend = True)
plt.show()
# Each district is given a color
boundaries_df.plot(column = 'Dist_Name', legend = True)
plt.show()<jupyter_output><empty_output><jupyter_text>The following code can be used if we would like to plot some other data points over the polygons<jupyter_code># Plot of the market yards within Telangana
boundaries_df.plot(column = 'Dist_Name', legend = False, cmap = 'Set2')
plt.scatter(market_yards.longitude, market_yards.latitude, marker = 'p', c = 'red', edgecolor = 'darkred')
plt.title('Unknown')
plt.show()
#GeoJSON and plotting with geopandas
leg_kwds = {'title': 'District Name', 'loc': 'upper left', 'bbox_to_anchor': (1,1.03), 'ncol': 3}
boundaries_df.plot(column = 'Dist_Name', cmap='Set3', legend=True, legend_kwds=leg_kwds)
plt.title('District names')
plt.show()<jupyter_output><empty_output><jupyter_text>Example colormaps:Pastel1, Pastel2, Paired, Accent, Dark2, Set1, summer<jupyter_code>soil_df.plot(column = 'Dist_Name', cmap='Set2', legend=True)
plt.show()
# Creating a GeoDataFrame from a DataFrame
# from shapely.geometry import Point

market_yards['geometry'] = market_yards.apply(lambda x: Point((x.longitude, x.latitude)), axis = 1)
print(market_yards.head(3))

# Now we have a geometry column, we can build a GeoDataFrame from this

marketyards_geo = gpd.GeoDataFrame(market_yards, geometry = market_yards.geometry, crs = {'init': 'epsg:4326'})
#marketyards_geo.crs = "EPSG:4326"
marketyards_geo.head(3)


#market_yards_crs = {'init:': 'epsg:4326'}
#marketyards_geo = gpd.GeoDataFrame(market_yards, crs = market_yards_crs, geometry = market_yards.geometry)


# So we just changed the datatype from a DataFrame to a GeoDataFrame
#gdf = gdf.set_crs(4326, allow_override=True)
#marketyards_geo.set_crs(epsg=4326)
# Decimal degrees are used to measure the distance from the reference points: the Prime Meridian and the Equator
# We can convert the geometry to measure distance in meters, using to_crs()
# The original longitude and latitude columns remain in decimal degree units
# Convert geometry from decimal degrees to meters
marketyards_geo.crs = "EPSG:4326"
marketyards_geo.geometry = marketyards_geo.geometry.to_crs(epsg = 3857)
marketyards_geo.head(2)

#Convert geometry from decimal degrees to meters
#fires_data_geo.geometry = fires_data_geo.geometry.to_crs(epsg = 3857)
#fires_data_geo.head()

#art['geometry'] = art.apply(lambda x: Point(float(x.lng), float(x.lat)), axis=1)
#art_geo = gpd.GeoDataFrame(art, crs = neighborhoods.crs, geometry = art.geometry)
#print(type(art_geo))<jupyter_output><empty_output><jupyter_text>Spatial joinsgpd.sjoin(name1_gdf, name2_gdf, op = ) op = intersects, contains or within<jupyter_code># Returns all observations where the region intersects the points
intersect_gdf = gpd.sjoin(boundaries_df, marketyards_geo, op = 'intersects') # First input POINT, second input POLYGON geometry
print('Points which are completely within this region', intersect_gdf.shape[0]) 
# Number of market yards that intersect a district
# Number of market yards within a district
# Returns all observations where the region completely contains points (so not on boundary anymore)
contains_gdf = gpd.sjoin(boundaries_df, marketyards_geo, op = 'contains')
print('Points which are completely within this region', contains_gdf.shape[0]) 
# Number of districts that contain market yards 
# No cases where the region is completely in a point, but reverse order to get points within the region
within_gdf = gpd.sjoin(marketyards_geo, boundaries_df, op = 'within')
# Adjust the number to get it for another region
print('Points which are completely within this region', within_gdf.shape[0]) 
# Number of market yards within a district
within_gdf.head()
# Index right: the position of the joined data in the original right-hand GDF (here boundaries_df)
# Row index values for the left-hand GDF are retained as the new row indices for within_gdf
print(within_gdf.columns)
print(within_gdf.shape)<jupyter_output><empty_output><jupyter_text>within_gdf.district_left = boundaries_districtwithin_gdf.district_right = market_yards_districtwithin_gdf[['boundaries_district', 'market_yards_district']].groupby('market_yards_district').agg('count').sort_values('boundaries_district', ascending = False) Get name and title from neighborhood_art and group by nameneighborhood_art_grouped = neighborhood_art[['name', 'title']].groupby('name') Aggregate the grouped data and count the artworks within each polygonprint(neighborhood_art_grouped.agg('count').sort_values(by = 'title', ascending = False)) Create urban_art from neighborhood_art where the neighborhood name is Urban Residentsurban_art = neighborhood_art.loc[neighborhood_art.name == "Urban Residents"] Get just the Urban Residents neighborhood polygon and save it as urban_polygonurban_polygon = neighborhoods.loc[neighborhoods.name == "Urban Residents"] Plot the urban_polygon as ax ax = urban_polygon.plot(color = 'lightgreen') Add a plot of the urban_art and show iturban_art.plot( ax = ax, column = 'type', legend = True);plt.show() You can think of a GeoSeries as the geometry column of a GeoDataFrame.GeoSeries.area: returns the area of each geometry in GeoSeriesGeoSeries.centroid: returns the center point of each geometry in a GeoSeriesGeoSeries.distance(other): returns the minimum distance to a location specified using the 'other'argument<jupyter_code># Area of first polygon in districts (units correspond to units given)
print(boundaries_df.head(4))
print(boundaries_df.geometry[0].area)
print('There are ', boundaries_df.shape[0], ' districts.')
# [0] for rows, [1] for columns
# Calculate area of each district
# EPSG: 4326 uses decimal degrees for distance (so areas are decimal degrees squared)
#boundaries_df.geometry.to_crs(4326)
district_area = boundaries_df.geometry.area
# Print the areas and crs used
print(district_area.sort_values(ascending = False))
print(boundaries_df.crs)<jupyter_output><empty_output><jupyter_text>NOTE: The coordinate references system decides upon the distance units<jupyter_code># Now, kilometers squared. Use CRS that uses meters for distance. Later convert to km.
boundaries_df_3857 = boundaries_df.to_crs(epsg = 3857)
sqm_to_sqkm = 10**6
district_area_km = boundaries_df_3857.geometry.area / sqm_to_sqkm
print(district_area_km.sort_values(ascending = False))
print(boundaries_df_3857.crs)
# NOTE: Use 3857 instead of 4326 (projected instead of geo)
boundaries_df_3857['center'] = boundaries_df_3857.geometry.centroid
# Create GDF with districts and centers
part = ['Dist_Name', 'center']
boundaries_df_3857_centers = boundaries_df_3857[part]
boundaries_df_3857_centers.head(3)
district_Maha = boundaries_df_3857.loc[boundaries_df_3857.Dist_Name == 'Mahabubabad']
district_Maha.head()
marketyards_geo.head() # GDF for the market yards (see previous derivation)
marketyards_geo_3857 = marketyards_geo.to_crs(epsg = 3857)
# Spatial join market yards within district Mahabubabad
# Finds market yards that are within district Mahabubabad
yards_in_Maha = gpd.sjoin(marketyards_geo_3857, district_Maha, op = 'within')
yards_in_Maha.shape # There are 7 schools in district Mahabubabad
# Distances dictionary that includes the distance from the center of the district to each market yard in the district 
distances = {}
for row in yards_in_Maha.iterrows():
    vals = row[1]
    key = vals['name']
    ctr = vals['center']
    distances[key] = vals['geometry'].distance(ctr)
pprint.pprint(distances)# NOTE: Use 3857 instead of 4326 (projected instead of geo)
boundaries_df_3857['center'] = boundaries_df_3857.geometry.centroid
# Create GDF with districts and centers
part = ['Dist_Name', 'center']
boundaries_df_3857_centers = boundaries_df_3857[part]
boundaries_df_3857_centers.head(3)

# NOTE: Use 3857 instead of 4326 (projected instead of geo)
district_Maha = boundaries_df_3857.loc[boundaries_df_3857.Dist_Name == 'Mahabubabad']
center_Maha = district_Maha.geometry.centroid

# Plot the urban_poly_3857 as ax and add the center point
ax = district_Maha.plot(color = 'lightgreen')
center_Maha.plot(ax = ax, color = 'black')
plt.xticks(rotation = 45)
# Show the plot
plt.show()
district_Maha = boundaries_df_3857.loc[boundaries_df_3857.Dist_Name == 'Mahabubabad']
center_Maha = district_Maha.geometry.centroid
# Finds market yards that are within district Mahabubabad
yards_in_Maha = gpd.sjoin(marketyards_geo_3857, district_Maha, op = 'within')
# Plot the urban_poly_3857 as ax and add the center point
ax = district_Maha.plot(color = 'lightgreen')
yards_in_Maha.plot(ax = ax, color = 'black')
plt.xticks(rotation = 45)
# Show the plot
plt.show()
# Create art_dist_meters using art and the geometry from art
district_exercise = boundaries_df.loc[boundaries_df.Dist_Name == 'Mahabubabad']
boundaries_exercise = gpd.GeoDataFrame(district_exercise, geometry = district_exercise.geometry, crs = {'init': 'epsg:4326'})
print(boundaries_exercise.head(2))
#Use 'GeoDataFrame.set_crs(crs, allow_override=True)' to overwrite CRS or 'GeoDataFrame.to_crs(crs)' to reproject geometries.
# Set the crs of art_dist_meters to use EPSG:3857
boundaries_exercise.geometry = boundaries_exercise.geometry.to_crs(epsg = 3857)
print(boundaries_df.head(2))

# Add a column to art_meters, center
boundaries_exercise['center'] = center_Maha
# Folium wants coordinates as an array with latitude first
boundaries_df.head()
district_Maha_fol = boundaries_df.loc[boundaries_df.Dist_Name == 'Mahabubabad']
print(district_Maha_fol.head())
district_Maha_fol['center'] = district_Maha_fol.geometry.centroid
type(district_Maha_fol.center)
# Create an array of the center of the district
center_Maha = district_Maha_fol['center'].iloc[0]
print(center_Maha)
print(type(center_Maha))
# Reverse the order of the coordinate pair for folium
district_center = [center_Maha.y, center_Maha.x]
print(district_center)
# Create a folium map centered on the Maha district
district_map = folium.Map(location = district_center, zoom_start = 12)
# Add the outline of the Maha district
folium.GeoJson(district_Maha_fol.geometry).add_to(district_map)
# Display the resulting map
display(district_map)

for row in yards_in_Maha.iterrows():
    row_values = row[1]
    #print(row_values)
    location = [row_values['latitude'], row_values['longitude']]
    print(location)
    # Create popups inside the loop you built to create the markers
    # Here we bould a popup message that includes the name
    popup = popup = '<strong>' + row_values['name'] + '</strong>'
    marker = folium.Marker(location = location, popup = popup)
    marker.add_to(district_map)
display(district_map)

# Other popup example
#popup = (str(row_values['title']) + ': ' + str(row_values['desc'])).replace("'", "`")

# Two common issues:
# missing values and single/double quotes in the data itself

# Replace Nan and ' values in description
#urban_art.desc.fillna('', inplace = True)
#urban_art.desc = urban_art.desc.str.replace("'", "`")
boundaries_df.Dist_Name = boundaries_df.Dist_Name.str.replace(" ", "")
marketyards_geo.district = marketyards_geo.district.str.replace(" ", "")
marketyards_geo.head()

# Count the number of yards within each district
yards_counts = marketyards_geo.groupby(['district']).size()
yards_counts
boundaries_df.head()
# Add these counts to the GDF
yards_counts_df = yards_counts.to_frame()
yards_counts_df = yards_counts_df.reset_index() # .reset_index(inplace=True)
yards_counts_df.columns = ['Dist_Name', 'yards_count'] # Assign a name to each column
districts_with_counts = pd.merge(boundaries_df, yards_counts_df, on = 'Dist_Name')
print(districts_with_counts.head(2))
# Normalize
# Store the area for each district
districts_with_counts['area'] = districts_with_counts.geometry.area
# Divide the counts by the areas (so we create a density column)
districts_with_counts['yards_density'] = districts_with_counts.apply(
    lambda row: row.yards_count/row.area, axis = 1)
districts_with_counts.head(2)
# 4326: uses deciaml degrees for distance
# 3857: uses meters for distance
print(boundaries_df.crs) # Check whether we have crs 3857
# If not, use the following: .to_crs(epsg = 3857)
print(boundaries_df_3857.crs) # Check whether we have crs 3857
sqm_to_sqkm = 10**6 # define a variable for m^2 to km^2
boundaries_df_3857['area'] = boundaries_df_3857.geometry.area / sqm_to_sqkm
boundaries_df_3857.head(5)
# Now also latitude and longitude measures are in meters instead of decimal degrees (traditional)
# Reconvert crs back to 4326
boundaries_df_back_4326 = boundaries_df_3857.to_crs(epsg = 4326)
print(boundaries_df_back_4326.crs)
# Now, area is in km^2 and geometry in decimal degrees
print(boundaries_df_back_4326.head(2))
marketyards_geo = marketyards_geo.to_crs(epsg = 4326)
print(marketyards_geo.crs)
# Spatial join to get the districts that contain the yards
yards_in_districts = gpd.sjoin(boundaries_df_back_4326, marketyards_geo, op = 'contains')
print(yards_in_districts.head(2))
#print(boundaries_df_back_4326.head(2))
#print(marketyards_geo.head(2))

# Aggregate to get counts
yards_counts_choropleths = yards_in_districts.groupby(['district']).size()
#print(yards_counts_choropleths)

# Convert this to a df
yards_counts_choropleths_df = yards_counts_choropleths.to_frame()
yards_counts_choropleths_df = yards_counts_choropleths_df.reset_index(level=0)
yards_counts_choropleths_df.columns = ['Dist_Name', 'yards_count']
yards_counts_choropleths_df.head(2)

# merge
# Because we list the GeoDF first, the resulting merged object is also a GeoDF. Otherwise DF.
districts_with_counts = pd.merge(boundaries_df_back_4326, yards_counts_choropleths_df, on = 'Dist_Name')
districts_with_counts.head()

# Create the yards_density
districts_with_counts['yards_density'] = districts_with_counts.apply(lambda row: row.yards_count / row.area, axis = 1)
districts_with_counts.head()

# Create choropleth plot
districts_with_counts.plot(column = 'yards_density', cmap = 'BuGn', edgecolor = 'black', legend = True)
plt.title('Market yards per kilometers squared')
plt.xlabel('longitude')
plt.ylabel('latitude')
plt.show()<jupyter_output><empty_output><jupyter_text>Arguments of the folium choroplethgeo_data: the source data for the polygons (geojson file or a GeoDF)name: the name of hte geometry column (or geojson property) for the polygonsdata: the source DF or Series for the normalized datacolumns: a list of columns: one that corresponds to the polygons and one that has the value to plotkey_on: a GeoJSON variable to bind the data to (always start with 'feature')fill_color: polygon fill color (default blue)fill_opacity: range between 0 (transparent) and 1 (completely opaque)line_color: color of polygon border lines (defaults to black)line_opacity: range between 0 (transparent) and 1 (completely opaque)<jupyter_code># Center point and map
center = [center_Maha.y, center_Maha.x]
m = folium.Map(location = center, zoom_start = 10)
print(districts_with_counts.head())
cds = ColumnDataSource(districts_with_counts.drop(columns=['geometry']))
# Define a choropleth layer for the map
m.choropleth(
    geo_data = districts_with_counts,
    name = 'geometry',
    data = districts_with_counts,
    columns = ['Dist_Name', 'yards_density'],
    key_on = 'feature.properties.Dist_Name',
    fill_color = 'YlGn',
    fill_opacity = 0.75,
    line_opacity = 0.5,
    legend_name = 'Yards per km squared by District'
)

# Add layer control and display
folium.LayerControl().add_to(m)
display(m)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09775414380268145, 'fraction_numerical': 0.021779713752333542, 'mean_word_length': 6.5773681954564935, 'pattern_counts': {'":': 0, '<': 31, '<?xml version=': 0, '>': 31, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '410304', 'license': 'MIT', 'n_tokens_mistral': 6305, 'n_tokens_neox': 5813, 'n_words': 1586, 'path': 'analytics/notebooks/crop_yield_prediction/data_exploration.ipynb', 'repo_name': 'UNDP-India/Data4Policy', 'size': '17677'}