CoCalc Public Files07. Jupyter-ostalo / basemap.ipynb
Author: Ivica Nakić

# Basemap

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


## Crtanje najkraće udaljenosti

In [2]:
fig = plt.figure(figsize=(14,10))
# 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');

/projects/anaconda3/lib/python3.5/site-packages/mpl_toolkits/basemap/__init__.py:3327: MatplotlibDeprecationWarning: The ishold function was deprecated in version 2.0. b = ax.ishold() /projects/anaconda3/lib/python3.5/site-packages/mpl_toolkits/basemap/__init__.py:3336: MatplotlibDeprecationWarning: axes.hold is deprecated. See the API Changes document (http://matplotlib.org/api/api_changes.html) for more details. ax.hold(b) /projects/anaconda3/lib/python3.5/site-packages/mpl_toolkits/basemap/__init__.py:1800: MatplotlibDeprecationWarning: The get_axis_bgcolor function was deprecated in version 2.0. Use get_facecolor instead. axisbgc = ax.get_axis_bgcolor()

## Dan i noć

In [3]:
from datetime import datetime
# Millerova projekcija
fig = plt.figure(figsize=(8,8))
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()
ax.set_title(u'Mapa dana/noći za %s (UTC)' % date.strftime("%d %b %Y %H:%M:%S"));

/projects/anaconda3/lib/python3.5/site-packages/mpl_toolkits/basemap/__init__.py:1800: MatplotlibDeprecationWarning: The get_axis_bgcolor function was deprecated in version 2.0. Use get_facecolor instead. axisbgc = ax.get_axis_bgcolor() /projects/anaconda3/lib/python3.5/site-packages/mpl_toolkits/basemap/__init__.py:3683: MatplotlibDeprecationWarning: The ishold function was deprecated in version 2.0. b = ax.ishold() /projects/anaconda3/lib/python3.5/site-packages/mpl_toolkits/basemap/__init__.py:3719: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future xx = x[x.shape[0]/2,:] /projects/anaconda3/lib/python3.5/site-packages/mpl_toolkits/basemap/__init__.py:3752: MatplotlibDeprecationWarning: axes.hold is deprecated. See the API Changes document (http://matplotlib.org/api/api_changes.html) for more details. ax.hold(b)

# Crtanje potresa

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:
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


### Crtanje

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");

/projects/anaconda3/lib/python3.5/site-packages/mpl_toolkits/basemap/__init__.py:1656: MatplotlibDeprecationWarning: The get_axis_bgcolor function was deprecated in version 2.0. Use get_facecolor instead. fill_color = ax.get_axis_bgcolor() /projects/anaconda3/lib/python3.5/site-packages/mpl_toolkits/basemap/__init__.py:1800: MatplotlibDeprecationWarning: The get_axis_bgcolor function was deprecated in version 2.0. Use get_facecolor instead. axisbgc = ax.get_axis_bgcolor() /projects/anaconda3/lib/python3.5/site-packages/mpl_toolkits/basemap/__init__.py:3289: MatplotlibDeprecationWarning: The ishold function was deprecated in version 2.0. b = ax.ishold() /projects/anaconda3/lib/python3.5/site-packages/mpl_toolkits/basemap/__init__.py:3298: MatplotlibDeprecationWarning: axes.hold is deprecated. See the API Changes document (http://matplotlib.org/api/api_changes.html) for more details. ax.hold(b)

## Crtanje SST (sea surface temperature)

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)

<class 'netCDF4._netCDF4.Dataset'> root group (NETCDF3_CLASSIC data model, file format NETCDF3): Conventions: CF-1.0 title: Daily-OI-V2, Final, Data (Ship, Buoy, AVHRR: NOAA19, METOP, NCEP-ice) History: Version 2.0 creation_date: 2017-04-14 11:10 Description: Reynolds, et al.(2007) Daily High-resolution Blended Analyses. Available at ftp://eclipse.ncdc.noaa.gov/pub/OI-daily/daily-sst.pdf Climatology is based on 1971-2000 OI.v2 SST, Satellite data: Navy NOAA19 METOP AVHRR, Ice data: NCEP ice Source: NOAA/National Climatic Data Center Contact: Dick Reynolds, email: [email protected] & Chunying Liu, email: [email protected] dimensions(sizes): time(1), zlev(1), lat(720), lon(1440) variables(dimensions): float32 time(time), float32 zlev(zlev), float32 lat(lat), float32 lon(lon), int16 sst(time,zlev,lat,lon), int16 anom(time,zlev,lat,lon), int16 err(time,zlev,lat,lon), int16 ice(time,zlev,lat,lon) groups: 
In [92]:
print(f.variables.keys()) # imena varijabli
temp = f.variables['sst']  # temperatura
print(temp)

odict_keys(['time', 'zlev', 'lat', 'lon', 'sst', 'anom', 'err', 'ice']) <class 'netCDF4._netCDF4.Variable'> int16 sst(time, zlev, lat, lon) long_name: Daily sea surface temperature units: degrees C _FillValue: -999 add_offset: 0.0 scale_factor: 0.01 valid_min: -300 valid_max: 4500 unlimited dimensions: current shape = (1, 1, 720, 1440) filling off
In [93]:
temp.dimensions

('time', 'zlev', 'lat', 'lon')
In [94]:
data = temp[0,0]
data.shape

(720, 1440)
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()

/projects/anaconda3/lib/python3.5/site-packages/mpl_toolkits/basemap/__init__.py:3683: MatplotlibDeprecationWarning: The ishold function was deprecated in version 2.0. b = ax.ishold() /projects/anaconda3/lib/python3.5/site-packages/mpl_toolkits/basemap/__init__.py:3752: MatplotlibDeprecationWarning: axes.hold is deprecated. See the API Changes document (http://matplotlib.org/api/api_changes.html) for more details. ax.hold(b) /projects/anaconda3/lib/python3.5/site-packages/mpl_toolkits/basemap/__init__.py:1656: MatplotlibDeprecationWarning: The get_axis_bgcolor function was deprecated in version 2.0. Use get_facecolor instead. fill_color = ax.get_axis_bgcolor() /projects/anaconda3/lib/python3.5/site-packages/mpl_toolkits/basemap/__init__.py:3363: MatplotlibDeprecationWarning: The ishold function was deprecated in version 2.0. b = ax.ishold() /projects/anaconda3/lib/python3.5/site-packages/mpl_toolkits/basemap/__init__.py:3372: MatplotlibDeprecationWarning: axes.hold is deprecated. See the API Changes document (http://matplotlib.org/api/api_changes.html) for more details. ax.hold(b)