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
|
>>> from liblas import srs
>>> from liblas import point
>>> from liblas import header
>>> import liblas
>>> s = srs.SRS()
>>> s.proj4
''
>>> s.set_proj4(b'+proj=utm +zone=15 +datum=WGS84 +units=m +no_defs')
True
>>> print s.proj4
+proj=utm +zone=15 +datum=WGS84 +units=m +no_defs
>>> s.proj4 == '+proj=utm +zone=15 +datum=WGS84 +units=m +no_defs '
True
>>> s = srs.SRS()
>>> s.set_userinput(b'EPSG:4326')
True
>>> s.proj4 == '+proj=longlat +datum=WGS84 +no_defs '
True
>>> from liblas import file
>>> f = file.File('../test/data/1.2_3.las',mode='r')
>>> s = f.header.srs
>>> s.wkt == b"""PROJCS["NAD83 / UTM zone 15N",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.2572221010002,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4269"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-93],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","26915"]]"""
True
>>> s2 = srs.SRS()
>>> s2.wkt = b"""GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]"""
>>> p = f.read(0)
>>> p.x
470692.44
>>> f.set_srs(s2)
True
>>> p = f.read(0)
>>> s2.vlr_count()
4
>>> s2.GetVLR(0).recordlength
64
>>> int(round(p.x))
-93
>>> int(round(p.y))
42
>>> del f
# -93.3515625902 41.5771483954
>>> def new_offset(old_scale, new_scale, old_offset, x):
... return (new_scale*(x - old_offset) - old_scale*x)/(-1.0*old_scale)
>>> utm_wkt = b"""PROJCS["NAD83 / UTM zone 15N",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.2572221010002,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4269"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-93],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","26915"]]"""
>>> dd_wkt = b"""GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]"""
>>> s_dd = srs.SRS()
>>> s_dd.wkt = dd_wkt
>>> s_utm = srs.SRS()
>>> s_utm.wkt = utm_wkt
>>> f = file.File('../test/data/1.2_3.las',mode='r')
>>> dd_header = liblas.header.Header(handle=f.header.handle,copy=True)
>>> utm_header = liblas.header.Header(handle=f.header.handle,copy=True)
>>> del f
>>> utm_header.offset
[0.0, 0.0, 0.0]
>>> utm_header.scale
[0.01, 0.01, 0.01]
# >>> new_offset(0.01, 0.0001, 0.0, 470692.44)
# >>> utm_header.offset = [offset+1.0/0.000001 for offset in utm_header.offset]
# >>> utm_header.offset
[1000000.0, 1000000.0, 1000000.0]
# >>> utm_header.scale = [0.000001,0.000001,0.000001]
>>> utm_header.srs = s_utm
# >>> dd_header.scale = [0.000001,0.000001,0.000001]
>>> dd_header.srs = s_dd
>>> f = file.File('../test/data/1.2_3.las',mode='r', header = utm_header)
>>> int(f.header.data_offset)
438
>>> f.header.scale
[0.01, 0.01, 0.01]
>>> p = f.read(0)
>>> origx, origy = p.x, p.y
>>> "{:.2f}, {:.2f}".format(origx, origy)
'470692.44, 4602888.90'
>>> f.set_srs(s_dd)
True
>>> p = f.read(0)
We only get truncated values because our header scale
values are 0.01
>>> "{:.2f}, {:.2f}".format(p.x, p.y)
'-93.35, 41.58'
#real values
# (-93.351562590199833, 41.577148395415108)
>>> f_project = file.File('junk_srs_project.las',mode='w',header=dd_header)
>>> p.header = dd_header
>>> "{:.2f}, {:.2f}".format(p.x, p.y)
'-93.35, 41.58'
>>> dd_header.srs.proj4
'+proj=longlat +datum=WGS84 +no_defs '
>>> f_project.write(p);f_project.write(p);f_project.write(p)
>>> f_project.close()
>>> del f_project
>>> f3 = file.File('junk_srs_project.las')
>>> int(f3.header.data_offset)
789
>>> s_utm = srs.SRS()
>>> s_utm.wkt = utm_wkt
>>> p3 = f3.read(1)
>>> int(round(p3.x)), int(round(p3.y))
(-93, 42)
>>> p3 = f3.read(0)
>>> int(round(p3.x)), int(round(p3.y))
(-93, 42)
>>> import os
>>> os.remove('junk_srs_project.las')
>>> f = file.File('../test/data/srs_vertcs.las',mode='r')
>>> s = f.header.srs
>>> str(s.get_wkt_compoundok().decode())
'COMPD_CS["unknown",PROJCS["WGS 84 / UTM zone 17N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-81],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32617"]],VERT_CS["NAVD88 height",VERT_DATUM["North American Vertical Datum 1988",2005,AUTHORITY["EPSG","5103"],EXTENSION["PROJ4_GRIDS","g2012a_conus.gtx,g2012a_alaska.gtx,g2012a_guam.gtx,g2012a_hawaii.gtx,g2012a_puertorico.gtx,g2012a_samoa.gtx"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Up",UP],AUTHORITY["EPSG","5703"]]]'
>>> s2 = srs.SRS()
>>> s2.wkt = b"""PROJCS["WGS 84 / UTM zone 17N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-81],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32617"]]"""
>>> s2.set_verticalcs( 5703, b'abc', 5103, 9001 )
True
>>> str(s2.get_wkt_compoundok().decode())
'COMPD_CS["unknown",PROJCS["WGS 84 / UTM zone 17N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-81],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32617"]],VERT_CS["NAVD88 height",VERT_DATUM["North American Vertical Datum 1988",2005,AUTHORITY["EPSG","5103"],EXTENSION["PROJ4_GRIDS","g2012a_conus.gtx,g2012a_alaska.gtx,g2012a_guam.gtx,g2012a_hawaii.gtx,g2012a_puertorico.gtx,g2012a_samoa.gtx"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Up",UP],AUTHORITY["EPSG","5703"]]]'
|
