In [1]:

    

# import the owslib package
import owslib.wms
import numpy as np
import matplotlib.pyplot as plt
from skimage import color

%matplotlib inline


    

In [2]:

    

# get the wms
wmsurl   = 'http://geodata1.nationaalgeoregister.nl/luchtfoto/wms'
wms      = owslib.wms.WebMapService(wmsurl)


    

In [3]:

    

# let's check out the contents of the wms
wms.contents


    

    
Out[3]:





{'luchtfoto': <owslib.wms.ContentMetadata instance at 0x000000000B0AD808>,
 'luchtfoto_jpeg': <owslib.wms.ContentMetadata instance at 0x000000000B0AF3C8>,
 'luchtfoto_png': <owslib.wms.ContentMetadata instance at 0x000000000B0AF608>}

In [4]:

    

# we will use the png layer
wmslayer = ['luchtfoto_png']


    

In [5]:

    

# set the GPS coordinates of Delft
# top, left, bottom, right = (52.02,4.30,51.97,4.40)
top, left, bottom, right = (52.013,4.357,52.011,4.364)
bbox = (left, bottom, right, top)


    

In [6]:

    

# Get the x and y proportionalities right
ypix = 500
xpix = ypix * np.cos(top/360.*2*np.pi) * np.sin(right-left) / np.sin(top-bottom)
size  = (int(round(xpix)), ypix)
print 'Size [pixels]:',size


    

    




Size [pixels]: (1077, 500)

In [7]:

    

# Set coordinate system
coor   = 'EPSG:4326'


    

In [8]:

    

# get map
f = wms.getmap(layers = wmslayer,
               srs    = coor,
               bbox   = bbox,
               size   = size,
               format = 'image/png')


    

In [9]:

    

# plot the map


img = plt.imread(f)

fig, ax = plt.subplots()
ax.imshow(img)
ax.set_title('Aerial photograph of Delft')


    

    
Out[9]:





<matplotlib.text.Text at 0xed5ac88>






    

In [10]:

    

# The map looks good, but a bit gray. Let's increase the saturation

# Let us increase the saturation with 50 percent
# Convert rgb to hsv (hue saturation value), increase saturation and convert
# back
imgsat = color.rgb2hsv(img)
imgsat[:,:,1] *= 1.5
imgsat = color.hsv2rgb(imgsat)


    

In [11]:

    

# Now plot the map again

fig = plt.figure()
ax = fig.add_subplot(111)
ax.imshow(imgsat)
ax.set_title('Aerial photograph of Delft, 150% saturation')
plt.show()


    

In [ ]: