onefan/UrbanFlow-DataAnalysis-Public

1 2import pandas as pd3import numpy as np4import math5import urllib6import io7from PIL import Image8import os9import sys10def searchfile(filename):11    fileroot = ''12    for i in sys.path:13        try:14            if filename in os.listdir(i):15                fileroot = i16        except:17            pass18    if fileroot == '':19        fileroot = sys.path[-1]20    return fileroot+'/'+filename21 22def set_mapboxtoken(mapboxtoken):23    filepath = searchfile('mapboxtoken.txt')24    f = open(filepath,mode = 'w')25    f.write(mapboxtoken)26    f.close()27    print('Success')28    29def read_mapboxtoken():30    filepath = searchfile('mapboxtoken.txt')31    try:32        f = open(filepath,mode = 'r')33        mapboxtoken = f.readline()34        f.close()35    except:36        raise Exception('Mapboxtoken not found, please use tbd.set_mapboxtoken() to set it first, see: https://transbigdata.readthedocs.io/en/latest/plot_map.html') 37    return mapboxtoken38 39 40def set_imgsavepath(imgsavepath):41    filepath = searchfile('imgsavepath.txt')42    f = open(filepath,mode = 'w')43    f.write(imgsavepath)44    f.close()45    print('Success')46def read_imgsavepath():47    filepath = searchfile('imgsavepath.txt')48    try:49        f = open(filepath,mode = 'r')50        imgsavepath = f.readline()51        f.close()52    except:53        raise Exception('Map base map storage path not found, please use tbd.set_imgsavepath() to set it first, see: https://transbigdata.readthedocs.io/en/latest/plot_map.html') 54    return imgsavepath55 56def deg2num(lat_deg, lon_deg, zoom):57    lat_rad = math.radians(lat_deg)58    n = 2.0 ** zoom59    xtile = int((lon_deg + 180.0) / 360.0 * n)60    ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)61    return (xtile, ytile)62 63def num2deg(xtile, ytile, zoom):64    n = 2.0 ** zoom65    lon_deg = xtile / n * 360.0 - 180.066    lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))67    lat_deg = math.degrees(lat_rad)68    return (lat_deg, lon_deg)69 70def getImageCluster( lon_deg,lat_deg, delta_long, delta_lat,zoom,printlog,imgsavepath,style = 4,access_token = ''):71    '''72    access_token - mapbox token73    '''74    if (style == 1 )|(style == 'streets' ):75        styleid = 'ckwinzgw581od14mpyfhka6nk'76        smurl = r'https://api.mapbox.com/styles/v1/ni1o1/'+styleid+r'/tiles/256/{0}/{1}/{2}?&access_token='+access_token77    if (style == 2 )|(style == 'outdoors' ):78        styleid = 'ckwinx7aj4y4a15p7ftfwq9dn'79        smurl = r'https://api.mapbox.com/styles/v1/ni1o1/'+styleid+r'/tiles/256/{0}/{1}/{2}?&access_token='+access_token80    if (style == 3 )|(style == 'satellite' ):81        styleid = 'cjv36cj9u4h1q1ftemjed4f2y'82        smurl = r'https://api.mapbox.com/styles/v1/ni1o1/'+styleid+r'/tiles/256/{0}/{1}/{2}?&access_token='+access_token83    if (style == 4 )|(style == 'light' ):84        styleid = 'ckwfx658z4dpb14ocnz6tky9d'85        smurl = r'https://api.mapbox.com/styles/v1/ni1o1/'+styleid+r'/tiles/256/{0}/{1}/{2}?&access_token='+access_token86    if (style == 5 )|(style == 'dark' ):87        styleid = 'cjetnd20i1vbi2qqxbh0by7p8'88        smurl = r'https://api.mapbox.com/styles/v1/ni1o1/'+styleid+r'/tiles/256/{0}/{1}/{2}?&access_token='+access_token89    if (style == 6 )|(style == 'light-ch' ):90        styleid = 'ckj9bhq7s9mvj19mq3e3fye35'91        smurl = r'https://api.mapbox.com/styles/v1/ni1o1/'+styleid+r'/tiles/256/{0}/{1}/{2}?&access_token='+access_token92    if (style == 7 )|(style == 'ice creem' ):93        styleid = 'cjv36iiz9243t1fo8mweb4z6r'94        smurl = r'https://api.mapbox.com/styles/v1/ni1o1/'+styleid+r'/tiles/256/{0}/{1}/{2}?&access_token='+access_token95    if (style == 8 )|(style == 'night' ):96        styleid = 'ck2o3fyvy0dch1cp6j2pkz2dv'97        smurl = r'https://api.mapbox.com/styles/v1/ni1o1/'+styleid+r'/tiles/256/{0}/{1}/{2}?&access_token='+access_token98    if (style == 9 )|(style == 'terrain' ):99        styleid = 'cjv36gyklf43q1fnuwibiuetl'100        smurl = r'https://api.mapbox.com/styles/v1/ni1o1/'+styleid+r'/tiles/256/{0}/{1}/{2}?&access_token='+access_token101    if (style == 10 )|(style == 'basic blue' ):102        styleid = 'ckwio2ze12fgk15p2alr5a4xj'103        smurl = r'https://api.mapbox.com/styles/v1/ni1o1/'+styleid+r'/tiles/256/{0}/{1}/{2}?&access_token='+access_token104 105    else:106        styleid = ''107    xmin, ymax =deg2num(lat_deg, lon_deg, zoom)108    xmax, ymin =deg2num(lat_deg + delta_lat, lon_deg + delta_long, zoom)109    110    def get_img(smurl,zoom, xtile, ytile,imgsize,imgsavepath):111        import os112        113        filename = str(style)+str(styleid)+'-'+str(zoom)+'-'+str(xtile)+'-'+str(ytile)+'-'+str(imgsize)+'.png'114        def savefig(filename,tile):115            try:116                if 'tileimg' in os.listdir(imgsavepath):117                    if filename in os.listdir(imgsavepath+'tileimg'):118                        pass119                    else:120                        tile.save(imgsavepath+'tileimg/'+filename)121                        if printlog:122                            print('figsaved:'+imgsavepath+'tileimg/'+filename)123                else:124                    os.mkdir(imgsavepath+'tileimg')125            except:126                pass127        def loadfig(filename):128            try:129                if 'tileimg' in os.listdir(imgsavepath):130                    if filename in os.listdir(imgsavepath+'tileimg'):131                        tile = Image.open(imgsavepath+'tileimg\\'+filename)132                        return tile133                    else:134                        return None135                else:136                    os.mkdir(imgsavepath+'tileimg')137                    return None138            except:139                return None140        tile = loadfig(filename)141        if tile is None:142            try:143                t = 0144                while t<10:145                    try:146                        imgurl=smurl.format(zoom, xtile, ytile)147                        #print("Opening: " + imgurl)148                        imgstr = urllib.request.urlopen(imgurl,timeout = 6).read()149                        tile = Image.open(io.BytesIO(imgstr))150                        savefig(filename,tile)151                        Cluster.paste(tile, box=((xtile-xmin)*imgsize ,  (ytile-ymin)*imgsize))152                        t = 10153                    except:154                        if printlog:155                            print('Get map tile failed, retry ',t)156                        t += 1157            except: 158                print("Couldn't download image")159                tile = None160        else:161            Cluster.paste(tile, box=((xtile-xmin)*imgsize ,  (ytile-ymin)*imgsize))162 163    imgsize = 256164    import threading165    threads = []166    Cluster = Image.new('RGB',((xmax-xmin+1)*imgsize-1,(ymax-ymin+1)*imgsize-1)) 167    for xtile in range(xmin, xmax+1):168        for ytile in range(ymin,  ymax+1):169            threads.append(threading.Thread(target=get_img,args = (smurl,zoom, xtile, ytile,imgsize,imgsavepath)))170    for t in threads:171        t.setDaemon(True)172        t.start()173    for t in threads:174        t.join()175    threads.clear()176    177    return Cluster178 179    180def plot_map(plt,bounds,zoom = 'auto',style=4,printlog = False):181    '''182    Plot the basemap183 184    Parameters185    -------186    plt : matplotlib.pyplot187        Where to plot188    bounds : List189        The drawing boundary of the base map, [lon1,lat1,lon2,lat2] (WGS84 coordinate system), where lon1 and lat1 are the coordinates of the lower left corner and lon2 and lat2 are the coordinates of the upper right corner190    zoom : number191        The larger the magnification level of the base map, the longer the loading time. Generally, the range for a single city is between 12 and 16192    printlog : bool193        Show log194    style : number195        The style of map basemap can be 1-10, as follows196    '''197    access_token = read_mapboxtoken()198    imgsavepath = read_imgsavepath()199    try:200        import os201        os.listdir(imgsavepath)202    except:203        print('imgsavepath do not exist, your tile map will not save')204    lon1 = bounds[0]205    lat1 =  bounds[1]206    lon2 =  bounds[2]207    lat2 =  bounds[3]208    if zoom == 'auto':209        zoom = 11-np.log(lon2-lon1)/np.log(2)210    zoom = min(18,int(zoom+0.5))211    a = getImageCluster(lon1, lat1, lon2-lon1,  lat2-lat1, zoom,style = style,printlog = printlog,imgsavepath = imgsavepath,access_token = access_token)212    x1, y1 =deg2num(lat1, lon1, zoom)213    x2, y2 =deg2num(lat2, lon2, zoom)214    x1,y1 = num2deg(x1, y1+1, zoom)215    x2,y2 = num2deg(x2+1, y2, zoom)216    plt.imshow(np.asarray(a),extent = (y1,y2,x1+0.00,x2+0.00))217 218 219def plotscale(ax,bounds,textcolor = 'k',textsize = 8,compasssize = 1,accuracy = 'auto',rect=[0.1,0.1],unit = "KM",style = 1,**kwargs):220    '''221    Add compass and scale for a map222 223    Parameters224    -------225    bounds : List226        The drawing boundary of the base map, [lon1,lat1,lon2,lat2] (WGS84 coordinate system), where lon1 and lat1 are the coordinates of the lower left corner and lon2 and lat2 are the coordinates of the upper right corner227    textsize : number228        size of the text229    compasssize : number230        size of the compass231    accuracy : number232        Length of scale bar (m)233    unit : str234        ‘KM’,’km’,’M’,’m’, the scale units235    style : number236        1 or 2, the style of the scale237    rect : List238        The approximate position of the scale bar in the figure, such as [0.9,0.9], is in the upper right corner239    '''240    import math241    lon1 = bounds[0]242    lat1 = bounds[1]243    lon2 = bounds[2]244    lat2 = bounds[3]245    latStart = min(lat1, lat2);246    lonStart = min(lon1, lon2);247    if accuracy == 'auto':248        accuracy = (int((lon2-lon1)/0.0003/1000+0.5)*1000)249    a,c=rect250    b = 1-a251    d = 1-c252    alon,alat = (b*lon1+a*lon2)/(a+b),(d*lat1+c*lat2)/(c+d)253 254    deltaLon = accuracy * 360 / (2 * math.pi * 6371004 * math.cos((lat1 + lat2) * math.pi / 360));255    256    from shapely.geometry import Polygon257    import geopandas as gpd258    if style == 1:259        scale = gpd.GeoDataFrame({'color':[(0,0,0),(1,1,1),(0,0,0),(1,1,1)],'geometry':260        [Polygon([(alon,alat),(alon+deltaLon,alat),(alon+deltaLon,alat+deltaLon*0.4),(alon,alat+deltaLon*0.4)]),261        Polygon([(alon+deltaLon,alat),(alon+2*deltaLon,alat),(alon+2*deltaLon,alat+deltaLon*0.4),(alon+deltaLon,alat+deltaLon*0.4)]),262        Polygon([(alon+2*deltaLon,alat),(alon+4*deltaLon,alat),(alon+4*deltaLon,alat+deltaLon*0.4),(alon+2*deltaLon,alat+deltaLon*0.4)]),263        Polygon([(alon+4*deltaLon,alat),(alon+8*deltaLon,alat),(alon+8*deltaLon,alat+deltaLon*0.4),(alon+4*deltaLon,alat+deltaLon*0.4)])264        ]})265        scale.plot(ax = ax,edgecolor= textcolor,facecolor = scale['color'],lw = 0.6,**kwargs)266 267        if (unit == 'KM')|(unit == 'km'):268            ax.text(alon+1*deltaLon,alat+deltaLon*0.5,str(int(1*accuracy/1000)),color = textcolor,fontsize = textsize,ha = 'center',va = 'bottom')269            ax.text(alon+2*deltaLon,alat+deltaLon*0.5,str(int(2*accuracy/1000)),color = textcolor,fontsize = textsize,ha = 'center',va = 'bottom')270            ax.text(alon+4*deltaLon,alat+deltaLon*0.5,str(int(4*accuracy/1000)),color = textcolor,fontsize = textsize,ha = 'center',va = 'bottom')271            ax.text(alon+8*deltaLon,alat+deltaLon*0.5,str(int(8*accuracy/1000)),color = textcolor,fontsize = textsize,ha = 'center',va = 'bottom')272            ax.text(alon+8.5*deltaLon,alat+deltaLon*0.5,unit,color = textcolor,fontsize = textsize,ha = 'left',va = 'top')273        if (unit == 'M')|(unit == 'm'):274            ax.text(alon+1*deltaLon,alat+deltaLon*0.5,str(int(1*accuracy)),color = textcolor,fontsize = textsize,ha = 'center',va = 'bottom')275            ax.text(alon+2*deltaLon,alat+deltaLon*0.5,str(int(2*accuracy)),color = textcolor,fontsize = textsize,ha = 'center',va = 'bottom')276            ax.text(alon+4*deltaLon,alat+deltaLon*0.5,str(int(4*accuracy)),color = textcolor,fontsize = textsize,ha = 'center',va = 'bottom')277            ax.text(alon+8*deltaLon,alat+deltaLon*0.5,str(int(8*accuracy)),color = textcolor,fontsize = textsize,ha = 'center',va = 'bottom')278            ax.text(alon+8.5*deltaLon,alat+deltaLon*0.5,unit,color = textcolor,fontsize = textsize,ha = 'left',va = 'top')279    if style == 2:280        scale = gpd.GeoDataFrame({'color':[(0,0,0),(1,1,1)],'geometry':281        [Polygon([(alon+deltaLon,alat),(alon+4*deltaLon,alat),(alon+4*deltaLon,alat+deltaLon*0.4),(alon+deltaLon,alat+deltaLon*0.4)]),282        Polygon([(alon+4*deltaLon,alat),(alon+8*deltaLon,alat),(alon+8*deltaLon,alat+deltaLon*0.4),(alon+4*deltaLon,alat+deltaLon*0.4)])283        ]})284        scale.plot(ax = ax,edgecolor= textcolor,facecolor = scale['color'],lw = 0.6,**kwargs)285 286        if (unit == 'KM')|(unit == 'km'):287            ax.text(alon+4*deltaLon,alat+deltaLon*0.5,str(int(4*accuracy/1000)),color = textcolor,fontsize = textsize,ha = 'center',va = 'bottom')288            ax.text(alon+8*deltaLon,alat+deltaLon*0.5,str(int(8*accuracy/1000)),color = textcolor,fontsize = textsize,ha = 'center',va = 'bottom')289            ax.text(alon+8.5*deltaLon,alat+deltaLon*0.5,unit,color = textcolor,fontsize = textsize,ha = 'left',va = 'top')290        if (unit == 'M')|(unit == 'm'):291            ax.text(alon+4*deltaLon,alat+deltaLon*0.5,str(int(4*accuracy)),color = textcolor,fontsize = textsize,ha = 'center',va = 'bottom')292            ax.text(alon+8*deltaLon,alat+deltaLon*0.5,str(int(8*accuracy)),color = textcolor,fontsize = textsize,ha = 'center',va = 'bottom')293            ax.text(alon+8.5*deltaLon,alat+deltaLon*0.5,unit,color = textcolor,fontsize = textsize,ha = 'left',va = 'top')294    deltaLon = compasssize*deltaLon295    alon = alon-deltaLon296    compass = gpd.GeoDataFrame({'color':[(0,0,0),(1,1,1)],'geometry':297    [Polygon([[alon,alat],[alon,alat+deltaLon],[alon+1/2*deltaLon,alat-1/2*deltaLon]]),298    Polygon([[alon,alat],[alon,alat+deltaLon],[alon-1/2*deltaLon,alat-1/2*deltaLon]])]})299    compass.plot(ax= ax, edgecolor= textcolor,facecolor = compass['color'],lw = 0.6,**kwargs)300    ax.text(alon,alat+deltaLon,'N',color = textcolor,fontsize = textsize,ha = 'center',va = 'bottom')301 302