-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathgs_utils.py
127 lines (92 loc) · 4.57 KB
/
gs_utils.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
# -*- coding: utf-8 -*-
"""gs_utils.ipynb
Automatically generated by Colaboratory.
Original file is located at
https://colab.research.google.com/drive/1vxpUHq43T3KFtz0mc7klguMKVmztulSM
"""
# -*- coding: utf-8 -*-
"""
Created on Tue Nov 30 11:58:48 2021
@author: amullen
"""
import os
os.environ['PROJ_LIB'] = '/adapt/nobackup/people/almullen/.conda/envs/smallsat/share/proj'
os.environ['GDAL_DATA'] = '/adapt/nobackup/people/almullen/.conda/envs/smallsat/share'
from osgeo import gdal, osr, gdalconst
import numpy as np
from skimage import filters
import rasterio
def gaussian_smoothing(image_path, output_path, sigma, input_ndv = None):
with rasterio.open(image_path) as src:
img_data = src.read(np.arange(1,src.count+1).tolist()).astype(float) # rescale bands for "analytic_sr" asset
ndv = src.nodata
if input_ndv != None:
ndv = input_ndv
img_data[img_data==input_ndv] = np.nan
else:
img_data[img_data==ndv] = np.nan
for band_num in range(0,src.count):
img_data[band_num] = filters.gaussian(img_data[band_num], sigma)
with rasterio.open(output_path,
'w',
driver='GTiff',
height=src.height,
width=src.width,
count=src.count,
dtype=rasterio.uint16,
crs=src.crs,
transform=src.transform,
compression = 'lzw',
nodata = ndv) as output:
output.write(img_data[0].astype(np.uint16), 1)
output.write(img_data[1].astype(np.uint16), 2)
output.write(img_data[2].astype(np.uint16), 3)
output.write(img_data[3].astype(np.uint16), 4)
return
def resample(original_source, desired_source, destination_path, nodata_value = -9999):
original_proj = original_source.GetProjection()
original_gt = original_source.GetGeoTransform()
desired_gt = desired_source.GetGeoTransform()
desired_proj = desired_source.GetProjection()
if original_proj != desired_proj or original_gt != desired_gt:
desired_cols = desired_source.RasterXSize
desired_rows = desired_source.RasterYSize
desired_band = desired_source.GetRasterBand(1)
minX = desired_gt[0]
maxX = desired_gt[0] + (desired_cols * desired_gt[1])
minY = desired_gt[3] + (desired_rows * desired_gt[5])
maxY = desired_gt[3]
warp_options = gdal.WarpOptions(format = 'GTiff',
outputBounds = [minX, minY, maxX, maxY],
width = desired_cols, height = desired_rows,
srcSRS = original_proj, dstSRS = desired_proj,
resampleAlg = gdalconst.GRA_Bilinear, outputType = gdalconst.GDT_Float32, creationOptions = ['COMPRESS=LZW'],
dstNodata=nodata_value)
out = gdal.Warp(destination_path, original_source, options = warp_options)
#out.GetRasterBand(1).SetNoDataValue(nodata_value)
desired_band = None
return out
return original_source
def write_geotiff(new_image_name, image_shape, geotransform, projection, dataset, ndv):
driver = gdal.GetDriverByName("GTiff")
outdata = driver.Create(new_image_name, image_shape[1], image_shape[0], 1,
gdalconst.GDT_Float32, ['COMPRESS=LZW'])
outdata.SetGeoTransform(geotransform) ##sets same geotransform as input
srs = osr.SpatialReference() # establish encoding
#srs.ImportFromEPSG(projection) # WGS84 lat/long
srs.ImportFromWkt(projection)
outdata.SetProjection(srs.ExportToWkt())
outdata.GetRasterBand(1).SetNoDataValue(ndv)
outdata.GetRasterBand(1).WriteArray(dataset)
#outdata.GetRasterBand(1).SetNoDataValue(10000)##if you want these values transparent
outdata.FlushCache() ##saves to disk!!
outdata = None
def stack_singleband_raster(raster, existing_stack):
raster_band = raster.GetRasterBand(1)
raster_array = raster_band.ReadAsArray().astype(np.float16)
raster_array[np.isinf(raster_array)] = np.nan
stackedBands = np.append(existing_stack, [raster_array], axis=0)
raster = None
raster_band = None
raster_array = None
return stackedBands