notebook.community

Edit and run



In [100]:

    
#!/usr/bin/env python
# -*- coding: utf-8 -*- 
from IPython.core.display import display, HTML #HTML output for ipython notebook
from osgeo import gdal
from geopy.geocoders import Nominatim
from collections import OrderedDict
from prettytable import PrettyTable
from prettytable import from_html
import calendar
import os
import numpy as np
from numpy import ma
from collections import Counter



In [101]:

    
path = 'data'
city = 'Partizansk'
bit_value = 4 # bit number for values
g = Nominatim().geocode(city, timeout=5)
position = (g.longitude, g.latitude)



In [102]:

    
def get_value_at_point(rasterfile, pos):
    # Extract raster data value at a point:
    gdata = gdal.Open(rasterfile)
    gt = gdata.GetGeoTransform()
    band = gdata.GetRasterBand(1)
    nodata = band.GetNoDataValue()     
    data = gdata.ReadAsArray().astype(np.float)
    gdata = None
    masked_data = ma.masked_values(data, nodata, copy=False) # mask no value data
    masked_data.fill_value = nodata
    d_values = np.around(masked_data, decimals=bit_value) # rounding values
    x = int((pos[0] - gt[0])/gt[1])
    y = int((pos[1] - gt[3])/gt[5])    
    x_value = round(d_values[y, x], bit_value) # desired quantity
    
    # search coordinates which corresponded by overlapped value:
    search_pix = zip(*np.where(d_values == x_value)) # search pixel cootdinates
    overlapped = map(lambda pix :(round(gt[3] + pix[1] * gt[4] + pix[0] * gt[5], 4), # lat
                  round(gt[0] + pix[1] * gt[1] + pix[0] * gt[2], 4)), search_pix) # # long
    
    return x_value, overlapped



In [103]:

    
def ext_data():
    result_data = []
    val_dataset = sorted(os.listdir(path))
    for val in val_dataset: # list of datasets by climate values
        res = []
        ovp = []
        for m in [x for x in range(1,13)]: #get data by month(from datasets like wc2.0_5m_tmax_01.tif) 
            result = get_value_at_point('data/%s/wc2.0_5m_%s_%02d.tif' % (val, val, m), position)
            res.append(result) 

        result_data.append((val, res))
    return result_data, val_dataset



In [104]:

    
def output_table(extdata):
    mgen_names = [calendar.month_name[x] for x in range(1,13)]
    table = PrettyTable()
    table.add_column("Month", mgen_names)
    for k in extdata[0]:
        vals = [j[0] for j in k[1]]
       # print vals
        table.add_column(k[0], vals)
    print(city)
    #print(table)
    #HTML output for ipython notebook:
    htmltable = table.get_html_string()
    display(HTML(htmltable))



In [105]:

    
extdata = ext_data() # extract data



In [106]:

    
output = output_table(extdata) # print table



In [107]:

    
def check_func(rasterfile, pos):
    gdata = gdal.Open(rasterfile)
    gt = gdata.GetGeoTransform()
    band = gdata.GetRasterBand(1)
    nodata = band.GetNoDataValue()     
    data = gdata.ReadAsArray().astype(np.float)
    gdata = None
    masked_data = ma.masked_values(data, nodata, copy=False) # mask no value data
    masked_data.fill_value = nodata
    d_values = np.around(masked_data, decimals=bit_value) # rounding values
    x = int((pos[0] - gt[0])/gt[1])
    y = int((pos[1] - gt[3])/gt[5]) 
    x_value = round(d_values[y, x], bit_value) 
    
    return rasterfile, (pos[1],pos[0]), x_value



In [108]:

    
def check_ext_data(extdata):
    result_data = []
    for i in extdata[0]:
        for idx, k in enumerate(i[1]): # i[1] - months
            print calendar.month_name[idx+1]
            for j in k[1]:
                #print j, k[0]
                result_data.append(j)
                result = check_func('data/tavg/wc2.0_5m_tavg_%02d.tif' % (idx+1), (j[1], j[0]) ) 
                print result
    return result_data



In [109]:

    
#check_ext_data(extdata)



In [110]:

    
import gmaps
import gmaps.datasets
gmaps.configure(api_key=" ")



In [111]:

    
def mapdata(extdata): # get array of overlapped coordinates
    result_data = []
    
    for i in extdata[0]:
        for idx, k in enumerate(i[1]): # i[1] - months
            result_data.append(k[1])
    
    #for i in extdata[0]:        
    #    for k in i[1][0][1]:
            
    return result_data



In [112]:

    
mapdata = mapdata(extdata)



In [113]:

    
#print mapdata[0]



In [ ]:

    
color_list = ['Black','Blue','Purple','Aqua','Lime','Teal','Green','Olive','Maroon','Red','Yellow','Gray']
#winter:'Black','Blue','Gray'
#spring: 'Purple','Aqua','Lime'
#summer: 'Teal','Green','Olive'
#fall: 'Maroon','Red','Yellow'

fig = gmaps.figure()
for idx, month_data in enumerate(mapdata):
    print color_list[idx]+": "+ calendar.month_name[idx+1]
    climate_layer = gmaps.symbol_layer(
                    month_data, fill_color=color_list[idx], stroke_color=color_list[idx], scale=3)

    fig.add_layer(climate_layer)

fig









    



Black: January



In [ ]:



In [65]:

    
geolocator = Nominatim()
for place in mapdata[8]:
    location = geolocator.reverse(place, timeout=100, language='en')
    print place, location.address









    



(50.5833, 47.25) Первопитерский, Pitersky District, Saratov Oblast, Volga Federal District, Russian Federation
(49.5, 62.8333) Тармак, Жангельдинский район, Kostanay Region, Kazakhstan
(48.25, 7.5833) D 705, Sundhouse, Sélestat-Erstein, Bas-Rhin, Great East, 67920, France
(47.8333, 16.3333) Güterweg Weichselbühel, Heutalhof, Gemeinde Lichtenwörth, Bezirk Wiener Neustadt, Lower Austria, 2493, Austria
(47.5833, 84.4167) Сарыбулак, Тарбагатайский район, East Kazakhstan Region, Kazakhstan
(47.3333, 27.9167) Ungheni-Chișinău, Hristoforovca, raionul Ungheni, 3641, Moldova
(47.3333, 36.5833) Т-08-03, Куйбишевська селищна рада, Більмацький район, Zaporizhia Oblast, Ukraine
(47.25, 76.1667) Актогайский район, Karaganda Region, Kazakhstan
(47.0, 26.9167) DN2, Traian, Neamț, 617398, Romania
(46.8333, 16.4167) Baksaszer, Őriszentpéter (Baksaszer), Baksaszer, Őriszentpéter, Körmendi járás, Vas, Western Transdanubia, Transdanubia, 9941, Hungary
(46.5833, 16.1667) 3, Prvomajska ulica, Veržej, Mura Statistical Region, 9241, Slovenia
(46.3333, 4.25) La Raie, Amanzé, Charolles, Saône-et-Loire, Bourgogne-Franche-Comté, 71800, France
(46.0, -104.0) Camp Crook Road, Bowman County, North Dakota, United States of America
(45.9167, -85.25) Mackinac County, Michigan, United States of America
(44.6667, -92.25) 370th Avenue, Ono, Pierce County, Wisconsin, 54750, United States of America
(44.4167, -87.6667) Pine Grove Road, Stangelville, Kewaunee County, Wisconsin, 54217, United States of America
(44.1667, -122.75) 16-1-26 Road, Lane County, Oregon, United States of America
(44.0, -85.8333) West 1 Mile Road, Wolf Lake, Lake County, Michigan, 49304, United States of America
(43.9167, -116.1667) Twilegar Lane, Boise County, Idaho, 83629, United States of America
(43.8333, -84.5) Lyle Road, Gladwin County, Michigan, 48612, United States of America
(43.8333, -84.25) Townline 16 Road, Mills Township, Midland County, Michigan, 48652, United States of America
(43.75, 131.8333) Школьная улица, Linevichi, Уссурийский городской округ, Primorsky Krai, Far Eastern Federal District, 692508, Russian Federation
(43.1667, 133.0833) 1-я Красноармейская улица, Партизанский городской округ, Primorsky Krai, Far Eastern Federal District, Russian Federation
(42.25, 98.0) Ejin Banner, Alxa League, Inner Mongolia, China
(42.0, -105.1667) Cooney Hills Road, Platte County, Wyoming, United States of America
(41.75, -74.75) Swan Lake Road, Ferndale, Town of Liberty, Sullivan County, New York, 12734, United States of America
(41.5, -79.1667) Francis Road, Blue Jay, Forest County, Pennsylvania, United States of America
(41.4167, 74.75) Ugut, Naryn Region, Kyrgyzstan
(38.5, 107.9167) Otog Front Banner, Ordos City, Inner Mongolia, China
(36.6667, 48.8333) امام, بخش مرکزی, Zanjan County, Zanjan Province, Iran
(35.25, 64.1667) R17, Jawand, Badghis, Afghanistan
(32.6667, 49.5) لبدعلیا, بخش بازفت, Kuhrang County, Chaharmahal and Bakhtiari Province, Iran
(9.0, 39.0833) Addis Ababa, Ethiopia
(-0.25, 36.1667) C69, Nakuru, Central, Kenya
(-19.1667, -65.0833) RN6: Sucre-Tarabuco, Municipio Yamparaez, Provincia Yamparaez, CHQ, Bolivia
(-27.0, -51.4167) Rua Caçador, Treze Tílias, Microrregião de Joaçaba, Western Santa Catarina Mesoregion, Santa Catarina, South Region, Brazil
(-27.1667, -49.3333) Apiúna, Microrregião de Blumenau, Mesorregião do Vale do Itajaí, Santa Catarina, South Region, Brazil
(-27.75, 15.75) Restricted, ǁKaras Region, Namibia
(-28.8333, 26.0833) Mangaung Ward 44, Bloemfontein, Mangaung Metropolitan Municipality, Free State, RSA
(-29.6667, -64.0) Ruta Provincial 18, Sobremonte, Córdoba, Argentina
(-30.5, 122.8333) Western Australia, Australia
(-31.3333, -55.5) Rivera, 40000, Uruguay
(-31.3333, 131.0833) Nullarbor, South Australia, 5690, Australia
(-31.75, -68.25) Divisoria, Caucete, SJ, Argentina
(-32.0, 136.0) Moonaree, South Australia, 5717, Australia
(-32.0833, 136.75) South Australia, Australia
(-32.8333, 142.5) Menindee, Central Darling Shire Council, New South Wales, 2879, Australia
(-33.0833, 142.1667) Pooncarie, Wentworth Shire Council, New South Wales, 2648, Australia



In [ ]:

    
#from osgeo import gdal
#import matplotlib.pyplot as plt
#import numpy as np
#from mpl_toolkits.basemap import Basemap
#
## Plotting 2070 projected August (8) precip from worldclim
#gdata = gdal.Open("data/tavg/wc2.0_10m_tavg_01.tif")
#geo = gdata.GetGeoTransform()
#data = gdata.ReadAsArray()
#
#xres = geo[1]
#yres = geo[5]
#
## A good LCC projection for USA plots
#m = Basemap(llcrnrlon=-119,llcrnrlat=22,urcrnrlon=-64,urcrnrlat=49,
#            projection='lcc',lat_1=33,lon_0=-95)
#
#xmin = geo[0] + xres * 0.5
#xmax = geo[0] + (xres * gdata.RasterXSize) - xres * 0.5
#ymin = geo[3] + (yres * gdata.RasterYSize) + yres * 0.5
#ymax = geo[3] - yres * 0.5
#
#x,y = np.mgrid[xmin:xmax+xres:xres, ymax+yres:ymin:yres]
#x,y = m(x,y)
#
#cmap = plt.cm.gist_rainbow
#cmap.set_under ('1.0')
#cmap.set_bad('0.8')
#
#im = m.pcolormesh(x,y, data.T, cmap=cmap, vmin=45, vmax=0)
#
#cb = plt.colorbar( orientation='vertical', fraction=0.10, shrink=0.7)
#plt.title('Tavg')
#plt.show()



In [ ]:

Month	prec
January	14.0
February	15.0
March	28.0
April	49.0
May	65.0
June	84.0
July	106.0
August	131.0
September	113.0
October	56.0
November	39.0
December	22.0