def draw_umap(data,n_neighbors=15, min_dist=0.1, n_components=2, metric='euclidean', title='' ,cc=0,rstate=42,dens=False): import umap import matplotlib.pyplot as plt fit = umap.UMAP( n_neighbors=n_neighbors, min_dist=min_dist, n_components=n_components, metric=metric, random_state=rstate, verbose=True, densmap=dens ) u = fit.fit_transform(data); plt.figure(figsize=(6, 5)) if n_components == 2: plt.scatter(u[:,0], u[:,1], c=cc,s=3,cmap=plt.cm.seismic) plt.clim(-5,5) plt.colorbar() plt.title(title, fontsize=18) return u; def TestHetero(DB1,DB2,UMAPMRK,CNum=20000,rs=42,fname=None): from scipy.spatial import distance from scipy.spatial import ConvexHull, convex_hull_plot_2d import pandas as pd import umap import numpy as np import matplotlib.pyplot as plt CAll=pd.concat([ DB1.sample(CNum,random_state=rs), DB2.sample(CNum,random_state=rs) ]).copy() X_2d=draw_umap(CAll[UMAPMRK],cc=CAll['H4'],min_dist=0.001,n_neighbors=60,rstate=42) L1=DB1.iloc[0].Line L2=DB2.iloc[0].Line print(L1,L2) m=CAll['Line']==L1 plt.scatter(X_2d[m,0],X_2d[m,1],c='b',s=1,label=L1); plt.scatter(X_2d[~m,0],X_2d[~m,1],c='r',s=1,label=L2); plt.legend(markerscale=10) if fname is not None: plt.savefig(fname) xmax=X_2d[:,0].max() xmin=X_2d[:,0].min() ymax=X_2d[:,0].max() ymin=X_2d[:,0].min() mx=np.max([xmax,ymax]) mn=np.min([xmin,ymin]) m=CAll.Line==L1 b=np.linspace(round(mn,0)-1,round(mx,0)+1,50) A,_,_=np.histogram2d(X_2d[m,0],X_2d[m,1],bins=b) # plt.figure(figsize=(5,5)) # plt.imshow(A>0) DD=distance.cdist(X_2d[m],X_2d[m]).flatten() hull1 = ConvexHull(X_2d[m]) print(L1," Local: ",(A>0).sum()," Global: ",np.round(np.quantile(DD,0.95),2)," Area: ",hull1.volume) m=CAll.Line==L2 b=np.linspace(round(mn,0)-1,round(mx,0)+1,50) A,_,_=np.histogram2d(X_2d[m,0],X_2d[m,1],bins=b) # plt.figure(figsize=(5,5)) # plt.imshow(A>0) DD=distance.cdist(X_2d[m],X_2d[m]).flatten() hull2 = ConvexHull(X_2d[m]) print(L2," Local: ", (A>0).sum()," Global: ",np.round(np.quantile(DD,0.95),2)," Area: ",hull2.volume) def TestHeteroSingle(DB1,UMAPMRK,CNum=20000,fname=None): from scipy.spatial import distance from scipy.spatial import ConvexHull, convex_hull_plot_2d import pandas as pd import umap import numpy as np import matplotlib.pyplot as plt CAll=DB1.sample(CNum,random_state=42,replace=False).copy() X_2d=draw_umap(CAll[UMAPMRK],cc=CAll['H4'],min_dist=0.001,n_neighbors=60,rstate=42) plt.scatter(X_2d[:,0],X_2d[:,1],c='b',s=1); if fname is not None: plt.savefig(fname) xmax=X_2d[:,0].max() xmin=X_2d[:,0].min() ymax=X_2d[:,0].max() ymin=X_2d[:,0].min() mx=np.max([xmax,ymax]) mn=np.min([xmin,ymin]) b=np.linspace(round(mn,0)-1,round(mx,0)+1,50) A,_,_=np.histogram2d(X_2d[:,0],X_2d[:,1],bins=b) # plt.figure(figsize=(5,5)) # plt.imshow(A>0) DD=distance.cdist(X_2d,X_2d).flatten() hull1 = ConvexHull(X_2d) print("Local: ",(A>0).sum()," Global: ",np.round(np.quantile(DD,0.95),2)," Area: ",hull1.volume) def TestHeteroSingle2D(X,MX,MN): b=np.linspace(round(MN,0)-1,round(MX,0)+1,50) A,_,_=np.histogram2d(X[:,0],X[:,1],bins=b) # plt.figure(figsize=(5,5)) # plt.imshow(A>0) DD=distance.cdist(X,X).flatten() hull1 = ConvexHull(X) print("Local: ",(A>0).sum()," Global: ",np.round(np.quantile(DD,0.95),2)," Area: ",hull1.volume) def HeteroMoran(DF, MRK, gridsize=10,k_min=4,k_max=5): import matplotlib.pyplot as plt from pointpats import PointPattern, PoissonPointProcess, as_window import pointpats from pysal.viz import splot from tqdm import tqdm import xarray # Surface data manipulation from pysal.explore import esda # Exploratory Spatial analytics from pysal.lib import weights from splot.esda import plot_moran import mapclassify as mc import numpy as np import geopandas as gpd import pandas as pd import matplotlib.pyplot as plt from shapely.geometry import Point, Polygon from matplotlib.patches import RegularPolygon import numpy as np, shapely.vectorized as sv from esda.geary import Geary MI=np.zeros((2,k_max-k_min)); GI=np.zeros((2,k_max-k_min)); def moran_by_col(y): mo = esda.Moran(y, w=w) mo_s = pd.Series( {"I": mo.I, "I-P value": mo.p_sim}, ) return mo_s DB=DF.copy() DB['Hex']=-1 DB.reset_index(inplace=True) Lns=DB.Line.unique() PP_NDD=[] X_2d=draw_umap(DB[MRK],cc=DB[DB.columns[0]],min_dist=0.001,n_neighbors=60,rstate=42) for j,L in enumerate(Lns): m=DB.Line==L pp=PointPattern(X_2d[m]); PP_NDD.append(pp.nnd.flatten()) print(L) plt.figure(figsize=(5,5)) A=plt.hexbin(X_2d[m,0],X_2d[m,1],gridsize=gridsize,cmap=plt.cm.seismic); hex_polys = A.get_paths()[0].vertices hex_array = [] for xs,ys in A.get_offsets(): hex_x = np.add(hex_polys[:,0], xs) hex_y = np.add(hex_polys[:,1], ys) hex_array.append(Polygon(np.vstack([hex_x, hex_y]).T)) Q=A.get_offsets() crs = {'init': 'epsg:4326'} hex_grid = gpd.GeoDataFrame({'geometry':hex_array, #'Cent':[Point(xy) for xy in Q] }, crs=crs) for i,G in tqdm(zip(hex_grid.index,hex_grid.geometry)): idx=sv.contains(G, x=X_2d[m,0], y=X_2d[m,1]) DB[m].loc[idx,'Hex']=i hex_grid.loc[i,'Num']=idx.sum() hex_grid[hex_grid.Num>0].plot(cmap=plt.cm.seismic,legend=True,column='Num') plt.show() for ck,k in enumerate(range(k_min,k_max)): wk = weights.distance.KNN.from_dataframe(hex_grid,k=k) moran = esda.moran.Moran(hex_grid["Num"], wk) MII=round(moran.I, 3) # print(k,MII) # print(j,ck) MI[j,ck]=MII geary=Geary(hex_grid["Num"], wk) GII=round(geary.C,3) GI[j,ck]=GII return MI,GI,PP_NDD # moran_stats = hex_grid[['Num']].apply(moran_by_col, axis=0).T # print(k,moran_stats) # moran_stats = DF[NamesAll].apply(moran_by_col, axis=0).T