-
Notifications
You must be signed in to change notification settings - Fork 0
/
read_kienholz
executable file
·103 lines (86 loc) · 3.44 KB
/
read_kienholz
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
#!/usr/bin/env python
import gdal
import math
from gdalconst import *
import struct
from netCDF4 import Dataset
import matplotlib.pylab as mpl
import numpy as np
import sys
import glob
increment = 30.
# big_rows = big_cols = 10010
# outfile = Dataset("big_file.nc", 'w', format='NETCDF4')
# #print outfile.file_format
# fdx = outfile.createDimension('x', big_cols)
# fdy = outfile.createDimension('y', big_rows)
# fvx = outfile.createVariable('x','f4',('x',))
# fvy = outfile.createVariable('y','f4',('y',))
# out_height = outfile.createVariable('height','f4',('y','x'), fill_value=-9.e9)
# out_height.coordinates = "y x" ;
# count = outfile.createVariable('count','i4',('y','x'), fill_value=-9.e9)
# count.coordinates = "y x" ;
# count[:] = 0
# # out_height[:]=0.
# x_offset = 1000010.8
# y_offset = 1000028.1
# fvx[:] = np.arange (big_cols) * increment + (x_offset + increment/2.)
# fvy[:] = np.arange (big_rows) * increment + (y_offset + increment/2.)
for filename in glob.glob("*.tif"): # .grid
dataset = gdal.Open(filename, GA_ReadOnly)
driver = dataset.GetDriver().LongName
geotransform = dataset.GetGeoTransform()
band = dataset.GetRasterBand(1)
bandtype = gdal.GetDataTypeName(band.DataType)
scanline = band.ReadRaster( 0, 0, band.XSize, 1,band.XSize, 1, band.DataType)
cols = dataset.RasterXSize
rows = dataset.RasterYSize
bands = dataset.RasterCount
fincrement = - geotransform[-1]
if fincrement != increment:
print "GRID INCREMENTS DON'T MATCH! EXPECTING " + str(increment) + " GOT " + str(fincrement)
print filename
continue
# sys.exit(666)
xll = geotransform [0]
xur = xll + cols * increment # remap to ur corner.
yur = geotransform [3]
yll = yur - rows * increment # remap to ll corner.
data = band.ReadAsArray(0, 0, cols, rows)
#mpl.imshow(data, interpolation = "nearest")
#mpl.colorbar()
#mpl.show()
# print geotransform
xpos = np.arange(cols) * increment + xll + increment *.5
ypos = np.arange(rows) * increment + yll + increment *.5
writesingle = True
if writesingle :
single_outfilename = "test.nc"
single_outfile = Dataset(single_outfilename, 'w', format='NETCDF4')
# print single_outfile.file_format
lon = single_outfile.createDimension('x', cols)
lat = single_outfile.createDimension('y', rows)
latitudes = single_outfile.createVariable('y','f4',('y',))
longitudes = single_outfile.createVariable('x','f4',('x',))
height = single_outfile.createVariable('height','f4',('y','x'))
print xpos.shape
print ypos.shape
print data.shape
longitudes[:] = xpos
latitudes[:] = ypos
height[:] = data[::-1,:]
single_outfile.close()
# print xpos.shape
# print ypos.shape
# print data.shape
# xmin = int(math.floor((xpos[0] - x_offset) / increment))
# xmax = int(xmin + cols)
# ymin = int(math.floor((ypos[0] - y_offset) / increment))
# ymax = int(ymin + rows)
# if xmin > 0 and xmax < big_cols and ymin > 0 and ymax < big_rows:
# # print (xmin, xmax, ymin, ymax)
# out_height[ymin:ymax,xmin:xmax] = np.where(data[::-1,:]> 0 , data[::-1,:], out_height[ymin:ymax,xmin:xmax] )
# count[ymin:ymax,xmin:xmax] += (data [::-1,:] > 0 )
# # out_height[:] = np.where(count[:] > 0 , ( out_height[:] / count[:] ), out_height[:])
# out_height[:]=np.where(count[:],out_height[:],0)
# outfile.close()