Commit 51a3052a69d1ecf916e45fed9f5b041dd2e5b922

Authored by michel-lepage
0 parents
Exists in master

init

dem_modcou.py 0 → 100644
  1 +++ a/dem_modcou.py
... ... @@ -0,0 +1,94 @@
  1 +#! /usr/bin/env python
  2 +# -*- coding: iso-8859-1 -*-
  3 +
  4 +# Michel Le Page, 06/04/2013
  5 +#
  6 +# Create a unregular matrix from a DEM file for Modcou model
  7 +# Input: A DEM file with a pixel size wihch is multiple of 1000 meters (eg 500, 250, 125, 100...)
  8 +# Outputs: Two filse. One is an index file that can be used within a GIS app' to be vectorize. The second one if the average height of each new grid.
  9 +
  10 +from math import *
  11 +import numpy as np
  12 +#import scipy as sp
  13 +from osgeo import gdal
  14 +import csv
  15 +from osgeo.gdalconst import *
  16 +
  17 +
  18 +if __name__ == "__main__":
  19 +
  20 + # Fichiers
  21 + infile="D:\DEM\GMTED2010\Gmted_250m_maroc.tif"
  22 + outfile1="D:\DEM\GMTED2010\Gmted_250m_maroc_surfexind.tif"
  23 + outfile2="D:\DEM\GMTED2010\Gmted_250m_maroc_surfexval.tif"
  24 +
  25 + resol=250 #resolution en metres du DEM d'entree, attention, doit etre un multiple de 1000m (500, 250, 125, 100...)
  26 +
  27 + pasmax=8000/resol # pas en pixel pour 8km
  28 + pasmin=1000/resol # pas en pixel pour 1km
  29 + stdevmin=30 # min stdev pour regrouper des pixels
  30 + resols=(pasmax/8,pasmax/4,pasmax/2,pasmax) # résolutions intermediaires à 1, 2, 4 et 8 km, attention l'ordre du + petit au + grand est important
  31 +
  32 + inDs = gdal.Open(infile)
  33 + cols = inDs.RasterXSize
  34 + rows = inDs.RasterYSize
  35 + bands = inDs.RasterCount
  36 + driver = inDs.GetDriver()
  37 +
  38 + print "Cols:", cols, "Ligs:", rows, "Bands:", bands,driver
  39 +
  40 +
  41 + band=inDs.GetRasterBand(1)
  42 +
  43 + mnt=band.ReadAsArray()
  44 +
  45 + matout=np.zeros([rows,cols], np.uint32)
  46 + mntout=np.zeros([rows/pasmin,cols/pasmin], np.float32)
  47 +
  48 + val=1L
  49 +
  50 + for i in range(0,rows,pasmax):
  51 + if i%256==0: print i
  52 + i1=i/pasmin
  53 + for j in range(0,cols,pasmax):
  54 + j1=j/pasmin
  55 + # pour les 4 résolutions intermédiaires, on teste de la plus petite à la plus grande resol
  56 + for k in resols:
  57 + k1=k/pasmin
  58 + for ii in range(0,pasmax,k):
  59 + ii1=ii/pasmin
  60 + for jj in range(0,pasmax, k):
  61 + stdev=np.std(mnt[i+ii:i+ii+k,j+jj:j+jj+k])
  62 + if stdev<stdevmin or (k==pasmin and stdev>=stdevmin):
  63 + #print "ok:", (i*pasmax)+ii,(i*pasmax)+ii+k,(j*pasmax)+jj,(j*pasmax)+jj+k, val
  64 + matout[i+ii:i+ii+k,j+jj:j+jj+k]=val
  65 + mean=np.mean(mnt[i+ii:i+ii+k,j+jj:j+jj+k])
  66 + mntout[i1+ii1:i1+ii1+k1,j1+jj/pasmin:j1+jj/pasmin+k1]=mean
  67 + val=val+1
  68 + #else:
  69 + # print i,j,stdev, mnt[i+ii:i+ii+k,j+jj:j+jj+k]
  70 +
  71 +
  72 + outDs=driver.Create(outfile1,cols,rows,bands,GDT_UInt32)
  73 + outDs.SetGeoTransform(inDs.GetGeoTransform())
  74 + outDs.SetProjection(inDs.GetProjection())
  75 + bando=outDs.GetRasterBand(1)
  76 + bando.WriteArray(matout)
  77 + outDs=None
  78 +
  79 + outDs=driver.Create(outfile2,cols/pasmin,rows/pasmin,bands,GDT_Float32)
  80 + geo=inDs.GetGeoTransform()
  81 + # changement de la résolution de sortie à 1km
  82 + lgeo=list(geo)
  83 + lgeo[1]*=pasmin
  84 + lgeo[5]*=pasmin
  85 + tgeo=tuple(lgeo)
  86 + outDs.SetGeoTransform(tgeo)
  87 + outDs.SetProjection(inDs.GetProjection())
  88 + bando=outDs.GetRasterBand(1)
  89 + bando.WriteArray(mntout)
  90 + outDs=None
  91 +
  92 + inDs=None
  93 +
  94 +
... ...
dem_shp.py 0 → 100644
  1 +++ a/dem_shp.py
... ... @@ -0,0 +1,94 @@
  1 +#! /usr/bin/env python
  2 +# -*- coding: iso-8859-1 -*-
  3 +
  4 +# Michel Le Page, 08/04/2013
  5 +#
  6 +# This code calculate some simple stats from an image file (in this case a DEM) related to a shapefile
  7 +# input: a rasterfile (repert and infile) and a shapefile (shapefile) with its ID column (attribut)
  8 +# output: a text file with, id,number,average,stdev,min,max
  9 +#
  10 +
  11 +
  12 +import os
  13 +from math import *
  14 +import numpy as np
  15 +from osgeo import gdal
  16 +import csv
  17 +from osgeo.gdalconst import *
  18 +from osgeo import gdal, gdal_array
  19 +import sys
  20 +
  21 +if __name__ == '__main__':
  22 +
  23 + gdaldem="C:\\OSGeo4W\\bin\\gdaldem"
  24 + gdal_rasterize="C:\\OSGeo4W\\bin\\gdal_rasterize"
  25 +
  26 + repert="D:\\DEM\\500m_SRTM (CGIAR_2008) Global et Maroc"
  27 + infile="srtm_500m_maroc.tif"
  28 +
  29 + shapefile="d:\\shp\\communes_maroc_geo.shp"
  30 + attribut="COMMUNE_ID"
  31 + layer="communes_maroc_geo"
  32 +
  33 + os.chdir(repert)
  34 +
  35 + inDs = gdal.Open(infile)
  36 + cols = inDs.RasterXSize
  37 + rows = inDs.RasterYSize
  38 + bands = inDs.RasterCount
  39 + driver = inDs.GetDriver()
  40 + geo=inDs.GetGeoTransform()
  41 + print geo
  42 + proj=inDs.GetProjection()
  43 + print proj
  44 + band=inDs.GetRasterBand(1)
  45 + maxi=band.GetMaximum()
  46 + mini=band.GetMinimum()
  47 + inDs=0
  48 +
  49 +
  50 + # sauvergarder la projection
  51 + #f = open('proj.txt', 'w')
  52 + #f.write(proj)
  53 + #f.close()
  54 + # rasteriser le fichier shape
  55 + commande=gdal_rasterize+" -a "+attribut+" -ts "+str(cols)+" "+str(rows)
  56 + commande+=" -l "+layer+" "+shapefile+" tmp.tif -ot Int16"
  57 + #commande+=" -a_srs proj.txt "
  58 + commande+=" -te "+str(geo[0])+" "+str(geo[3]+geo[5]*rows)+" "+str(geo[0]+geo[1]*cols)+" "+str(geo[3])
  59 +
  60 + print commande
  61 +
  62 + #os.system(commande)
  63 +
  64 + inDs = gdal.Open(infile)
  65 + band=inDs.GetRasterBand(1)
  66 + mnt=band.ReadAsArray()
  67 +
  68 + inDs2 = gdal.Open("tmp.tif")
  69 + driver2 = inDs2.GetDriver()
  70 + band2=inDs2.GetRasterBand(1)
  71 + shp=band2.ReadAsArray()
  72 + maxishp=shp.max()
  73 + minishp=shp.min()
  74 + print minishp,maxishp
  75 +
  76 + f = open('stats.txt', 'w')
  77 + out= "i,n,moy,std,min,max\n"
  78 + f.write(out)
  79 +
  80 + aa=np.unique(shp)
  81 +
  82 + for i in aa:
  83 + a=np.where(shp==i)
  84 + if mnt[a].size>0:
  85 + moy=mnt[a].mean()
  86 + std=mnt[a].std()
  87 + min=mnt[a].min()
  88 + max=mnt[a].max()
  89 + out= str(i)+","+str(mnt[a].size)+","+str(moy)+","+str(std)+","+str(min)+","+str(max)+"\n"
  90 + f.write(out)
  91 + if (i%100)==0:print(i)
  92 +
  93 + f.close()
  94 +
... ...
dem_winstral.py 0 → 100644
  1 +++ a/dem_winstral.py
... ... @@ -0,0 +1,119 @@
  1 +#! /usr/bin/env python
  2 +# -*- coding: iso-8859-1 -*-
  3 +
  4 +# Michel Le Page, 20/04/2013
  5 +# Winstral, Spatial Snow Modeling of Wind-Redistributed Snow Using Terrain-Based Parameters
  6 +# Journal of Hydrometeorology, 2002
  7 +# seules les eq 1 et 2 sont calculées.
  8 +# Input : un DEM
  9 +# Output: average of n maximum upwind slopes relative to seasonally averaged winds (TIF Float16)
  10 +#
  11 +# Python 3.3
  12 +
  13 +
  14 +from math import *
  15 +import numpy as np
  16 +from osgeo import gdal
  17 +from osgeo.gdalconst import *
  18 +import sys
  19 +
  20 +
  21 +if __name__ == "__main__":
  22 +
  23 + # ======= DEBUT CONFIG ==============
  24 + # Fichiers
  25 + infile="testDEM.tif"
  26 + outfile1="testDEM_winstral.tif"
  27 +
  28 + resol=30. #resolution en metres du DEM d'entree
  29 +
  30 + azimuth_vent=90. # azimuth du vent dominant en degrés
  31 + az1=azimuth_vent-15. #azimuth dominant - un certain angle (-15° par défaut)
  32 + az2=azimuth_vent+15. #azimuth dominant + un certain angle (+15° par défaut)
  33 + inc=5. # incrément de calcul des droites de Vent (5° par défaut)
  34 + dmax=500 # distance dmax
  35 +
  36 + # ======= FIN CONFIG ==============
  37 +
  38 + nbvect=int(((az2-az1)/inc)+1) # Nombre de vecteurs= nombre de directions de calcul (fig3, p528)
  39 + longvect=int(dmax/resol) # Longueur des vecteurs en fonction du param dmax.
  40 +
  41 + # lecture du fichier d'entree
  42 + inDs = gdal.Open(infile)
  43 + cols = inDs.RasterXSize
  44 + rows = inDs.RasterYSize
  45 + bands = inDs.RasterCount
  46 + driver = inDs.GetDriver()
  47 + print("Cols:", cols, "Ligs:", rows, "Bands:", bands,driver)
  48 + band=inDs.GetRasterBand(1)
  49 + mnt=band.ReadAsArray()
  50 + mnt=np.float32(mnt)
  51 +
  52 + # on passe tout en radians
  53 + azimuth_vent=azimuth_vent/180*pi
  54 + az1=az1/180.*pi
  55 + az2=az2/180.*pi
  56 + inc=inc/180.*pi
  57 +
  58 + # vect[][0]: les x, vect[][1]=mles y, vect[][2]=les dists
  59 + vects=np.zeros((nbvect,3,longvect),np.float32)
  60 +
  61 + # Calcul des coordonnées pixel des droites.
  62 +
  63 + for n in range(0,nbvect):
  64 + az=az1+(n*inc)
  65 + print("Azimuth=",az/pi*180)
  66 + x=dmax*cos(az)
  67 + y=dmax*sin(az)
  68 + pixx=x/dmax
  69 + pixy=y/dmax
  70 + print(x,y,pixx,pixy)
  71 +
  72 + for xx in range(1,longvect+1):
  73 + vects[n][0][xx-1]=round(xx*pixx)
  74 + vects[n][1][xx-1]=round(xx*pixy)
  75 + vects[n][2][xx-1]=sqrt(pow(float(vects[n][0][xx-1])*resol,2)+pow(float(vects[n][1][xx-1])*resol,2))
  76 + print(vects[n][0][xx-1],vects[n][1][xx-1],vects[n][2][xx-1])
  77 +
  78 + maxx=max(0,int(np.max(vects[:,0,:]))+1)
  79 + maxy=max(0,int(np.max(vects[:,1,:]))+1)
  80 + minx=max(0,(int(np.min(vects[:,0,:]))-1)*-1)
  81 + miny=max(0,(int(np.min(vects[:,1,:]))-1)*-1)
  82 + print(minx,miny,maxx,maxy)
  83 +
  84 + # matrice de sortie
  85 + mntout=np.zeros((rows,cols),np.float32)
  86 +
  87 + # stockage temporaire des pentes
  88 + sx=np.zeros(nbvect)
  89 +
  90 + #calculs pour chaque pixel
  91 + hash10=range(0,rows,int(rows/10))
  92 + hash100=range(0,rows,int(rows/100))
  93 + maxyy=rows-maxy-1
  94 + maxxx=cols-maxx-1
  95 + for j in range(miny,maxyy,1):
  96 + if j in hash10: print(hash10.index(j)*10,'%',end="")
  97 + if j in hash100: print('.',end="")
  98 + sys.stdout.flush()
  99 + sx=sx*0.
  100 + for i in range(minx,maxxx,1):
  101 + z=mnt[j,i]
  102 + # calcule les sx pour les 7 (nbvect) droites.
  103 + for n in range(0,nbvect):
  104 + alt=mnt[np.int32(vects[n][1]+j),np.int32(vects[n][0]+i)] # numérateur de eq1 p528
  105 + sx[n]=np.max(np.tan((alt-z)/vects[n][2])) # eq 1 p528
  106 + mntout[j,i]=np.average(sx) # eq 2 p 529 (ou comment écrire quelque chose de simple de manière compliquée!)
  107 +
  108 + mntout=np.arctan(mntout)/pi*180
  109 + #ecriture de mnt out
  110 + outDs=driver.Create(outfile1,cols,rows,bands,GDT_Float32)
  111 + outDs.SetGeoTransform(inDs.GetGeoTransform())
  112 + outDs.SetProjection(inDs.GetProjection())
  113 + bando=outDs.GetRasterBand(1)
  114 + bando.WriteArray(mntout)
  115 + outDs=None
  116 +
  117 + # done! print "done!"
  118 +
  119 +
... ...
et0_hargeaves.py 0 → 100644
  1 +++ a/et0_hargeaves.py
... ... @@ -0,0 +1,33 @@
  1 +# ------------------------------------------------------------------
  2 +# Compute Reference Evapotranspiration ET0 with Hargreaves equation on a daily time step
  3 +# Michel Le Page, 2011
  4 +#
  5 +# When solar radiation data, relative humidity data and/or wind speed data are missing, they should be estimated using the procedures presented in this section. As an alternative, ETo can be estimated using the Hargreaves ETo equation
  6 +#
  7 +# Hargreaves, G.L., Hargreaves, G.H., and Riley, J.P. 1985. Agricultural benefits for Senegal River Basin. J. Irrigation and Drainage Engr., ASCE 111:113-124.
  8 +#
  9 +# day1: julian day from 0 to 366
  10 +# alt : altitude of the meteo station (meters)
  11 +# tmoy: average temperature of the day (Celsius Degrees)
  12 +# tmin: minimum temperature of the day (Celsius Degrees)
  13 +# tmax: maximum temperature of the day (Celsius Degrees)
  14 +#
  15 +# ------------------------------------------------------------------
  16 +
  17 +from math import *
  18 +
  19 +def et0_hargreaves(day1,lat,tmoy,tmin,tmax):
  20 + if (tmin==tmax): return(-1)
  21 + try:
  22 + conv_rad=lat * 3.1416 / 180
  23 + const_sol=0.082
  24 + p=3.142
  25 + dr= 1+0.033*cos(2*pi*day1/365)
  26 + d= 0.409*sin((2*pi*day1/365)-1.39)
  27 + ws= acos(-tan(conv_rad)*tan(d))
  28 + ra= (24*60/pi)*const_sol*dr*(ws*sin(conv_rad)*sin(d)+cos(conv_rad)*cos(d)*sin(ws))
  29 + ra=ra/2.45
  30 + et0=0.0023*(tmoy+17.8)*sqrt(tmax-tmin)*ra
  31 + return(et0)
  32 + except:
  33 + return(-1)
0 34 \ No newline at end of file
... ...
et0_pm.py 0 → 100644
  1 +++ a/et0_pm.py
... ... @@ -0,0 +1,53 @@
  1 +# ------------------------------------------------------------------
  2 +# Compute Reference Evapotranspiration ET0 with Penman-Monteith equation on a daily time step
  3 +# Michel Le Page, 2011
  4 +#
  5 +# The FAO-56 Penman-Monteith equation (Allen et al., 1998) estimates the reference crop evapotranspiration, ET0 (mm d-1)
  6 +#
  7 +# Monteith, J.L., 1965. Evaporation and Environment. 19th Symposia of the Society for Experimental Biology, University Press, Cambridge, 19:205-234.
  8 +#
  9 +# day1: julian day from 0 to 366
  10 +# alt : altitude of the meteo station (meters)
  11 +# hmes: height of measurement (meters)
  12 +# lat : latitude of the meteo station (decimal degrees)
  13 +# tmoy: average temperature of the day (Celsius Degrees)
  14 +# tmin: minimum temperature of the day (Celsius Degrees)
  15 +# tmax: maximum temperature of the day (Celsius Degrees)
  16 +# vv: Average Wind Speed of the day (m/s)
  17 +# hrmoy: average Relative Humidity of the day (%)
  18 +# hrmin: minimum Relative Humidity of the day (%)
  19 +# hrmax: maximum Relative Humidity of the day (%)
  20 +# rs: Solar Radiation averaged for the day (W/m2)
  21 +#
  22 +# ------------------------------------------------------------------
  23 +
  24 +from math import *
  25 +
  26 +def et0_pm(day1,alt,hmes,lat,tmoy,tmin,tmax,vv,hrmoy,hrmin,hrmax,rs):
  27 + if ((hrmax-hrmin)<.1 or tmin==tmax): return(-1)
  28 + try:
  29 + conv_rad=lat * pi / 180.
  30 + # -- constantes --
  31 + const_sol=0.082
  32 + const_stef=4.9*pow(10,-9)
  33 + g=0.063
  34 +
  35 + u2= vv * 4.87 / log(67.8*hmes-5.42)
  36 + rs_mj= rs*24*3600*0.000001
  37 + dr= 1+0.033*cos(2*pi*day1/365)
  38 + d= 0.409*sin((2*pi*day1/365)-1.39)
  39 + ws= acos(-tan(conv_rad)*tan(d))
  40 + ra= (24*60/pi)*const_sol*dr*(ws*sin(conv_rad)*sin(d)+cos(conv_rad)*cos(d)*sin(ws))
  41 + rso= ra*(0.75+0.00002*alt)
  42 + f= 1.35*(rs_mj/rso)-0.35
  43 + es= (0.6108*exp(17.27*tmin/(tmin+237.3))+0.6108*exp(17.27*tmax/(tmax+237.3)))/2
  44 + ea= (hrmin*0.6108*exp(17.27*tmax/(tmax+237.3))+hrmax*0.6108*exp(17.27*tmin/(tmin+237.3)))/(2*100)
  45 + rnl= const_stef*0.5*(pow((tmin+273),4)+pow((tmax+273),4))*(0.34-0.14*sqrt(ea))*f
  46 + rn=(1-0.23)*rs_mj-rnl
  47 +
  48 + gflux=0
  49 + delta= 4098*0.6108*exp(17.27*tmoy /(tmoy +237.3))/pow((tmoy +237.3),2)
  50 + et0= (0.408*delta*(rn)+(900*g/(tmoy+273))*u2*(es-ea))/(delta+g*(1+0.34*u2))
  51 + return(et0)
  52 + except:
  53 + return(-1)
0 54 \ No newline at end of file
... ...
et0_priestley_taylor.py 0 → 100644
  1 +++ a/et0_priestley_taylor.py
... ... @@ -0,0 +1,48 @@
  1 +# ------------------------------------------------------------------
  2 +# Compute Reference Evapotranspiration ET0 with Priestley Taylor method on a daily time step
  3 +# Michel Le Page, 2011
  4 +
  5 +# The Priestley-Taylor method (Priestley and Taylor, 1972) for the calculation of daily ET0 (mm d-1) replaces the aerodynamic term of Penman-Monteith equation by a dimensionless empirical multiplier (a, Priestley-Taylor coefficient)
  6 +#
  7 +# day1: julian day from 0 to 366
  8 +# alt : altitude of the meteo station (meters)
  9 +# lat : latitude of the meteo station (decimal degrees)
  10 +# tmoy: average temperature of the day (Celsius Degrees)
  11 +# tmin: minimum temperature of the day (Celsius Degrees)
  12 +# tmax: maximum temperature of the day (Celsius Degrees)
  13 +# hrmin: minimum Relative Humidity of the day (%)
  14 +# hrmax: maximum Relative Humidity of the day (%)
  15 +# rs: Solar Radiation averaged for the day (W/m2)
  16 +#
  17 +# ------------------------------------------------------------------
  18 +
  19 +from math import *
  20 +
  21 +def et0_priestley_taylor(day1,alt,lat,tmoy,tmin,tmax,hrmin,hrmax,rs):
  22 + if ((hrmax-hrmin)<.1 or tmin==tmax): return(-1)
  23 + try:
  24 + lambda=2.45
  25 + conv_rad=lat * 3.1416 / 180
  26 + pt_coeff=1.26
  27 + const_sol=0.082
  28 + g=0.063
  29 + gflux=0
  30 + p=3.142
  31 + const_stef=4.9*pow(10,-9)
  32 +
  33 + rs_mj= rs*24*3600*0.000001
  34 + dr= 1+0.033*cos(2*p*day1/365)
  35 + d= 0.409*sin((2*p*day1/365)-1.39)
  36 + ws= acos(-tan(conv_rad)*tan(d))
  37 + ra= (24*60/p)*const_sol*dr*(ws*sin(conv_rad)*sin(d)+cos(conv_rad)*cos(d)*sin(ws))
  38 + rso= ra*(0.75+0.00002*alt)
  39 + ea= (hrmin*0.6108*e(17.27*tmax/(tmax+237.3))+hrmax*0.6108*e(17.27*tmin/(tmin+237.3)))/(2*100)
  40 + f= 1.35*(rs_mj/rso)-0.35
  41 + rnl= const_stef*0.5*(pow((tmin+273),4)+pow((tmax+273),4))*(0.34-0.14*sqrt(ea))*f
  42 + rn= (1-0.23)*rs_mj-rnl
  43 + delta= 4098*0.6108*e(17.27*tmoy /(tmoy +237.3))/pow((tmoy +237.3),2)
  44 +
  45 + et0=pt_coeff*(delta/(delta+g))*(rn/lambda)
  46 + return(et0)
  47 + except:
  48 + return(-1)
0 49 \ No newline at end of file
... ...
netcdf2tif_GRACE.py 0 → 100644
  1 +++ a/netcdf2tif_GRACE.py
... ... @@ -0,0 +1,122 @@
  1 +# -*- coding: utf-8 -*-
  2 +"""
  3 +Created on Tue Nov 12 15:29:18 2013
  4 +
  5 +@author: Michel Le Page
  6 +
  7 +GRACE Total Water Content data was downloaded from http://www.esrl.noaa.gov/psd/data/gridded/data.gpcc.html
  8 +Three different solutions (GFZ, CSR, JPL) are given into netcdf files at a 1 degree resolution (360x180 pixels) for the period 2003 to present, at a monthly timestep. The matrix is oriented from 90S à 90N and 0 to 360 degrees
  9 +
  10 +This code extract a sub-area and a sub-period, and save it to a TIF file
  11 +
  12 +INPUT: infiles (netcdf GRACE files)
  13 +PARAMS: year_begin, year_end, lat_sud, lat_nord, lon_ouest, lon_est,
  14 +OUTPUT: Three TIF files
  15 +"""
  16 +
  17 +from math import *
  18 +import numpy as np
  19 +from osgeo import gdal
  20 +from osgeo.gdalconst import *
  21 +import sys
  22 +import scipy.ndimage
  23 +
  24 +
  25 +if __name__ == "__main__":
  26 +
  27 + # params des fichiers d'entree
  28 +
  29 + os.chdir("C:\\Users\\michel\\Documents\\GRACE\\med_dataset\\GRACE\\netcdf")
  30 +
  31 + infile_sf="CLM4.RL05.SCALE_FACTOR.DS.G200KM.nc"
  32 + inDs_sf=gdal.Open('NETCDF:"'+infile_sf+'":SCALE_FACTOR')
  33 + cols = inDs_sf.RasterXSize
  34 + rows = inDs_sf.RasterYSize
  35 + bands = inDs_sf.RasterCount
  36 + bandi = inDs_sf.GetRasterBand(1)
  37 + scale_factor=bandi.ReadAsArray(0,0, cols, rows)
  38 + inDs_sf=None
  39 +
  40 + infiles=["GRACE.CSR.LAND.RL05.DS.G200KM.nc","GRACE.GFZ.LAND.RL05.DS.G200KM.nc","GRACE.JPL.LAND.RL05.DS.G200KM.nc"]
  41 + outfiles=["..\\GRACE.CSR.LAND.RL05.DS.G200KM_MEDIT_2003_2012.tif","..\\GRACE.GFZ.LAND.RL05.DS.G200KM_MEDIT_2003_2012.tif","..\\GRACE.JPL.LAND.RL05.DS.G200KM_MEDIT_2003_2012.tif"]
  42 + rows=180
  43 + cols=360
  44 + resolution=1
  45 +
  46 + # facteur du resampling
  47 +
  48 + resize_factor=1
  49 +
  50 + # params de la zone à extraire
  51 + year_begin=2003
  52 + year_end=2013
  53 + lat_sud=24.
  54 + lat_nord=48.
  55 + lon_ouest=-15
  56 + lon_est=40.
  57 +
  58 + # attention la matrice GRACE va de 90S à 90N et de 0 à 360°
  59 +
  60 + x0_ouest=int((360+lon_ouest)/resolution)
  61 + x0_est=cols
  62 +
  63 + x1_ouest=0
  64 + x1_est=int(lon_est/resolution)
  65 +
  66 + y_sud=int((90-lat_sud)/resolution)
  67 + y_nord=int((90-lat_nord)/resolution)
  68 +
  69 + rows_out=int((y_sud-y_nord)*resize_factor)
  70 + cols_out=int(((x0_est-x0_ouest)+(x1_est-x1_ouest))*resize_factor)
  71 + nbband_out=int((year_end-year_begin)*12)
  72 +
  73 +
  74 + for ifile in range(0,3):
  75 + infile=infiles[ifile]
  76 + outfile=outfiles[ifile]
  77 + print(infile)
  78 +
  79 + inDs=gdal.Open('NETCDF:"'+infile+'":lwe_thickness')
  80 +
  81 + cols = inDs.RasterXSize
  82 + rows = inDs.RasterYSize
  83 + bands = inDs.RasterCount
  84 +
  85 + driver = gdal.GetDriverByName("GTiff")
  86 + #outDs=driver.Create("GPCC_precip.mon.total.v6.tif",37+98,58,bands,GDT_Float32)
  87 + outDs=driver.Create(outfile,cols_out,rows_out,nbband_out,GDT_Float32)
  88 + #outDs.SetGeoTransform(inDs.GetGeoTransform())
  89 + #outDs.SetProjection(inDs.GetProjection())
  90 +
  91 + band_begin=(year_begin-2003)*12
  92 + band_end=(year_end-2003)*12
  93 +
  94 + for band in range(band_begin,band_end):
  95 + bandi=inDs.GetRasterBand(band+1)
  96 + raster=bandi.ReadAsArray(0,0, cols, rows)
  97 +
  98 + raster=raster*scale_factor
  99 +
  100 + rasterout=np.zeros((rows_out/resize_factor,cols_out/resize_factor),np.float32)
  101 +
  102 + for j in range(y_nord,y_sud):
  103 + for i in range(x0_ouest,x0_est):
  104 + rasterout[j-y_nord,i-x0_ouest]=raster[j,i]
  105 + for i in range(x1_ouest,x1_est):
  106 + rasterout[j-y_nord,i+(x0_est-x0_ouest)]=raster[j,i]
  107 +
  108 + bando=outDs.GetRasterBand(band-band_begin+1)
  109 +
  110 + if resize_factor > 1:
  111 + raster_resample=scipy.ndimage.zoom(rasterout,resize_factor,order=3)
  112 + bando.WriteArray( raster_resample )
  113 + if resize_factor < 1:
  114 + raster_resample=block_reduce(rasterout, block_size=(1./resize_factor, 1./resize_factor), func=np.mean)
  115 + bando.WriteArray( raster_resample )
  116 + if resize_factor == 1:
  117 + bando.WriteArray( rasterout )
  118 +
  119 +
  120 + outDs=None
  121 +
  122 +
0 123 \ No newline at end of file
... ...