#!/usr/bin/python # ./generate-worldmap.py -i mapping/f7.txt -o public_html/maps/testing/7/f7.png # import sys import getopt import GeoIP import random import matplotlib matplotlib.use('Agg2') from pylab import * from matplotlib.numerix import ma from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas from matplotlib.toolkits.basemap import Basemap from matplotlib.colors import LogNorm from matplotlib.ticker import LogFormatterMathtext as JefFormatter gi = GeoIP.open("/usr/share/GeoIP/GeoLiteCity.dat", GeoIP.GEOIP_MEMORY_CACHE) random.seed() def draw_client_density(r={}): if not r.has_key('name'): print "Region Error Aborting" sys.exit(2) m = Basemap(llcrnrlon=r['ll_lon'],llcrnrlat=r['ll_lat'],urcrnrlon=r['ur_lon'],urcrnrlat=r['ur_lat'],\ resolution='l',projection='cyl') # plot them as filled circles on the map. # first, create a figure. dpi=100 dimx=800/dpi dimy=400/dpi fig=figure(figsize=(dimx,dimy), dpi=dpi, frameon=True, facecolor='blue',edgecolor='white') ax=fig.add_axes([0.05,0.1,0.8,0.8],axisbg=(0.05,0.65,0.05)) canvas = FigureCanvas(fig) X,Y,Z,countries = parse_clients(m) # draw coasts and fill continents. m.drawcoastlines(linewidth=0.5) m.drawcountries(linewidth=0.5) m.drawlsmask([100,100,100,0],[100,210,210,255]) palette = cm.YlOrRd # palette = cm.cool m.imshow(Z,palette,extent=(m.xmin,m.xmax,m.ymin,m.ymax), norm=LogNorm(),interpolation='gaussian') l,b,w,h = ax.get_position() cax = axes([l+w+0.005, b, 0.03, h]) colorbar(cax=cax,format=JefFormatter()) # draw colorbar figtext(l+w+0.1,0.5,"Client Density: clients per km^2",va="center",ha="center",rotation="vertical",fontsize=10) figtext(0.05,0.05,footer,backgroundcolor='white',fontsize="smaller",va="bottom") figtext(0.05,0.95,header,backgroundcolor='white',fontsize="smaller",va="top") canvas.print_figure(outputfile, dpi=100,facecolor='white',edgecolor='white') if stats : for country,count in countries.items(): print country," :: ",count def parse_clients(m,latpixels=180,lonpixels=360,lat_smooth=1,lon_smooth=1): rad_deg=pi/180.0 r0=6378.1 latscale=(m.ymax-m.ymin)/latpixels lonscale=(m.xmax-m.xmin)/lonpixels lat_array=arange(m.ymin,m.ymax+latscale,latscale) lon_array=arange(m.xmin,m.xmax+lonscale,lonscale) maxlat=len(lat_array)-1 maxlon=len(lon_array)-1 Z=zeros((len(lat_array),len(lon_array)),dtype='float') seen_countries={} f = open(inputfile, 'r') for line in f: gir=None try: gir = gi.record_by_addr(line.strip()) except: continue if gir != None: country_code=str(gir['country_code']) client = (line.strip(), gir['country_code'], gir['latitude'], gir['longitude']) # Deal with client if stats == True: if seen_countries.has_key(country_code) : seen_countries[country_code]+=1 else : seen_countries[country_code]=1 lat=client[2] lon=client[3] if ( lat == 0.0 ) and ( lon == 0.0 ) : # print "Lat/Lon 0.0:",country_code lat = -89.0 i_lat=int(float((lat-lat_array[0]))/float(latscale)) i_lon=int(float((lon-lon_array[0]))/float(lonscale)) for i in xrange(-int(lat_smooth),int(lat_smooth+1),1): for j in xrange(-int(lon_smooth),int(lon_smooth+1),1): if ( i_lat+i >= 0 ) and (i_lat+i < maxlat) : if ( i_lon+j >= 0) and ( i_lon+j < maxlon) : Z[i_lat+i,i_lon+j]+=1.0 f.close() # End of file, now do summary density calculation for i in xrange(len(lat_array)): area=r0*r0*rad_deg*lonscale*abs( sin(rad_deg*(lat_array[i]-latscale/2.0) )-sin(rad_deg*(lat_array[i]+latscale/2.0) ) ) for j in xrange(len(lon_array)): if area == 0.0 : Z[i,j]=0.0 else: Z[i,j]=Z[i,j]/area/(2.0*lon_smooth+1)/(2.0*lat_smooth+1) Lon,Lat=meshgrid(lon_array,lat_array) X,Y=m(Lon,Lat) Z = where(Z <= 0.,1.e10,Z) Z = ma.masked_values(Z, 1.e10) return X,Y,Z,seen_countries def parse_regions(region=None): r={} if region is None: r={'name':'World','ll_lat':-90,'ur_lat':90,'ll_lon':-180,'ur_lon':180} else: if not (regionsfile is None): try: f = open(regionsfile, 'r') for line in f: if region==line.strip().split(":")[0].strip(): r['name']=line.strip().split(":")[0].strip() r['ll_lat']=float(line.strip().split(":")[1]) r['ur_lat']=float(line.strip().split(":")[2]) r['ll_lon']=float(line.strip().split(":")[3]) r['ur_lon']=float(line.strip().split(":")[4]) break f.close() except: print "Error parsing regionsfile %s" % regionsfile sys.exit(2) return r def list_regions(): if not (regionsfile is None): try: f = open(regionsfile, 'r') for line in f: print line.strip().split(":")[0].strip() f.close() except: print "Error parsing regionsfile %s" % regionsfile sys.exit(2) def usage(): print "%s --header=header --footer=footer --input=inputfile --output=outputfile --indir=indir" % sys.argv[0] +\ " --region=region --regionsfile=regionsfile --outdir=outdir -v" def main(): try: opts, args = getopt.getopt(sys.argv[1:], "h:f:i:o:v:r:s",\ ["help","verbose","stats","list-regions","indir=","outdir=","input=","output=","footer=","header=","region=","regionsfile="]) except getopt.GetoptError: # print help information and exit: usage() sys.exit(2) global header, footer, inputfile,outputfile,indir,outdir, stats, verbose,regionsfile verbose = False stats = False listregions=False inputfile=None outputfile=None regionsfile=None region=None header="" footer="" indir='./' outdir='./' for o, a in opts: if o in ("-h","--header"): header = a if o in ("-f","--footer"): footer = a if o in ("-i","--input"): inputfile = a if o in ("-o","--output"): outputfile = a if o == "--indir": indir = a if o == "--outdir": outdir = a if o == "--regionsfile": regionsfile = a if o in ("-r","--region"): region = a if o in ("-v","--verbose"): verbose = True if o in ("-s","--stats"): stats = True if o in ("--help"): usage() sys.exit() if o in ("--list-regions"): listregions=True if outputfile is None : outputfile = outdir+'/clientmap.png' if inputfile is None : inputfile = indir+'/ips.txt' if regionsfile is None : regionsfile = indir+'/regions.txt' if listregions: list_regions() sys.exit() r=parse_regions(region) draw_client_density(r) if __name__ == "__main__": sys.exit(main())