In [42]:
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
% matplotlib inline
In [2]:
def make_basemap(ax, resolution='i'):
bmap = Basemap(projection='lcc',
lat_0=46, lon_0=-62,
resolution=resolution, area_thresh=0.1,
llcrnrlon=-64.0, llcrnrlat=44.8,
urcrnrlon=-59.0, urcrnrlat=46.6,
ax=ax)
return bmap
In [60]:
def draw_basemap(bmap):
bmap.drawcoastlines(color='#9caf9c')
bmap.drawcountries()
bmap.fillcontinents(color='#d8e3d8')
bmap.drawmapboundary(color = 'gray')
bmap.drawmeridians(np.arange(0, 360, 1), color='gray')
bmap.drawparallels(np.arange(-90, 90, 1), color='gray')
return bmap
In [61]:
fig = plt.figure(figsize = (10,8))
ax = fig.add_subplot(111)
m = make_basemap(ax, 'i')
draw_basemap(m)
plt.show()
In [62]:
import fiona
from shapely.geometry import shape, Point, LineString
from shapely.ops import transform
from functools import partial
import pyproj as pp
import os
env = %env
In [63]:
fname = os.path.join(env['HOME'], 'Dropbox/geotransect_test/data/transects/Transects.shp')
In [64]:
ll_nad83 = pp.Proj("+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs")
utm_nad83 = pp.Proj("+init=EPSG:26920")
project = partial(pp.transform, utm_nad83, ll_nad83)
In [65]:
with fiona.open(fname) as f:
print "bounds:", f.bounds
print f.crs['init']
l3 = shape(f[3]['geometry'])
print l3
In [66]:
def plot_element(m, element, colour='b'):
# Plots a line or point given an
# element with lon,lat coordinates.
# Requires a Basemap object.
lo, la = element.xy
x, y = m(lo, la)
return m.plot(x, y, color=colour, linewidth=3, solid_capstyle='round')
Let's reorganize things a bit.
In [67]:
fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(111)
pad = 20000 # Interior padding in m
with fiona.open(fname) as c:
# Get & convert bounds of shapefile.
llx, lly, urx, ury = c.bounds
ll = transform(project, Point(llx-pad, lly-pad))
ur = transform(project, Point(urx+pad, ury+pad))
mid = Point(ll.x + 0.5*(ur.x-ll.x), ll.y + 0.5*(ur.y - ll.y))
# Generate Basemap object.
m = Basemap(projection='tmerc',
lon_0=mid.x, lat_0=mid.y,
resolution='i', # c, l, i, h, f
llcrnrlon=ll.x, llcrnrlat=ll.y,
urcrnrlon=ur.x, urcrnrlat=ur.y)
# Step over shapes in collection.
for s in c:
line = shape(s['geometry'])
line_t = transform(project, line)
plot_element(m, line_t)
# Finish drawing the basemap.
draw_basemap(m)
plt.show()
In [81]:
def draw_basemap(bmap):
bmap.drawcoastlines(color='#9caf9c')
bmap.drawcountries()
# bmap.fillcontinents(color='#d8e3d8')
bmap.drawmapboundary(color = 'gray')
bmap.drawmeridians(np.arange(0, 360, 1), color='gray')
bmap.drawparallels(np.arange(-90, 90, 1), color='gray')
return bmap
In [50]:
import rasterio
from affine import Affine
In [68]:
shade = os.path.join(env['HOME'], 'Dropbox/geotransect_test/data/elevation/DEM_slope.tif')
with rasterio.drivers(CPL_DEBUG=True):
with rasterio.open(shade) as src:
v = src.read()
meta = src.meta
aff = meta['affine']
h = meta['height']
w = meta['width']
meta
Out[68]:
In [52]:
aff
Out[52]:
In [53]:
print aff * (0, 0) # Upper left corner
print aff * (0, h) # Lower left
print aff * (w, 0) # Upper right
# Raster is UTM
# ll = transform(project, Point(aff * (0, h)))
# ur = transform(project, Point(aff * (w, 0)))
ll = Point(aff * (0, 0))
ur = Point(aff * (w, 0))
In [54]:
extent = [ll.coords.xy[0][0], ur.coords.xy[0][0], ll.coords.xy[1][0], ur.coords.xy[1][0]] # L, R, B, T
In [69]:
image = np.repeat(np.swapaxes(np.swapaxes(v, 1, 2), 0, 2), 3, axis=-1)
image.shape
Out[69]:
This is big. It takes ages to plot just an imshow, and is just black...
In [76]:
dem = os.path.join(env['HOME'], 'Dropbox/geotransect_test/data/elevation/NS_DEM_ascii.asc')
It's a monochrome raster so we'll need rasterio.
In [77]:
with rasterio.drivers(CPL_DEBUG=True):
with rasterio.open(shade) as src:
v = src.read()
meta = src.meta
aff = meta['affine']
h = meta['height']
w = meta['width']
In [78]:
image = np.repeat(np.swapaxes(np.swapaxes(v, 1, 2), 0, 2), 3, axis=-1)
image.shape
Out[78]:
In [79]:
image = np.flipud(image)
In [82]:
fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(111)
pad = 20000 # Interior padding in m
with fiona.open(fname) as c:
# Get & convert bounds of shapefile.
llx, lly, urx, ury = c.bounds
ll = transform(project, Point(llx-pad, lly-pad))
ur = transform(project, Point(urx+pad, ury+pad))
mid = Point(ll.x + 0.5*(ur.x-ll.x), ll.y + 0.5*(ur.y - ll.y))
# Generate Basemap object.
m = Basemap(projection='tmerc',
lon_0=mid.x, lat_0=mid.y,
resolution='i', # c, l, i, h, f
llcrnrlon=ll.x, llcrnrlat=ll.y,
urcrnrlon=ur.x, urcrnrlat=ur.y)
# # Step over shapes in collection.
# for s in c:
# line = shape(s['geometry'])
# line_t = transform(project, line)
# plot_element(m, line_t)
m.imshow(image)
# Finish drawing the basemap.
draw_basemap(m)
plt.show()
In [83]:
from osgeo import gdal
In [88]:
fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(111)
pad = 20000 # Interior padding in m
with fiona.open(fname) as c:
# Get & convert bounds of shapefile.
llx, lly, urx, ury = c.bounds
ll = transform(project, Point(llx-pad, lly-pad))
ur = transform(project, Point(urx+pad, ury+pad))
mid = Point(ll.x + 0.5*(ur.x-ll.x), ll.y + 0.5*(ur.y - ll.y))
# Generate Basemap object.
m = Basemap(projection='tmerc',
lon_0=mid.x, lat_0=mid.y,
resolution='i', # c, l, i, h, f
llcrnrlon=ll.x, llcrnrlat=ll.y,
urcrnrlon=ur.x, urcrnrlat=ur.y)
ds = gdal.Open(dem)
data = ds.ReadAsArray()
np.flipud(data)
x = np.linspace(0, m.urcrnrx, data.shape[1])
y = np.linspace(0, m.urcrnry, data.shape[0])
xx, yy = np.meshgrid(x, y)
m.pcolormesh(xx, yy, data)
# Finish drawing the basemap.
draw_basemap(m)
plt.show()
And we should make a 'highlight' map.
In [15]:
with fiona.open(fname) as f:
print "bounds:", f.bounds
print f.crs['init']
print "properties:"
for i in f:
print 'Transect '+str(i['properties']['id'])+': ', i['properties']['basin']
l = list(f)
l[0]
Out[15]:
In [16]:
transect = 3
In [17]:
fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(111)
pad = 20000 # Interior padding in m
res = 'i'
with fiona.open(fname) as c:
bounds = c.bounds
lst = list(c)
if transect:
pad = 5000
res = 'h'
for s in lst:
if s['properties']['id'] == transect:
bounds = shape(s['geometry']).bounds
print bounds
llx, lly, urx, ury = bounds
w, h = urx - llx, ury - lly
x_adj, y_adj = 0, 0
if h > (w*3/8):
x_adj = ((8*h/3.) - w) / 2. # Aspect is hard-coded in uberplot
else: # TODO Fix that!
y_adj = ((3*w/8.) - h) / 2.
print x_adj, y_adj
ll = transform(project, Point(llx-pad-x_adj, lly-pad-y_adj))
ur = transform(project, Point(urx+pad+x_adj, ury+pad+y_adj))
mid = Point(ll.x + 0.5*(ur.x-ll.x), ll.y + 0.5*(ur.y - ll.y))
# Generate Basemap object.
m = Basemap(projection='tmerc',
lon_0=mid.x, lat_0=mid.y,
resolution=res, # c, l, i, h, f
llcrnrlon=ll.x, llcrnrlat=ll.y,
urcrnrlon=ur.x, urcrnrlat=ur.y)
# Step over shapes in collection.
for s in lst:
line = shape(s['geometry'])
line_t = transform(project, line)
if s['properties']['id'] == transect:
col = 'r'
else:
col = 'b'
plot_element(m, line_t, colour=col)
# Finish drawing the basemap.
draw_basemap(m)
plt.show()
In [18]:
transect = 5
In [19]:
fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(111)
pad = 20000 # Interior padding in m
res = 'i'
with fiona.open(fname) as c:
bounds = c.bounds
lst = list(c)
if transect:
pad = 5000
res = 'h'
for s in lst:
if s['properties']['id'] == transect:
bounds = shape(s['geometry']).bounds
print bounds
llx, lly, urx, ury = bounds
w, h = urx - llx, ury - lly
print w, h
x_adj, y_adj = 0, 0
if h > (w*3/8):
print "adjusting x"
x_adj = ((8*h/3.) - w) / 2. # Aspect is hard-coded in uberplot
else: # TODO Fix that!
print "adjusting y"
y_adj = ((3*w/8.) - h) / 2.
print x_adj, y_adj
ll = transform(project, Point(llx-pad-x_adj, lly-pad-y_adj))
ur = transform(project, Point(urx+pad+x_adj, ury+pad+y_adj))
mid = Point(ll.x + 0.5*(ur.x-ll.x), ll.y + 0.5*(ur.y - ll.y))
# Generate Basemap object.
m = Basemap(projection='tmerc',
lon_0=mid.x, lat_0=mid.y,
resolution=res, # c, l, i, h, f
llcrnrlon=ll.x, llcrnrlat=ll.y,
urcrnrlon=ur.x, urcrnrlat=ur.y)
# Step over shapes in collection.
for s in lst:
line = shape(s['geometry'])
line_t = transform(project, line)
if s['properties']['id'] == transect:
col = 'r'
else:
col = 'b'
plot_element(m, line_t, colour=col)
# Finish drawing the basemap.
draw_basemap(m)
plt.show()
In [17]:
with fiona.open(fname) as c:
bounds = c.bounds
lst = [shape(l['geometry']) for l in c]
lst[2]
Out[17]:
In [19]:
basins = os.path.join(env['HOME'], 'Dropbox/geotransect_test/data/basins/basins.shp')
with fiona.open(basins) as f:
print "bounds:", f.bounds
print f.crs['init']
print "properties:"
for i in f.items():
print i[1]['properties']
We will try using the QGIS Composer's interface.
In [6]:
from qgis.core import *
import os
In [7]:
QgsApplication.setPrefixPath("/Applications/QGIS.app/Contents/MacOS")
QgsApplication.initQgis()
In [8]:
env = %env
In [9]:
print "Loading",
data_source = os.path.join(env['HOME'], 'Dropbox/geotransect_test/data/seismic/matt/seismic_lines.shp')
print data_source, os.path.isfile(data_source)
layer_name = "seismic"
provider_name = "ogr"
layer = QgsVectorLayer(data_source, layer_name, provider_name)
if layer.isValid():
print "Success!"
else:
print "Layer failed to load!"
If you see Layer failed to load! then there's not much point in carrying on.
In [10]:
layer.id()
Out[10]:
In [13]:
from PyQt4.QtGui import QImage, QPainter
from PyQt4.QtCore import QSize
In [16]:
# Create image
img = QImage(QSize(800, 600), QImage.Format_ARGB32_Premultiplied)
In [17]:
# Create painter
p = QPainter(img)
#p.begin(img)
p.setRenderHint(QPainter.Antialiasing)
render = QgsMapRenderer()
In [18]:
lst = [] # Layer set
In [20]:
QgsMapLayerRegistry.instance().addMapLayer(layer)
lst.append(layer.id())
render.setLayerSet(lst)
In [ ]:
comp = QgsComposition(render)
comp.setPlotStyle(QgsComposition.Print)
composerMap = QgsComposerMap(comp, 5,5,200,200)
comp.addItem(composerMap)
printer = QPrinter()
printer.setOutputFormat(QPrinter.PdfFormat)
printer.setOutputFileName("out.pdf")
printer.setPaperSize(QSizeF(comp.paperWidth(), comp.paperHeight()), QPrinter.Millimeter)
printer.setFullPage(True)
printer.setColorMode(QPrinter.Color)
printer.setResolution(comp.printResolution())
pdfPainter = QPainter(printer)
paperRectMM = printer.pageRect(QPrinter.Millimeter)
paperRectPixel = printer.pageRect(QPrinter.DevicePixel)
comp.render(pdfPainter, paperRectPixel, paperRectMM)
pdfPainter.end()
app.exitQgis()
In [21]:
# set extent
rect = QgsRectangle(render.fullExtent())
rect.scale(1.1)
render.setExtent(rect)
# set output size
render.setOutputSize(img.size(), img.logicalDpiX())
# do the rendering
render.render(p)
p.end()
# save image
img.save("render.png","png")
Out[21]:
In [66]:
ls -l render.png
The following example is from GIS StackExchange. It doesn't work either. In fact, it often crashes Python.
In [1]:
# Uncomment these to run
# from qgis.core import *
# from qgis.gui import *
# from PyQt4.QtCore import *
# from PyQt4.QtGui import *
# from PyQt4.QtXml import *
QgsApplication.setPrefixPath("/usr", True)
QgsApplication.initQgis()
app = QgsApplication([], True)
fh = open("eg.csv","w")
fh.write("""
x,y,name
153.0278, -27.4679, Brisbane
144.2500, -23.4500, Longreach
145.7753, -16.9256, Cairns
""")
fh.close()
uri = "eg.csv?delimiter=%s&xField=%s&yField=%s" % (",", "x", "y")
layer = QgsVectorLayer(uri, "eglayer", "delimitedtext")
print "layer isValid", layer.isValid()
layerset = []
QgsMapLayerRegistry.instance().addMapLayer(layer)
layerset.append(layer.id())
myMapRenderer = QgsMapRenderer()
myMapRenderer.setLayerSet(layerset)
mapRectangle = QgsRectangle(140,-28,155,-15)
myMapRenderer.setExtent(mapRectangle)
comp = QgsComposition(myMapRenderer)
comp.setPlotStyle(QgsComposition.Print)
composerMap = QgsComposerMap(comp, 5,5,200,200)
comp.addItem(composerMap)
printer = QPrinter()
printer.setOutputFormat(QPrinter.PdfFormat)
printer.setOutputFileName("out.pdf")
printer.setPaperSize(QSizeF(comp.paperWidth(), comp.paperHeight()), QPrinter.Millimeter)
printer.setFullPage(True)
printer.setColorMode(QPrinter.Color)
printer.setResolution(comp.printResolution())
pdfPainter = QPainter(printer)
paperRectMM = printer.pageRect(QPrinter.Millimeter)
paperRectPixel = printer.pageRect(QPrinter.DevicePixel)
comp.render(pdfPainter, paperRectPixel, paperRectMM)
pdfPainter.end()
app.exitQgis()
In [1]:
In [ ]: