dem_derivates.py 7.64 KB
#! /usr/bin/env python
# -*- coding: iso-8859-1 -*-

# Michel Le Page, 06/04/2013

import os
from math import *
import numpy as np
from osgeo import gdal
import csv
from osgeo.gdalconst import *
from osgeo import gdal, gdal_array
import sys

# =============================================================================
# rgb_to_hsv()
#
# rgb comes in as [r,g,b] with values in the range [0,255].  The returned
# hsv values will be with hue and saturation in the range [0,1] and value
# in the range [0,255]
#
def rgb_to_hsv( r,g,b ):

    maxc = np.maximum(r,np.maximum(g,b))
    minc = np.minimum(r,np.minimum(g,b))

    v = maxc

    minc_eq_maxc = np.equal(minc,maxc)

    # compute the difference, but reset zeros to ones to avoid divide by zeros later.
    ones = np.ones((r.shape[0],r.shape[1]))
    maxc_minus_minc = np.choose( minc_eq_maxc, (maxc-minc,ones) )

    s = (maxc-minc) / np.maximum(ones,maxc)
    rc = (maxc-r) / maxc_minus_minc
    gc = (maxc-g) / maxc_minus_minc
    bc = (maxc-b) / maxc_minus_minc

    maxc_is_r = np.equal(maxc,r)
    maxc_is_g = np.equal(maxc,g)
    maxc_is_b = np.equal(maxc,b)

    h = np.zeros((r.shape[0],r.shape[1]))
    h = np.choose( maxc_is_b, (h,4.0+gc-rc) )
    h = np.choose( maxc_is_g, (h,2.0+rc-bc) )
    h = np.choose( maxc_is_r, (h,bc-gc) )

    h = np.mod(h/6.0,1.0)

    hsv = np.asarray([h,s,v])
    
    return hsv

# =============================================================================
# hsv_to_rgb()
#
# hsv comes in as [h,s,v] with hue and saturation in the range [0,1],
# but value in the range [0,255].

def hsv_to_rgb( hsv ):

    h = hsv[0]
    s = hsv[1]
    v = hsv[2]

    #if s == 0.0: return v, v, v
    i = (h*6.0).astype(int)
    f = (h*6.0) - i
    p = v*(1.0 - s)
    q = v*(1.0 - s*f)
    t = v*(1.0 - s*(1.0-f))

    r = i.choose( v, q, p, p, t, v )
    g = i.choose( t, v, v, q, p, p )
    b = i.choose( p, p, t, v, v, q )

    rgb = np.asarray([r,g,b]).astype(np.uint8)
    
    return rgb

# =============================================================================
# Usage()

def Usage():
    print("""Usage: hsv_merge.py [-q] [-of format] src_color src_greyscale dst_color

where src_color is a RGB or RGBA dataset,
      src_greyscale is a greyscale dataset (e.g. the result of gdaldem hillshade)
      dst_color will be a RGB or RGBA dataset using the greyscale as the
      intensity for the color dataset.
""")
    sys.exit(1)
    
# =============================================================================
# 	Mainline
# =============================================================================
def hsv_merge(src_color_filename,src_greyscale_filename,dst_color_filename):
    
    datatype = GDT_Byte
    format = 'GTiff'

    hilldataset = gdal.Open( src_greyscale_filename, GA_ReadOnly )
    colordataset = gdal.Open( src_color_filename, GA_ReadOnly )

    #check for 3 or 4 bands in the color file
    if (colordataset.RasterCount != 3 and colordataset.RasterCount != 4):
        print('Source image does not appear to have three or four bands as required.')
        sys.exit(1)

    #define output format, name, size, type and set projection
    out_driver = gdal.GetDriverByName(format)
    outdataset = out_driver.Create(dst_color_filename, colordataset.RasterXSize, \
                       colordataset.RasterYSize, colordataset.RasterCount, datatype)
    outdataset.SetProjection(hilldataset.GetProjection())
    outdataset.SetGeoTransform(hilldataset.GetGeoTransform())

    #assign RGB and hillshade bands
    rBand = colordataset.GetRasterBand(1)
    gBand = colordataset.GetRasterBand(2)
    bBand = colordataset.GetRasterBand(3)
    if colordataset.RasterCount == 4:
        aBand = colordataset.GetRasterBand(4)
    else:
        aBand = None

    hillband = hilldataset.GetRasterBand(1)
    hillbandnodatavalue = hillband.GetNoDataValue()

    #check for same file size
    if ((rBand.YSize != hillband.YSize) or (rBand.XSize != hillband.XSize)):
        print('Color and hilshade must be the same size in pixels.')
        sys.exit(1)

    #loop over lines to apply hillshade
    for i in range(hillband.YSize):
        #load RGB and Hillshade arrays
        rScanline = rBand.ReadAsArray(0, i, hillband.XSize, 1, hillband.XSize, 1)
        gScanline = gBand.ReadAsArray(0, i, hillband.XSize, 1, hillband.XSize, 1)
        bScanline = bBand.ReadAsArray(0, i, hillband.XSize, 1, hillband.XSize, 1)
        hillScanline = hillband.ReadAsArray(0, i, hillband.XSize, 1, hillband.XSize, 1)

        #convert to HSV
        hsv = rgb_to_hsv( rScanline, gScanline, bScanline )
        
        # if there's nodata on the hillband, use the v value from the color
        # dataset instead of the hillshade value.
        if hillbandnodatavalue is not None:
            equal_to_nodata = np.equal(hillScanline, hillbandnodatavalue)
            v = np.choose(equal_to_nodata,(hillScanline,hsv[2]))
        else:
            v = hillScanline

        #replace v with hillshade
        hsv_adjusted = np.asarray( [hsv[0], hsv[1], v] )
        
        #convert back to RGB
        dst_color = hsv_to_rgb( hsv_adjusted )

        #write out new RGB bands to output one band at a time
        outband = outdataset.GetRasterBand(1)
        outband.WriteArray(dst_color[0], 0, i)
        outband = outdataset.GetRasterBand(2)
        outband.WriteArray(dst_color[1], 0, i)
        outband = outdataset.GetRasterBand(3)
        outband.WriteArray(dst_color[2], 0, i)
        if aBand is not None:
            aScanline = aBand.ReadAsArray(0, i, hillband.XSize, 1, hillband.XSize, 1)
            outband = outdataset.GetRasterBand(4)
            outband.WriteArray(aScanline, 0, i)
        
        #update progress line
        #if not quiet:
        #    gdal.TermProgress_nocb( (float(i+1) / hillband.YSize) )




if __name__ == '__main__':

    gdaldem="C:\\OSGeo4W\\bin\\gdaldem"

    repert="D:\\DEM\\500m_SRTM (CGIAR_2008) Global et Maroc"
    infile="srtm_500m_maroc.tif"


    os.chdir(repert)

    inDs = gdal.Open(infile)
    cols = inDs.RasterXSize
    rows = inDs.RasterYSize
    bands =  inDs.RasterCount
    driver = inDs.GetDriver()
    geo=inDs.GetGeoTransform()
    band=inDs.GetRasterBand(1)
    maxi=band.GetMaximum()
    mini=band.GetMinimum()
    inDs=0
    #mnt=band.ReadAsArray(xoffset,yoffset,width,height)

    #  ----- Hillshade ------------	
    commande1=gdaldem+" hillshade "+infile+" shade.tif -z 5 -s 111120 -az 90"

    # ------- coloriage -----------
    colors=["65535,255,255,255"]
    val=str(maxi)+",254,254,254"
    colors.append(val)
    colors.append("2500, 121, 117 ,10")
    colors.append("1200, 151, 106 ,47")
    colors.append("800, 127, 166 ,122")
    colors.append("400, 213, 213 ,149")
    colors.append("1, 218, 179 ,122")
    colors.append("0, 20, 120 ,200")
    f = open('ramp.txt', 'w')
    for i in range(0,8):
        val=colors[i]+'\n'
        f.write(val)
    f.close()
    commande2=gdaldem+" color-relief "+infile+"  ramp.txt relief.tif"

    # -------- slope et aspect----------

    commande4=gdaldem+" slope "+infile+"  slope.tif -p use percent"
    commande5=gdaldem+" aspect "+infile+"  aspect.tif"

    # -------- os calls...
    print "Calcul du Hillshade : ",commande1
    os.system(commande1)
    print "Coloriage du Hillshade : ",commande2
    os.system(commande2)
    print "Melange des deux : hsvmerge"
    hsv_merge("relief.tif","shade.tif","colour_shade.tif")
    print "Pentes: ",commande4
    os.system(commande4)
    print "Aspects : ",commande5
    os.system(commande5)
    os.system("del ramp.txt")