工具投影没反应_从空间格网开始，谈投影坐标，空间索引与快速空间连接

import geopandas as gpd
import pandas as pd
import numpy as np
from geopandas import GeoDataFrame
from shapely.geometry import Point
import pyproj
import utm
import math
from geopandas.tools import sjoin
import datetime
import os
import matplotlib.pyplot as plt
import mapclassify
import scipy.spatial# Get the UTM code
def convert_wgs_to_utm(lon, lat):utm_band &＃61; str((math.floor((lon &＃43; 180) / 6) % 60) &＃43; 1)if len(utm_band) &＃61;&＃61; 1:utm_band &＃61; &＃39;0&＃39; &＃43; utm_bandif lat >&＃61; 0:epsg_code &＃61; &＃39;326&＃39; &＃43; utm_band # lat>0: N;else:epsg_code &＃61; &＃39;327&＃39; &＃43; utm_bandreturn epsg_code# Read PGCOUNTY
# &＃39;D:COVID19RNN_COVID_DataShapefileNational_County_4326.shp&＃39;
poly &＃61; gpd.GeoDataFrame.from_file(r&＃39;D:COVID19RNN_COVID_DataShapefileNational_County_4326.shp&＃39;, crs&＃61;{&＃39;init&＃39;: &＃39;EPSG:4326&＃39;})
# Get the bounds
poly_bound &＃61; poly.geometry.bounds
# Get the utm code
utm_code &＃61; convert_wgs_to_utm(poly_bound[&＃39;minx&＃39;][0], poly_bound[&＃39;miny&＃39;][0])
# Project the lat/lon to meters
transformed_min &＃61; pyproj.transform(pyproj.Proj(init&＃61;&＃39;EPSG:4326&＃39;), pyproj.Proj(init&＃61;&＃39;EPSG:{0}&＃39;.format(utm_code)),poly_bound[&＃39;minx&＃39;][0], poly_bound[&＃39;miny&＃39;][0]) # lon lat
transformed_max &＃61; pyproj.transform(pyproj.Proj(init&＃61;&＃39;EPSG:4326&＃39;), pyproj.Proj(init&＃61;&＃39;EPSG:{0}&＃39;.format(utm_code)),poly_bound[&＃39;maxx&＃39;][0], poly_bound[&＃39;maxy&＃39;][0]) # lon lat
# Make a grid
Buffer_Value &＃61; 100
xmin, xmax, ymin, ymax &＃61; transformed_min[0], transformed_max[0], transformed_min[1], transformed_max[1]
xx, yy &＃61; np.meshgrid(np.arange(xmin, xmax, Buffer_Value), np.arange(ymin, ymax, Buffer_Value)) # meter
# Output
Grid_Point &＃61; GeoDataFrame({&＃39;geometry&＃39;: [Point(x, y) for x, y in zip(xx.flatten(), yy.flatten())]})
Grid_Point.crs &＃61; {&＃39;init&＃39;: &＃39;EPSG:{0}&＃39;.format(utm_code)}
Grid_Point &＃61; Grid_Point.to_crs({&＃39;init&＃39;: &＃39;EPSG:4326&＃39;})
# Within
pointInPolys &＃61; sjoin(Grid_Point, poly, how&＃61;&＃39;inner&＃39;, op&＃61;&＃39;within&＃39;).reset_index(drop&＃61;True)
pointInPolys .to_file("result2.shp") # output shp to test
# Buffer
Buffer_S &＃61; pointInPolys.buffer(Buffer_Value / 2, cap_style&＃61;3).reset_index(drop&＃61;True) # cap_style&＃61;3
Buffer_S &＃61; Buffer_S.reset_index()
Buffer_S.columns &＃61; [&＃39;B_ID&＃39;, &＃39;geometry&＃39;]
Buffer_S &＃61; gpd.GeoDataFrame(Buffer_S, geometry&＃61;&＃39;geometry&＃39;)
Buffer_S.crs &＃61; {&＃39;init&＃39;: &＃39;EPSG:{0}&＃39;.format(utm_code)}
Buffer_S &＃61; Buffer_S.to_crs({&＃39;init&＃39;: &＃39;EPSG:4326&＃39;})

# The points of grids
pointInPolys &＃61; pointInPolys.to_crs({&＃39;init&＃39;: &＃39;EPSG:{0}&＃39;.format(utm_code)})
pointInPolys[&＃39;lon&＃39;] &＃61; pointInPolys.geometry.apply(lambda p: p.x)
pointInPolys[&＃39;lat&＃39;] &＃61; pointInPolys.geometry.apply(lambda p: p.y)
ref_points &＃61; pointInPolys[[&＃39;lon&＃39;, &＃39;lat&＃39;]].values# The points we need to join with spatial grids
table &＃61; pd.read_pickle(&＃39;MDPG_device_new_data_%04s&＃39; % date)
# table &＃61; table.sample(10000)
table[&＃39;Time_index&＃39;] &＃61; round((table[&＃39;DateTime&＃39;] - Today).dt.total_seconds() / (Time_Bandwith * 60))
points &＃61; np.dstack(pyproj.transform(pyproj.Proj(init&＃61;&＃39;EPSG:4326&＃39;), pyproj.Proj(init&＃61;&＃39;EPSG:{0}&＃39;.format(utm_code)),table[&＃39;Longitude&＃39;].values, table[&＃39;Latitude&＃39;].values))[0] # lon lat
dist_ckd1, indexes_ckd1 &＃61; scipy.spatial.cKDTree(ref_points).query(points)
table[&＃39;Min_Distance&＃39;] &＃61; dist_ckd1
table[&＃39;B_ID&＃39;] &＃61; indexes_ckd1
# plt.plot(table[&＃39;Min_Distance1&＃39;] - table[&＃39;Min_Distance&＃39;], &＃39;o&＃39;)
table &＃61; table[table[&＃39;Min_Distance&＃39;] <&＃61; Buffer_Value / 2]

TjInPolys &＃61; gpd.GeoDataFrame(table, geometry&＃61;gpd.points_from_xy(x&＃61;table.Longitude, y&＃61;table.Latitude))
TjInPolys.crs &＃61; {&＃39;init&＃39;: &＃39;EPSG:4326&＃39;}
TjInPolys &＃61; sjoin(TjInPolys, Buffer_S, how&＃61;&＃39;inner&＃39;, op&＃61;&＃39;within&＃39;).reset_index(drop&＃61;True)
print(len(TjInPolys))
TjInPolys &＃61; TjInPolys.merge(table, on&＃61;[&＃39;Cuebiq_ID&＃39;, &＃39;DateTime&＃39;], how&＃61;&＃39;outer&＃39;)
plt.plot(TjInPolys[&＃39;B_ID&＃39;], TjInPolys[&＃39;B_ID_right&＃39;], &＃39;o&＃39;, alpha&＃61;0.1,color&＃61;&＃39;k&＃39;)

Encounter &＃61; table.groupby([&＃39;Time_index&＃39;, &＃39;B_ID&＃39;])[&＃39;Cuebiq_ID&＃39;].nunique().reset_index()
# Count and output
Slice &＃61; Encounter.groupby([&＃39;index_right&＃39;]).sum()[&＃39;Cuebiq_ID&＃39;].reset_index()
Slice.columns &＃61; [&＃39;B_ID&＃39;, &＃39;Count&＃39;]
Buffer_S_N &＃61; Buffer_S.merge(Slice, on&＃61;&＃39;B_ID&＃39;, how&＃61;&＃39;left&＃39;)
Buffer_S_N &＃61; Buffer_S_N.fillna(0)
fig, ax &＃61; plt.subplots(nrows&＃61;1, ncols&＃61;1, figsize&＃61;(12, 6))
# scheme &＃61; mapclassify.Quantiles(Buffer_S_N[&＃39;Count&＃39;], k&＃61;6)
Buffer_S_N.plot(column&＃61;&＃39;Count&＃39;, ax&＃61;ax, legend&＃61;True, scheme&＃61;&＃39;quantiles&＃39;, k&＃61;5,legend_kwds&＃61;dict(loc&＃61;4, frameon&＃61;False), linewidth&＃61;0.1, edgecolor&＃61;&＃39;white&＃39;,cmap&＃61;&＃39;OrRd&＃39;)