Skip to content

Commit e54b0af

Browse files
wajuqiwajuqi
committed
python geoprocessing tools
0 parents  commit e54b0af

18 files changed

+785
-0
lines changed

‎call_gdal.py

+29
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,29 @@
1+
#from osgeo import gdal
2+
#import os
3+
##clip the image by shp (call gdalwarp utility)
4+
##import subprocess
5+
#shp = "Barrow_two.shp"
6+
#INPUT = 'Barrow_msk_float32.tif'
7+
#OUTPUT = 'gdal output5.tif'
8+
##os.system ('gdalwarp -cutline %s -crop_to_cutline -dstalpha %s %s' %(shp, INPUT, OUTPUT))
9+
#subprocess.call(['gdalwarp', '-cutline', shp, '-crop_to_cutline','-dstalpha', INPUT, OUTPUT]) #-dstalpha:add an alpha band to the output tiff which masks out the area falling outside the cutline.
10+
11+
from osgeo import gdal
12+
import ogr
13+
from numpy import *
14+
import os
15+
import subprocess
16+
17+
#get lake boundary extent
18+
polygon = 'Barrow.shp'
19+
source_ds = ogr.Open(polygon)
20+
source_layer = source_ds.GetLayer()
21+
x_min, x_max, y_min, y_max = source_layer.GetExtent()
22+
print("lake boundary extent:", x_min, x_max, y_min, y_max)
23+
24+
25+
#subset the image
26+
INPUT = 'subset_1_of_S1B_EW_GRDM_1SDH_20161206T032144_20161206T032244_003270_005937_8206.tif'
27+
OUTPUT = "subset2.tif"
28+
translate = 'gdal_translate -projwin %s %s %s %s %s %s' %(x_min, y_max, x_max, y_min, INPUT, OUTPUT)
29+
os.system(translate)

‎gdal_clip.py

+22
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,22 @@
1+
import os
2+
3+
# gdalwarp https://gdal.org/programs/gdalwarp.html
4+
def clip(shp, INPUT, OUTPUT):
5+
gdalwarp = 'gdalwarp -cutline %s -dstalpha %s %s' % (shp, INPUT, OUTPUT)
6+
os.system(gdalwarp)
7+
8+
def main():
9+
rootdir = r'H:\FINLAND_MODIS\Data\IM'
10+
output_path = rootdir + '\\clip'
11+
if not os.path.exists(output_path):
12+
os.makedirs(output_path)
13+
# shp = r'G:\MODIS_Lake_boundary\lake-polygons-PML-wkt_to_shp\GLWD00000677_BEAR-Lake.shp'
14+
shp = r'G:\MODIS_Lake_boundary\CCI_Lakes_fixed.shp'
15+
16+
for file in os.listdir(rootdir):
17+
if file.endswith(".tif"):
18+
print(os.path.join(rootdir, file))
19+
clip(shp, os.path.join(rootdir, file), os.path.join(output_path, file))
20+
21+
if __name__ == "__main__":
22+
main()

‎gdal_convert_format.py

+11
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,11 @@
1+
## Converting formats for geospatial rasters
2+
from osgeo import gdal
3+
4+
in_image = gdal.Open(r"C:\Users\Junqian Wang\Downloads\T10SFF_20180917T185019_B03.jp2")
5+
6+
driver = gdal.GetDriverByName("GTiff")
7+
8+
out_image = driver.CreateCopy(r"C:\Users\Junqian Wang\Downloads\test\T10SFF_20180917T185019_B03.tif", in_image, 0)
9+
10+
in_image = None
11+
out_image = None

‎gdal_read_write_geotiff/read_write_geotiff.py

+35
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,35 @@
1+
from osgeo import gdal
2+
3+
filename = r'D:\Python scripts\Threshold\subset_1_of_S1A_IW_GRDH_1SDV_20171205T165930_20171205T165959_019570_0213BE_1251_40.tif'
4+
ds = gdal.Open(filename, gdal.GA_ReadOnly)
5+
cols = ds.RasterXSize
6+
rows = ds.RasterYSize
7+
driver = ds.GetDriver().LongName
8+
geotransform = ds.GetGeoTransform()
9+
originX = geotransform[0] #top left x
10+
originY = geotransform[3] # top left y
11+
pixelWidth = geotransform[1] #w-e pixel resolution
12+
pixelHeight = -geotransform[5] #n-s pixel resolution
13+
14+
#read in the first two bands into arrays
15+
band1 = ds.GetRasterBand(1)
16+
band2 = ds.GetRasterBand(2)
17+
band1_array = band1.ReadAsArray()
18+
band2_array = band2.ReadAsArray()
19+
20+
#apply some band math
21+
result_array = band1_array - band2_array #difference of the two bands
22+
23+
##write the image
24+
outfile = r'D:\Python scripts\Threshold\output.tif'
25+
driver = gdal.GetDriverByName("GTiff")
26+
outdata = driver.Create(outfile, cols, rows, 1, gdal.GDT_Float32)
27+
outdata.SetGeoTransform(geotransform)
28+
outdata.SetProjection(ds.GetProjection())##sets same projection as input
29+
outband = outdata.GetRasterBand(1)
30+
outband.WriteArray(result_array)
31+
outband.FlushCache() ##saves to disk!!
32+
outdata = None
33+
outband = None
34+
ds = None
35+
band_theta = None

‎gdal_read_write_geotiff/subset_1_of_S1A_IW_GRDH_1SDV_20171205T165930_20171205T165959_019570_0213BE_1251_40.tif

22.1 MB
Binary file not shown.

‎gdal_subset.py

+49
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,49 @@
1+
import ogr
2+
from osgeo import gdal
3+
import os
4+
from pyproj import Proj, transform
5+
6+
def subset(INPUT, OUTPUT, xmin, xmax, ymin, ymax):
7+
translate = 'gdal_translate -projwin %s %s %s %s %s %s' %(xmin, ymax, xmax, ymin, INPUT, OUTPUT)
8+
os.system(translate)
9+
10+
def main():
11+
rootdir = r'H:\FINLAND_MODIS\Data\IM'
12+
output_path = rootdir + '\\Subset'
13+
if not os.path.exists(output_path):
14+
os.makedirs(output_path)
15+
##get lake boundary extent
16+
# shp = 'Barrow.shp'
17+
# shp = 'Kytalyk_SiberiaC2.shp'
18+
# shp = 'Lena_delta2c.shp'
19+
# shp = 'Mackenzie Delta_NWTcmodify.shp'
20+
# shp = 'Teshekpuk_Lake2c.shp'
21+
# shp = 'Yamal_SiberiaC2.shp'
22+
shp = 'H:\FINLAND_MODIS\shp\Finland_lakes.shp'
23+
shp_ds = ogr.Open(shp)
24+
source_layer = shp_ds.GetLayer()
25+
x_min, x_max, y_min, y_max = source_layer.GetExtent()
26+
print("lake boundary extent:", x_min, x_max, y_min, y_max)
27+
28+
# ## if shp and image are in different coordinate system
29+
# inProj = Proj(init='epsg:4230')
30+
# outProj = Proj(init='epsg:32662')
31+
# x_min, y_max = transform(inProj, outProj, x_min, y_max)
32+
# x_max, y_min = transform(inProj, outProj, x_max, y_min)
33+
34+
for file in os.listdir(rootdir):
35+
if file.endswith(".tif"):
36+
print (os.path.join(rootdir, file))
37+
# ds = gdal.Open(os.path.join(rootdir, file), gdal.GA_ReadOnly)
38+
# ulx, xres, xskew, uly, yskew, yres = ds.GetGeoTransform()
39+
# lrx = ulx + (ds.RasterXSize * xres)
40+
# lry = uly + (ds.RasterYSize * yres) # ulx, uly is the upper left corner, lrx, lry is the lower right corner
41+
# xmin = max(x_min, ulx)
42+
# xmax = min(x_max, lrx)
43+
# ymin = max(y_min, lry)
44+
# ymax = min(y_max, uly)
45+
# subset(os.path.join(rootdir, file), os.path.join(output_path, file), xmin, xmax, ymin, ymax)
46+
subset(os.path.join(rootdir, file), os.path.join(output_path, file), x_min, x_max, y_min, y_max)
47+
48+
if __name__== "__main__":
49+
main()

‎gdal_to_utm.py

+27
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,27 @@
1+
#GDAL Utilities is only installed for Python 3.4 (not anaconda) on this machine
2+
#reproject the images to UTM OR UTM to WGS84
3+
# output path has to be different from rootdir (bc output won't overwrite)
4+
import os
5+
import subprocess
6+
7+
rootdir = r'D:\GlobPermafrost\raw images\C12_1_Cape_Espenberg_Lowland\EW\DOY_THRESH'
8+
output_path = r'D:\GlobPermafrost\raw images\C12_1_Cape_Espenberg_Lowland\EW\DOY_THRESH\UTM'
9+
if not os.path.exists(output_path):
10+
os.makedirs(output_path)
11+
12+
##UTM Zone
13+
#C01_Barrow: zone 4
14+
#C03_Teshekpuk: zone 5
15+
#C04_Mackenzie: zone 8
16+
#C06_Kytalyk: zone 55
17+
#C07_LenaDelta: zone 52
18+
#C08_Yamal: zone 42
19+
#CS12_1: Cape Espenberg Lowland (Alaska, USA):3
20+
#CS12_2: Central Seward Peninsula (Alaska, USA):3
21+
for file in os.listdir(rootdir):
22+
if file.endswith(".tif"):
23+
print (os.path.join(rootdir, file))
24+
subprocess.call(['gdalwarp','-t_srs', "+proj=utm +zone=3 +datum=WGS84 +units=m +no_defs", os.path.join(rootdir, file), os.path.join(output_path, file)], shell=True)
25+
# subprocess.call(['gdalwarp', '-of', 'GTiff', '-t_srs', 'EPSG:32604', os.path.join(rootdir, file),os.path.join(output_path, file)], shell=True)
26+
# subprocess.call(['gdalwarp', os.path.join(rootdir, file), os.path.join(output_path, file),'-t_srs', "+proj=longlat +zone=4 +ellps=WGS84"], shell=True)
27+

‎georeference.py

+110
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,110 @@
1+
##georeference output from MAGIC
2+
#IRGS result colour: grounded ice (0,0,254), floating ice (0,254,254) (MAGIC automatically saves 255 to 254)
3+
#IRGS result (.bmp) should have the same name as the SAR image (.tif).
4+
#IRGS results and SAR images are in the same folder
5+
#output in a separte folder
6+
from osgeo import gdal
7+
import numpy as np
8+
import os
9+
from PIL import Image
10+
import subprocess
11+
12+
def georeference(filename,output_path,code,percent):
13+
ds = gdal.Open(filename, gdal.GA_ReadOnly)
14+
cols = ds.RasterXSize
15+
rows = ds.RasterYSize
16+
17+
bmp = filename[:-4] + '.bmp'
18+
im = Image.open(bmp)
19+
p = np.array(im)
20+
data = np.zeros((rows, cols), dtype=np.float)
21+
landmask = filename[:-4] + '_landmask.bmp'
22+
mask_array = gdal.Open(landmask).ReadAsArray()
23+
24+
for i in range(0, rows):
25+
for j in range(0, cols):
26+
if mask_array[i,j] == 0:
27+
data[i,j] = np.NAN
28+
elif np.all(p[i,j] == [0,254,254]):
29+
data[i,j] = 2
30+
elif np.all(p[i,j] == [0,0,254]):
31+
data[i,j] = 1
32+
33+
grounded = np.count_nonzero(data == 1)
34+
floating = np.count_nonzero(data == 2)
35+
percent_g = round(grounded/(grounded+floating), 4)
36+
percent_f = round(floating/(grounded+floating), 4)
37+
print ("grounded ice%: ", percent_g)
38+
print ("floating ice%: ", percent_f)
39+
40+
##write the image
41+
##DATE
42+
d_start = filename.find('_EW_GRDM') + 14
43+
d_end = filename.find('T', d_start)
44+
date = filename[d_start:d_end]
45+
#percent = output_path + "\\Percentage.txt"
46+
with open(percent, "a") as file:
47+
#file.write('IRGSEM\n'+'Date\tGrounded ice\tFloating ice'+'\n')
48+
file.write(date +'\t'+str(percent_g)+'\t'+str(percent_f)+'\n')
49+
outfile = output_path + '\\' + 'H2O.CGI.IRGSEM.SENT1.V01.' + date + '.' + code + '.tif'
50+
driver = gdal.GetDriverByName("GTiff")
51+
outdata = driver.Create(outfile, cols, rows, 1, gdal.GDT_Byte)
52+
outdata.SetGeoTransform(ds.GetGeoTransform())
53+
outdata.SetProjection(ds.GetProjection())##sets same projection as input
54+
outband = outdata.GetRasterBand(1)
55+
ct = gdal.ColorTable()
56+
#ct.SetColorEntry(0, (255,255,255,255)) #white background for NaN
57+
ct.SetColorEntry(0, (0,0,0,255))
58+
ct.SetColorEntry(1, (0,0,255,255)) #grounded ice color(0,0,255), used to be (0,0,102)
59+
ct.SetColorEntry(2, (0,255,255,255)) #floating ice color(0,255,255)
60+
outband.SetColorTable(ct)
61+
outband.WriteArray(data)
62+
#outband.SetNoDataValue(0)
63+
outband.FlushCache() ##saves to disk!!
64+
outdata = None
65+
outband = None
66+
ds = None
67+
68+
69+
def main():
70+
rootdir = r'D:\GlobPermafrost\raw images\C12_Central_Seward_Peninsula\EW\IRGSEM'
71+
output_georef = rootdir + '\\' + 'georef'
72+
if not os.path.exists(output_georef):
73+
os.makedirs(output_georef)
74+
code = 'CS12_2'
75+
#C01: Barrow (Alaska, USA)
76+
#C03: Teshekpuk (Alaska, USA)
77+
#C04: Mackenzie Delta (Canada)
78+
#C06: Kytalyk (Russia)
79+
#C07: Lena Delta, Russia
80+
#C08: Yamal (Russia)
81+
#CS12_1: Cape Espenberg Lowland (Alaska, USA)
82+
#CS12_2: Central Seward Peninsula (Alaska, USA)
83+
percent = output_georef + "\\Percentage.txt"
84+
with open(percent, "a") as file:
85+
file.write('IRGSEM\n'+'Date\tGrounded ice\tFloating ice'+'\n')
86+
for file in os.listdir(rootdir):
87+
if file.endswith(".tif"):
88+
print (os.path.join(rootdir, file))
89+
georeference(os.path.join(rootdir, file), output_georef, code, percent)
90+
output_path = output_georef + '\\' + 'UTM'
91+
if not os.path.exists(output_path):
92+
os.makedirs(output_path)
93+
##UTM Zone
94+
# C01_Barrow: zone 4
95+
# C03_Teshekpuk: zone 5
96+
# C04_Mackenzie: zone 8
97+
# C06_Kytalyk: zone 55
98+
# C07_LenaDelta: zone 52
99+
# C08_Yamal: zone 42
100+
# CS12_1: Cape Espenberg Lowland (Alaska, USA):3
101+
# CS12_2: Central Seward Peninsula (Alaska, USA):3
102+
for file in os.listdir(output_georef):
103+
if file.endswith(".tif"):
104+
print(os.path.join(output_georef, file))
105+
subprocess.call(
106+
['gdalwarp', '-t_srs', "+proj=utm +zone=3 +datum=WGS84 +units=m +no_defs", os.path.join(output_georef, file),
107+
os.path.join(output_path, file)], shell=True)
108+
109+
if __name__== "__main__":
110+
main()

‎georeference_cci.py

+95
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,95 @@
1+
##georeference output from MAGIC
2+
#IRGS result colour: water (0,0,254), Ice (254,254,0) (MAGIC automatically saves 255 to 254)
3+
#IRGS result (.bmp) should have the same name as the SAR image (.tif).
4+
#IRGS results and SAR images are in the same folder
5+
#output in a separte folder
6+
from osgeo import gdal
7+
import numpy as np
8+
import os
9+
from PIL import Image
10+
import ntpath
11+
12+
13+
Image.MAX_IMAGE_PIXELS = 1000000000
14+
def georeference(filename,output_path):
15+
ds = gdal.Open(filename, gdal.GA_ReadOnly)
16+
cols = ds.RasterXSize
17+
rows = ds.RasterYSize
18+
19+
bmp = filename[:-4] + '.bmp'
20+
im = Image.open(bmp)
21+
p = np.array(im)
22+
data = np.zeros((rows, cols), dtype=np.float)
23+
landmask = filename[:-4] + '_landmask.bmp'
24+
mask_array = gdal.Open(landmask).ReadAsArray()
25+
26+
# for i in range(0, rows):
27+
# for j in range(0, cols):
28+
# if mask_array[i,j] == 0:
29+
# data[i,j] = np.NAN
30+
# elif np.all(p[i,j] == [254,254,0]):
31+
# data[i,j] = 2 #ice
32+
# elif np.all(p[i,j] == [0,0,254]):
33+
# data[i,j] = 1 #water
34+
for i in range(0, rows):
35+
for j in range(0, cols):
36+
if mask_array[i,j] == 0:
37+
data[i,j] = np.NAN
38+
elif np.all(p[i,j] == [254,254,0]):
39+
data[i,j] = 1
40+
elif np.all(p[i,j] == [0,0,254]):
41+
data[i,j] = 2
42+
elif np.all(p[i,j] == [254,152,0]):
43+
data[i,j] = 3
44+
elif np.all(p[i,j] == [0,254,254]):
45+
data[i,j] = 4
46+
elif np.all(p[i,j] == [127,0,127]):
47+
data[i,j] = 5
48+
elif np.all(p[i,j] == [0,127,0]):
49+
data[i,j] = 6
50+
51+
# water = np.count_nonzero(data == 1)
52+
# ice = np.count_nonzero(data == 2)
53+
# percent_w = round(water/(water+ice), 4)
54+
# percent_i = round(ice/(water+ice), 4)
55+
# print ("water%: ", percent_w)
56+
# print ("ice%: ", percent_i)
57+
58+
##write the image
59+
# with open(percent, "a") as file:
60+
# file.write(date +'\t'+str(percent_w)+'\t'+str(percent_i)+'\n')
61+
outfile = output_path + '\\' + ntpath.basename(os.path.normpath(filename))[:-4] + '_IRGS.tif'
62+
driver = gdal.GetDriverByName("GTiff")
63+
outdata = driver.Create(outfile, cols, rows, 1, gdal.GDT_Byte)
64+
outdata.SetGeoTransform(ds.GetGeoTransform())
65+
outdata.SetProjection(ds.GetProjection())##sets same projection as input
66+
outband = outdata.GetRasterBand(1)
67+
# ct = gdal.ColorTable()
68+
# #ct.SetColorEntry(0, (255,255,255,255)) #white background for NaN
69+
# ct.SetColorEntry(0, (0,0,0,255))
70+
# ct.SetColorEntry(1, (0,112,255,255)) #water color(0,112,255)
71+
# ct.SetColorEntry(2, (255,255,0,255)) #ice color(255,255,0)
72+
# outband.SetColorTable(ct)
73+
outband.WriteArray(data)
74+
#outband.SetNoDataValue(0)
75+
outband.FlushCache() ##saves to disk!!
76+
outdata = None
77+
outband = None
78+
ds = None
79+
80+
81+
def main():
82+
rootdir = r'C:\Users\Junqian Wang\Downloads\New folder'
83+
output_path = r'C:\Users\Junqian Wang\Downloads\New folder\IRGS'
84+
if not os.path.exists(output_path):
85+
os.makedirs(output_path)
86+
# percent = output_path + "\\Percentage.txt"
87+
# with open(percent, "a") as file:
88+
# file.write('IRGSEM\n'+'Date\tWATER\tICE'+'\n')
89+
for file in os.listdir(rootdir):
90+
if file.endswith(".tif"):
91+
print (os.path.join(rootdir, file))
92+
georeference(os.path.join(rootdir, file), output_path)
93+
94+
if __name__== "__main__":
95+
main()

0 commit comments

Comments
 (0)