qgis_script_mario.ipy

# set opacity of blayer 0 and redraw
def setopacity(layer,opacity):
     layer.renderer().setOpacity(opacity)
     layer.triggerRepaint()

# Load some basic configuration
# %load https://gist.githubusercontent.com/kidpixo/2ec078d09834b5aa7869/raw/c8812811211dc7cd62f5b530b51b0104f39263ff/ipython%20inizialization
#import pandas as pd
#import numpy as np
%matplotlib qt4 
import qgis
import matplotlib
from matplotlib import pyplot as plt

#######################################################################################################################
# Load files 
path = '/home/esacgis2015/Desktop/qgis_data/hrsc'

file = !ls -1 $path/H0360*.IMG
for f in file: qgis.utils.iface.addRasterLayer(f)
qgis.utils.iface.addRasterLayer("/home/esacgis2015/Desktop/DATA/global-data/Mars_MGS_MOLA_DEM_mosaic_global_463m-lsscale4-reduced.cub","global_mola")
# see http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/raster.html
canvas = qgis.utils.iface.mapCanvas()

#######################################################################################################################
# print layer props:

# layers name
for i in canvas.layers(): print i.name()

#see http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/raster.html

blayer = canvas.layer(0)
print blayer.width(), blayer.height()
print blayer.extent().toString()
print blayer.rasterType()
# 0 = GrayOrUndefined (single band), 1 = Palette (single band), 2 = Multiban
print blayer.bandCount()

print blayer.metadata()

#######################################################################################################################
# add pesudocolor to MOLA: 

mystat =  blayer.dataProvider().bandStatistics(1, QgsRasterBandStats.Min | QgsRasterBandStats.Max )
# see http://gis.stackexchange.com/a/134589
print mystat.minimumValue,mystat.maximumValue

min, max = mystat.minimumValue,mystat.maximumValue
range =max-min
add=range//2
int =min + add

colDic={'red':'#ff0000', 'yellow':'#ffff00','blue':'#0000ff'}

valueList =[min, int, max]

lst = [ QgsColorRampShader.ColorRampItem(valueList[0], QColor(colDic['red'])), \
    QgsColorRampShader.ColorRampItem(valueList[1], QColor(colDic['yellow'])), \
    QgsColorRampShader.ColorRampItem(valueList[2], QColor(colDic['blue']))]

myRasterShader = QgsRasterShader()
myColorRamp = QgsColorRampShader()

myColorRamp.setColorRampItemList(lst)
myColorRamp.setColorRampType(QgsColorRampShader.INTERPOLATED)
myRasterShader.setRasterShaderFunction(myColorRamp)

myPseudoRenderer = QgsSingleBandPseudoColorRenderer(\
    blayer.dataProvider(), blayer.type(),  myRasterShader)

blayer.setRenderer(myPseudoRenderer)

iface.mapCanvas().refresh()

setopacity(blayer,0.4)

#######################################################################################################################
# Extract data array

import gdal
# from http://www.maths.lancs.ac.uk/~rowlings/Software/Qgis/Plugins/rasterlang/python.html
def layerAsArray(layer):
    """ read the data from a single-band layer into a numpy/Numeric array.
    Only works for gdal layers!
    """
    gd = gdal.Open(str(layer.source()))
    array = gd.ReadAsArray()
    return array

layers = {}
for i in canvas.layers():
   layers[i.name().split('_')[-1][0:2]] = i

#######################################################################################################################
# zooming around
# to full extent
canvas.zoomToFullExtent()

#######################################################################################################################
# to extent of first layer
canvas.setExtent(layers['DT'].extent()) ; iface.mapCanvas().refresh()

#######################################################################################################################
# load data as array:
dt = layerAsArray(layers['DT'])
print dt.dtype,dt.shape

fig3 = plt.figure()
plt.subplot(221)
plt.imshow(dt)
plt.subplot(222)
plt.plot(dt[:,1500])
plt.subplot(223)
hist, bins = np.histogram(dt, bins=128)
width = 0.7 * (bins[1] - bins[0])
center = (bins[:-1] + bins[1:]) / 2
plt.bar(center, hist, align='center', width=width)
plt.subplot(224)
plt.plot(dt[1500,:])

#######################################################################################################################
# to do!!!
dt = dt.astype(float)
dt [dt == -32768.] = np.nan
ir = layerAsArray(layers['IR'])
ir = ir.astype(float)
ir [ir == 0.] = np.nan

# To Big!!!
plt.scatter(dt.reshape(-1),ir[:,1:].reshape(-1),s=10,alpha=0.75)