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())

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):
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
if hillbandnodatavalue is not None:
equal_to_nodata = np.equal(hillScanline, hillbandnodatavalue)
v = np.choose(equal_to_nodata,(hillScanline,hsv[2]))
else:
v = hillScanline

hsv_adjusted = np.asarray( [hsv[0], hsv[1], v] )

#convert back to RGB

#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

# ------- 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"