-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathCreate_NDWI_Coastline.py
150 lines (131 loc) · 5.51 KB
/
Create_NDWI_Coastline.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
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
import numpy as np
import scipy.ndimage
import geopandas as gpd
from osgeo import gdal,gdalconst,osr
import argparse
import subprocess
from shapely.geometry import Polygon, MultiPolygon
def load_image(img_path):
'''
Loads a geotiff file as a numpy array.
'''
img = gdal.Open(img_path,gdalconst.GA_ReadOnly)
img_array = np.array(img.GetRasterBand(1).ReadAsArray())
return img_array
def flood_fill(img_binary,erosion_struc=np.ones((3,3))):
'''
Performs a flood-fill operation on a binary image.
'''
img_binary_eroded = scipy.ndimage.binary_dilation(img_binary)
img_binary_floodfill = scipy.ndimage.binary_erosion(img_binary_eroded,structure=erosion_struc)
return img_binary_floodfill
def write_file(raster,file_name,src):
'''
Given a numpy array and a file name, writes the array to a geotiff file, using the information in src.
'''
driver = gdal.GetDriverByName('GTiff')
out_raster = driver.Create(file_name,raster.shape[1],raster.shape[0],1,gdal.GDT_Byte)
out_raster.SetGeoTransform(src.GetGeoTransform())
out_raster.SetProjection(src.GetProjection())
out_raster.GetRasterBand(1).WriteArray(raster)
out_raster.FlushCache()
out_raster = None
def get_lonlat_bounds_gdf(gdf):
'''
Returns the lon/lat boundarys of an entire GeoDataFrame.
'''
lon_min = np.min(gdf.bounds.minx)
lon_max = np.max(gdf.bounds.maxx)
lat_min = np.min(gdf.bounds.miny)
lat_max = np.max(gdf.bounds.maxy)
return lon_min,lon_max,lat_min,lat_max
def lonlat2epsg(lon,lat):
'''
Finds the EPSG code for a given lon/lat coordinate.
'''
if lat >= 0:
NS_code = '6'
elif lat < 0:
NS_code = '7'
EW_code = f'{int(np.floor(lon/6.0))+31:02d}'
epsg_code = f'32{NS_code}{EW_code}'
return epsg_code
def remove_interior_holes_polygon(geom,threshold):
'''
Removes interior holes smaller than a given threshold from the given geometry.
'''
if geom.geom_type == 'Polygon':
interior_list = []
for interior in geom.interiors:
p = Polygon(interior)
if p.area > threshold:
interior_list.append(interior)
new_geom = Polygon(geom.exterior.coords,holes=interior_list)
elif geom.geom_type == 'MultiPolygon':
polygon_list = []
for polygon in geom.geoms:
interior_list = []
for interior in geom.interiors:
p = Polygon(interior)
if p.area > threshold:
interior_list.append(interior)
polygon_list.append(Polygon(polygon.exterior.coords, holes=interior_list))
new_geom = MultiPolygon(geom.exterior.coords,holes=interior_list)
else:
new_geom = geom
return new_geom
def simplify_polygon(geom,threshold):
'''
Simplifies the polygon
'''
new_geom = geom.simplify(threshold)
return new_geom
def main():
'''
Write by Eduard Heijkoop, University of Colorado
Update June 2022: Created script
This script creates a coastlines shapefile from an NDWI image (e.g. from Sentinel-2 or Landsat-8/-9).
'''
parser = argparse.ArgumentParser()
parser.add_argument('--input_file',help='Path to input NDWI file.')
args = parser.parse_args()
input_file = args.input_file
NDWI_THRESHOLD = 0.0
SHP_POLYGON_THRESHOLD = 0.005
INTERIOR_THRESHOLD = 50000
SIMPLIFY_RADIUS = 20
output_file_floodfill = input_file.replace('.tif','_floodfill.tif')
output_file_floodfill_nodata = output_file_floodfill.replace('.tif','_nodata_0.tif')
output_file_floodfill_shp = output_file_floodfill.replace('.tif','.shp')
output_file_shp = input_file.replace('.tif','.shp')
output_file_shp_simplified = output_file_shp.replace('.shp','_simplified.shp')
src = gdal.Open(input_file,gdalconst.GA_ReadOnly)
ndwi = np.array(src.GetRasterBand(1).ReadAsArray())
ndwi_lt_0p0 = ndwi < NDWI_THRESHOLD
ndwi_lt_0p0 = ndwi_lt_0p0.astype(int)
ndwi_lt_0p0_floodfill = flood_fill(ndwi_lt_0p0)
write_file(ndwi_lt_0p0_floodfill,output_file_floodfill,src)
nodata_command = f'gdal_translate -q -a_nodata 0 -of GTiff -co "COMPRESS=LZW" {output_file_floodfill} {output_file_floodfill_nodata}'
polygonize_command = f'gdal_polygonize.py -q {output_file_floodfill_nodata} -f "ESRI Shapefile" {output_file_floodfill_shp}'
subprocess.run(nodata_command,shell=True)
subprocess.run(polygonize_command,shell=True)
subprocess.run(f'rm {output_file_floodfill}',shell=True)
subprocess.run(f'rm {output_file_floodfill_nodata}',shell=True)
gdf = gpd.read_file(output_file_floodfill_shp)
lon_min,lon_max,lat_min,lat_max = get_lonlat_bounds_gdf(gdf)
lon_center = np.mean((lon_min,lon_max))
lat_center = np.mean((lat_min,lat_max))
epsg_code = lonlat2epsg(lon_center,lat_center)
gdf = gdf.to_crs(f'EPSG:{epsg_code}')
idx_area = np.argwhere(np.asarray(gdf.area) > SHP_POLYGON_THRESHOLD * np.sum(gdf.area)).squeeze()
gdf = gdf.iloc[idx_area].reset_index(drop=True)
gdf['geometry'] = gdf.apply(lambda x : remove_interior_holes_polygon(x.geometry,INTERIOR_THRESHOLD),axis=1)
gdf_simplified = gdf
gdf_simplified['geometry'] = gdf_simplified.apply(lambda x : simplify_polygon(x.geometry,SIMPLIFY_RADIUS),axis=1)
gdf = gdf.to_crs('EPSG:4326')
gdf.to_file(output_file_shp)
gdf_simplified = gdf_simplified.to_crs('EPSG:4326')
gdf_simplified.to_file(output_file_shp_simplified)
subprocess.run(f'rm {output_file_floodfill_shp.replace(".shp",".*")}',shell=True)
if __name__ == '__main__':
main()