-
Notifications
You must be signed in to change notification settings - Fork 12
/
get_crs_urn.py
69 lines (47 loc) · 1.45 KB
/
get_crs_urn.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
#!/usr/bin/python3
import sys
import re
import json
from os import path
try:
from osgeo import ogr, osr, gdal
except:
sys.exit('ERROR: cannot find GDAL/OGR modules')
GEO_TIFF_CRS_PATTERN = re.compile(
'(?:WGS 84 / (UTM|Pseudo-Mercator)(?: zone ([0-9]{1,2})([SN])))'
)
def get_crs_urn(raster_file):
crs_tag = get_tif_crs_tag(raster_file)
epsg_code = get_epsg_code(crs_tag)
crs_urn = create_crs_urn(epsg_code)
return crs_urn
def get_tif_crs_tag(raster_file):
ds = gdal.Open(raster_file)
prj = ds.GetProjection()
srs = osr.SpatialReference(wkt=prj)
return srs.GetAttrValue('projcs')
def get_epsg_code(crs_tag):
m = GEO_TIFF_CRS_PATTERN.match(crs_tag)
if m:
if m.group(1) is None:
return 4326
if m.group(1) == 'Pseudo-Mercator':
return 3857
if m.group(1) == 'UTM':
zone = int(m.group(2))
if m.group(3) == 'N':
return 32600 + zone
elif m.group(3) == 'S':
return 32700 + zone
raise ValueError(f'No EPSG code in tag {crs_tag}')
def create_crs_urn(epsg_code):
return f'urn:ogc:def:crs:EPSG::{epsg_code}'
if __name__ == '__main__':
if len(sys.argv) != 2:
print('Provide a path to GeoTIFF as the first argument')
sys.exit(1)
file = sys.argv[1]
if not path.exists(file):
print(f'File "{file}" does not exist')
sys.exit(1)
print(get_crs_urn(file))