-
Notifications
You must be signed in to change notification settings - Fork 1
/
detect_cloud_threshold.py
69 lines (63 loc) · 2.19 KB
/
detect_cloud_threshold.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
#from PIL import Image
from osgeo import gdal
from osgeo import osr
import numpy as np
#setting
src_filename = 'cloud_c_22.tif'
dst_filename = 'mask_cloud/cloud_c_22.tif'
ds = gdal.Open(src_filename)
band = ds.GetRasterBand(1) #matrix band for tif file
gt = ds.GetGeoTransform() #info about geo (used for metadata)
array_band = np.array(band.ReadAsArray())
ysize, xsize = np.shape(array_band)
channel_out = np.zeros(np.shape(array_band)) #same size matrix with value 0
channel_out[array_band >= 240] = 1 #threshold
fileformat = "GTiff"
driver = gdal.GetDriverByName(fileformat)
dst_ds = driver.Create(dst_filename, xsize=xsize, ysize=ysize, bands=1, eType=gdal.GDT_Byte)
#write data for metadata and georefrencing
dst_ds.SetGeoTransform(gt)
wgs84_wkt = """
PROJCRS["WGS_1984_UTM_Zone_33N",
BASEGEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["Transverse Mercator",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",15,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
ID["EPSG",32633]]
"""
srs = osr.SpatialReference()
srs.ImportFromWkt(wgs84_wkt)
srs.SetUTM(33, 1)
srs.SetWellKnownGeogCS("WGS84")
dst_ds.SetProjection(srs.ExportToWkt())
dst_ds.GetRasterBand(1).WriteArray(channel_out)
# Once we're done, close properly the dataset
dst_ds = None