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

# Michel Le Page, 06/04/2013

from math import *
import numpy as np
#import scipy as sp
from osgeo import gdal
import csv
from osgeo.gdalconst import *


if __name__ == "__main__":

	# Fichiers
	infile="D:\DEM\SRTM_500m1"
	outfile1="D:\DEM\srtm_500m_maroc.tif"

	xmin=-11.
	ymin=36.
	larg=9.
	haut=8.

	inDs = gdal.Open(infile)
	cols = inDs.RasterXSize
	rows = inDs.RasterYSize
	bands =  inDs.RasterCount
	driver = inDs.GetDriver()
	geo=inDs.GetGeoTransform()

	xoffset=int((xmin-geo[0])/geo[1])
	width=int(larg/geo[1])
	yoffset=int((geo[3]-ymin)/(-geo[5]))
	height=int(haut/(-geo[5]))
	
	print "Input Cols:", cols, "Ligs:", rows,driver
	print "Output Cols:", width, "Ligs:", height

	band=inDs.GetRasterBand(1)
	  
	mnt=band.ReadAsArray(xoffset,yoffset,width,height)

	format = "GTiff"
	driver = gdal.GetDriverByName( format )
	outDs=driver.Create(outfile1,width,height,bands,GDT_UInt16)
	geo=inDs.GetGeoTransform()
	# changement des coordonées
	lgeo=list(geo)
	lgeo[0]=geo[0]+(xoffset*geo[1])
	lgeo[3]=geo[3]+(yoffset*geo[5])
	tgeo=tuple(lgeo)
	print tgeo
	outDs.SetGeoTransform(tgeo)
	outDs.SetProjection(inDs.GetProjection())
	print inDs.GetProjection()
	bando=outDs.GetRasterBand(1)
	bando.WriteArray(mnt)
	outDs=None
	
	inDs=None