-
Notifications
You must be signed in to change notification settings - Fork 2
/
mgdal.py
132 lines (105 loc) · 3.76 KB
/
mgdal.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
from osgeo import gdal, ogr, osr
import sys,time
import numpy as np
class read_raster():
def __init__(self, raster):
self.raster = raster
self.graster = gdal.Open(self.raster)
self.band = self.graster.GetRasterBand(1)
self.stats = self.band.GetStatistics( True, True )
def stat(self):
return {'min':self.stats[0],'max':self.stats[1],'mean':self.stats[2],'std':self.stats[3]}
def tonp(self):
return self.band.ReadAsArray().astype('float32')
def getnodata(self):
ds = gdal.Open(self.raster)
band = ds.GetRasterBand(1)
nodata = band.GetNoDataValue()
myarray1 = np.array(band.ReadAsArray()).astype('float32')
myarray1[myarray1 == nodata] = np.nan
return myarray1
def array2raster(self, rasterfn,newRasterfn,array,nodata):
raster = gdal.Open(rasterfn)
geotransform = raster.GetGeoTransform()
originX = geotransform[0]
originY = geotransform[3]
pixelWidth = geotransform[1]
pixelHeight = geotransform[5]
cols = raster.RasterXSize
rows = raster.RasterYSize
driver = gdal.GetDriverByName('GTiff')
outRaster = driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Float32) #gdal.GDT_Byte or gdal.GDT_Float32
outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
outband = outRaster.GetRasterBand(1)
outband.SetNoDataValue(nodata)
outband.WriteArray(array)
outRasterSRS = osr.SpatialReference()
outRasterSRS.ImportFromWkt(raster.GetProjectionRef())
outRaster.SetProjection(outRasterSRS.ExportToWkt())
outband.FlushCache()
def raster2array(raster):
return myraster_to_array(raster)
def myraster_to_array(raster_path1):
ds = gdal.Open(raster_path1)
band = ds.GetRasterBand(1)
nodata = band.GetNoDataValue()
myarray1 = np.array(band.ReadAsArray()).astype('float32')
myarray1[myarray1 == nodata] = np.nan
return myarray1
def array2raster(rasterfn,newRasterfn,array):
raster = gdal.Open(rasterfn)
band = raster.GetRasterBand(1)
nodata = band.GetNoDataValue()
geotransform = raster.GetGeoTransform()
originX = geotransform[0]
originY = geotransform[3]
pixelWidth = geotransform[1]
pixelHeight = geotransform[5]
cols = raster.RasterXSize
rows = raster.RasterYSize
driver = gdal.GetDriverByName('GTiff')
outRaster = driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Float32) #gdal.GDT_Byte or gdal.GDT_Float32
outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
outband = outRaster.GetRasterBand(1)
outband.SetNoDataValue(nodata)
outband.WriteArray(array)
outRasterSRS = osr.SpatialReference()
outRasterSRS.ImportFromWkt(raster.GetProjectionRef())
outRaster.SetProjection(outRasterSRS.ExportToWkt())
outband.FlushCache()
def write_info(file):
import gdalinfo2
gdalinfo2.main(file)
def reclassify(file,classes):
rArray = myraster_to_array(file)
a = rArray
classes[-1] = [classes[-1][0],classes[-1][1] + 1,classes[-1][2]]
for x in classes:
#a[np.where((a >= x[0]) & (a < x[1]))] =x[2]
a[((a >= x[0]) & (a < x[1]))] = x[2]
return a
if __name__ == '__main__':
file = r"bb.tif"
classes =[[0,500,1.0],[2500,3000,2.0],[3000,4000,3.0]]
a = reclassify(file,classes)
array2raster(file,"rrr.tif",a)
write_info("rrr.tif")
input()
greader = read_raster('a.tif')
min = greader.stat()['min']
max = greader.stat()['max']
print min,max
m='''
rasterfn = 'Slope.tif'
newValue = 0
newRasterfn = 'SlopeNew.tif'
# Convert Raster to array
rasterArray = raster2array(rasterfn)
# Get no data value of array
noDataValue = getNoDataValue(rasterfn)
# Updata no data value in array with new value
rasterArray[rasterArray == noDataValue] = newValue
# Write updated array to new raster
array2raster(rasterfn,newRasterfn,rasterArray)
'''
input()