Basemap je dio Matplotlib biblioteke s funkcijama i metodama za pravljenje geografskih karata. Korištenje Basemapa ilustrirat ćemo na nekim primjerima.
In [19]:
# osnovne biblioteke koje ćemo koristiti
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
%matplotlib inline
import numpy as np
In [2]:
fig = plt.figure(figsize=(14,10))
ax = fig.add_axes([0.1,0.1,0.9,0.9])
# Mercatorova projekcija
m = Basemap(projection='merc',llcrnrlat=-80,urcrnrlat=80,\
llcrnrlon=-180,urcrnrlon=180,lat_ts=20,resolution='c')
# geografske koordinate Zagreba
zaglat = 45.48; zaglon = 15.58
# geografske koordinate Pekinga
pelat = 39.55; pelon = 116.25
m.drawgreatcircle(zaglon,zaglat,pelon,pelat,linewidth=2,color='b')
m.drawcoastlines()
m.fillcontinents(color='coral',lake_color='aqua')
m.drawmapboundary(fill_color='aqua')
m.drawparallels(np.arange(-90.,91.,30.),labels=[1,1,0,1])
m.drawmeridians(np.arange(-180.,181.,60.),labels=[1,1,0,1])
zx,zy=m(zaglon,zaglat)
px,py=m(pelon,pelat)
ax.text(zx,zy,'Zagreb',fontsize=14,fontweight='bold',ha='center',va='center',color='b')
ax.text(px,py,'Peking',fontsize=14,fontweight='bold',ha='center',va='center',color='r')
ax.set_title('Zagreb - Peking');
In [3]:
from datetime import datetime
# Millerova projekcija
fig = plt.figure(figsize=(8,8))
ax = fig.add_axes([0,0,1,1])
map = Basemap(projection='mill',lon_0=90)
map.drawcoastlines()
map.drawparallels(np.arange(-90,90,30),labels=[1,0,0,0])
map.drawmeridians(np.arange(map.lonmin,map.lonmax+30,60),labels=[0,0,0,1])
map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='coral',lake_color='aqua')
date = datetime.utcnow()
CS=map.nightshade(date)
ax.set_title(u'Mapa dana/noći za %s (UTC)' % date.strftime("%d %b %Y %H:%M:%S"));
In [4]:
# Standardne biblioteke
from dateutil.tz import tzutc, tzlocal
# Ostalo što nam treba za crtanje
from matplotlib.colors import Normalize
from matplotlib.cm import get_cmap
Podaci se dohvaćaju sa stranice earthquake.usgs.gov. U datoteku query.csv spremamo podatke o zemljotresima u zadnjih tjedan dana.
In [85]:
url = "https://earthquake.usgs.gov/fdsnws/event/1/query.csv"
import time
import datetime as dt
today = dt.datetime.now()
now = today - dt.timedelta(hours=1)
endtime= now.isoformat()
starttime = (now - dt.timedelta(days=7)).isoformat()
params = dict(minmagnitude='2.5',orderby='time',starttime=starttime,endtime=endtime)
In [86]:
import requests
r = requests.get(url, params = params)
with open('query.csv','w') as file:
file.write(r.text)
In [87]:
# Parsiranje
import csv
lat, lon, mag = [], [], []
with open('query.csv', 'r') as csvfile:
reader = csv.DictReader(csvfile)
for redak in reader:
lon.append(float(redak['longitude']))
lat.append(float(redak['latitude']))
mag.append(float(redak['mag']))
In [88]:
# Normalizacija magnituda
norm = Normalize()
mag_norms = norm(np.array(mag))
z = (mag_norms * 10.0)**2.0
In [89]:
fig = plt.figure(figsize=(14,10))
m = Basemap(projection='moll', lat_0=0, lon_0=0, resolution='l', area_thresh=1000.0)
m.drawcoastlines()
m.fillcontinents(color='coral',lake_color='aqua')
m.drawmapboundary(fill_color='aqua')
m.drawmeridians(np.arange(0, 360, 30))
m.drawparallels(np.arange(-90, 90, 30))
cmap = get_cmap('jet')
x,y = m(lon, lat)
sc = m.scatter(x, y, s=z, cmap=cmap, c=mag)
plt.colorbar(sc, orientation='horizontal')
plt.title("Potresi");
Podaci se preuzimaju sa stranice CISL Research Data Archive (http://rda.ucar.edu/). Podaci su u netCDF formatu, koji je u verziji 4 u biti podskup HDF formata (Hierarchical Data Format). Opis formata možete naći ovdje.
Opis konkretnih podataka je ovdje, a sami podaci se nalaze ovdje.
In [90]:
from netCDF4 import Dataset
#f = Dataset('avhrr-only-v2.20150324.nc')
f = Dataset('avhrr-only-v2.20170330.nc')
In [91]:
# Kako su strukturirani podaci
print (f)
In [92]:
print(f.variables.keys()) # imena varijabli
temp = f.variables['sst'] # temperatura
print(temp)
In [93]:
temp.dimensions
Out[93]:
In [94]:
data = temp[0,0]
data.shape
Out[94]:
In [95]:
lonvals = f.variables['lon'][:]
latvals = f.variables['lat'][:]
In [96]:
plt.rcParams['figure.figsize'] = (12.0, 8.0)
plt.contourf(lonvals, latvals, data, 20, cmap=plt.get_cmap('YlGnBu_r'))
plt.colorbar()
plt.show()
Sad ćemo nacrati iste podatke na zemaljskoj kugli.
In [97]:
m = Basemap(projection='ortho', lon_0=-50, lat_0=40, resolution='l')
In [98]:
X,Y = np.meshgrid(lonvals, latvals)
x,y = m(X,Y)
In [99]:
pc = m.contourf(x, y, data, 30, cmap=plt.get_cmap('YlGnBu_r'))
m.bluemarble()
m.drawmapboundary()
m.drawcoastlines()
plt.title('Ortografska projekcija')
plt.colorbar(pc, orientation='vertical')
plt.show()
In [18]:
from verzije import *
from IPython.display import HTML
HTML(print_sysinfo()+info_packages('matplotlib,basemap'))
Out[18]: