Forum QGIS

Pełna wersja: NMT z CODGiK. Jak zamienić punkty na linie
Aktualnie przeglądasz uproszczoną wersję forum. Kliknij tutaj, by zobaczyć wersję z pełnym formatowaniem.
Stron: 1 2
Witam,

Proszę o pomoc jak część danych zakupionych w CODGiK (.asc) zamienić na linie. Rozdział jest za pomocą znaczników START/END. Udaje mi się jedynie zaimportować te dane jako punkty, ale łączenie zawsze źle wychodzi - nie potrafię otrzymać z jednego pliku szeregu linii złożonych z punktów pomiędzy START/END.
Problem ten był sygnalizowany na tym forum, ale nie znalazłem rozwiązania.

Struktura pliku wygląda tak:
Start
592793.63 778717.72 136.55
592788.07 778711.53 136.54
592782.29 778705.60 136.54
592773.54 778702.50 136.53
592768.93 778702.05 136.53
592762.83 778703.01 136.53
592758.77 778705.79 136.52
592756.16 778710.65 136.52
592755.39 778719.99 136.52
592755.94 778725.27 136.52
592757.95 778728.65 136.52
592763.65 778731.20 136.53
592767.54 778731.58 136.53
592791.59 778726.25 137.11
592794.58 778723.81 137.00
592795.40 778721.23 137.00
592795.03 778719.20 137.00
592793.63 778717.72 136.55
End
Start
593235.02 778420.16 133.73
593238.89 778430.51 132.83
593239.51 778435.76 132.39
593239.56 778438.78 132.39
593239.01 778441.71 132.39
593219.02 778467.30 132.37
593201.31 778488.96 132.58
593199.93 778488.63 132.81
593198.20 778487.88 132.80
593194.27 778484.85 132.80
593182.30 778470.43 131.89
593178.97 778467.24 131.89
593174.93 778464.22 131.89
593169.50 778460.31 131.89
End
Start
593164.49 778462.76 132.56
593169.36 778463.87 132.78
593173.89 778466.86 132.79
593177.60 778470.39 132.79
593194.41 778489.11 133.70
593196.80 778490.66 133.70
593199.47 778491.85 134.15
593202.69 778491.18 134.15
593215.43 778476.92 134.61
593233.22 778455.09 134.63
593240.61 778445.36 134.63
593242.09 778441.77 134.63
593244.78 778434.25 134.64
End
Start
596810.72 777312.02 142.63
596743.74 777305.73 141.74 (...)

Będę bardzo wdzięczny za pomoc.

Pozdrawiam,
Emeryt
Zakładając, że ten tekstowy plik wygląda jak wkleiłeś...

Trochę brzydki skryp, ale chyba robiący to co potrzebujesz.
Odpal konsole Pythona, edytor i wklej poniższy skrypt.
Ustaw ścieżkę do swojego pliku  (u mnie to:     'h:/python/asci/linie.txt') i odpal skrypt.

Ustawiłem domyślny układ współrzędnych na 2176   (2000/5 - jeśli dane są w innym to zmień na właściwy numer)

Kod:
from osgeo import ogr
from itertools import takewhile
from qgis.core import *

#sciezka do pliku z tekstem
plik = 'h:/python/asci/linie.txt'

warstwa_liniowa = QgsVectorLayer("LineString?crs=epsg:2176", 'test', "memory")
warstwa_liniowa.startEditing()
feat = QgsFeature(warstwa_liniowa.pendingFields())

linie = []
with open(plik) as f:
       array = []
       for line in f:
           if line.startswith('Start'):
               data = takewhile(lambda x: not x.startswith("End"),f)
               linia = ogr.Geometry(ogr.wkbLineString)
               for wiersz in data:
                   wsp = wiersz.split(' ')
                   linia.AddPoint(float(wsp[1]),float(wsp[0]) )
               gotowa_geom = linia.ExportToWkt()
               geom = QgsGeometry.fromWkt(gotowa_geom)
               feat=QgsFeature(warstwa_liniowa.pendingFields())
               feat.setGeometry(geom)
               warstwa_liniowa.addFeatures([feat])
warstwa_liniowa.commitChanges()
warstwa_liniowa.updateExtents()
QgsMapLayerRegistry.instance().addMapLayer(warstwa_liniowa, True)
Dziękuję xmaziax za poświęcony czas. Jednak z jakiegoś powodu skrypt nie działa.
Otrzymuję błąd jak poniżej:

Cytat:from osgeo import ogr
from itertools import takewhile
from qgis.core import *
#sciezka do pliku z tekstem
plik = 'D:\Testy\NMT-N34107Cb1_c.asc'
warstwa_liniowa = QgsVectorLayer("LineString?crs=epsg:2180", 'test', "memory")
warstwa_liniowa.startEditing()
True
feat = QgsFeature(warstwa_liniowa.pendingFields())
linie = []
with open(plik) as f:
      array = []
      for line in f:
          if line.startswith('Start'):
              data = takewhile(lambda x: not x.startswith("End"),f)
              linia = ogr.Geometry(ogr.wkbLineString)
              for wiersz in data:
                  wsp = wiersz.split(' ')
                  linia.AddPoint(float(wsp[1]),float(wsp[0]) )
              gotowa_geom = linia.ExportToWkt()
              geom = QgsGeometry.fromWkt(gotowa_geom)
              feat=QgsFeature(warstwa_liniowa.pendingFields())
              feat.setGeometry(geom)
              warstwa_liniowa.addFeatures([feat])
warstwa_liniowa.commitChanges()
 File "<input>", line 15
   warstwa_liniowa.commitChanges()
                 ^
SyntaxError: invalid syntax
Jeżeli chciałbyś sprawdzić jak dokładnie wygląda plik i jego struktura, a to może decydować o błędzie, to wrzuciłem tutaj przykładowy https://ufile.io/c5zyl
Dziękuję za pomoc.
Emeryt
Zmień ukośniki w ścieżce: z  'D:\Testy\NMT-N34107Cb1_c.asc' na  'D:/Testy/NMT-N34107Cb1_c.asc'

U mnie to zadziałało, w linku masz plik z zapisem wygenerowanej warstwy, który zapisałem do shp.
https://www.dropbox.com/s/dttx1hsqulnt5mw/NMT.zip?dl=0

[img]data:application/x-zip-compressed;base64,UEsDBBQAAAAIAO5IhkyO288lcAAAACICAAAHAAAATk1ULmRiZm3NOwoCQRBF0QbFRN2AkaFh/7snFEQwcTmuXwzEg8wLigMV3M1ru7uEEK7hENZ2f9y+fH7O/u9/PP8WccIZF1xxwx0PPPFCKuKEMy644oY7HnhiujnihDMuuOKGOx54Yrol4oQzLrjihjseeOLl9AZQSwMEFAAAAAgA7kiGTC3vY0IHAQAAkAEAAAcAAABOTVQucHJqbZDRa4MwEMb/lzyHksRmmsfShs6BVdQ+SQnBnjbgIsRsg/31O9nDJvQeQvjudx/fXVWXb8emI7qtm0yZap6sv5tjIxhj5nv2YCShZ12eVwgfs4KGq0wRejq016Ijp/9aU73qusxPCP9qjNCXJM14klKhsp2QqRCCM3670arOC10gGQD8l+sfhLIbvV7yFk1hRBWFHePpXiZCCcmV2uNH4ugaWx/bvLx0pA3WL58QFjAFhN7GORAkDvWh0K2uOzLZ6OLHHcw8mDm40Xn03RA9+BjsZN4huLuz2OdyAyy9ncAMtkfzNZTCEskGGeyECcAu0fmRUCnZWk8QP4f4WJm/bQuIgL54lB9QSwMEFAAAAAgA7kiGTLyHMh5hAQAAbwIAAAcAAABOTVQucXBqbZBfb4IwFMXf9ylInxvXFhH6aJQgJggpNctCDGnwqiRYllK3bJ9+Zcv+GL19as/v9N57CpGvF2WFYinKiHuPXtF3Su+9RckIId5Hr8ELEE7iPPnDEF7O5TZz14vpX0DpWoIxMFjTqq4WcAADuoG6fB8snGvKR0tZrGKRp8sKJaL03BtBeOaHEfVDzHg0YUHIGKOE4vlWrnKRymfXoCgThFFIKEe7HZb5U1JG04rgf2d3xzBjQTQaCpFmceZaGgD91jYn5By3eMQJHfHtJpUV2sPR4Y6cEBpOA59xFlDOp75/x8opY6P1Vpn+zpCv44VM802FpFF6eAUzQJ2BaZTtDXLEXMyzWMaiQp2yrb3soe4PdW/aY6vdHFdEA9oaF/MZTLtvldNpcAUMjeqgPqjGfT4uwV0x/wo5qM5NAGqwrT4iHARkrDuI7o09fTHkJ50zWAOu6b0syHeMtwqj4cwpD59QSwMEFAAAAAgA7kiGTHKu7vNhCgAAlBsAAAcAAABOTVQuc2hwjVgFiFRRFB272xlHV2fsVTHHblexW8HCBBUTAwW7W7FrbVEMbGeNWWt1wNq1u7E7UbFA733nvv/ffhVddvjcefecd/u9Py5XZFrXH/8ynH6WzOWif9eZhISEbNsio7Jvr9/toDt/VPhTkY6rEiKjXrXcv6p2xfxRrn//JaFPLeZqcH9E9tQ7IqPmHPtUpHwhhR1Xl+RB1dNeG0AyK6YQUEdS8lprtm7aa+vedSR5WXR0dDmDS3N/LtIxUwOS/4BPCjtsn0Zm334iS/78Ua32+ybnJZmeKfz5TTts3fUDq6fNQWtEVyCX1rWxar+cJJMZJXMZ3E48USbTdtB3M1qFYfO3iPxRZcuWDY4luWfPnk+a5TbssHTttckvW1bqRDLzp8+tsTYX6w8LK7uGVqZ1J54ok2s7tte//+3osciooh3n1L5OaxzjOyTz84Nph62r+G7Tmp8MuaF1bazKxxtZP2FwO/FMre0g3Qtnj6na8p0TH58fQ2yfOuzITkQPjTXBcl3kvEwypaXv5d+5OD/bH5LsxBNlSm0H+5AtTvmwcIhP5bBJjThwRfsMOyxdVWdtx/vQD8W0roVFPcwimeNQ1OD+Az4VfVokEx8LH41UnMfzKB/2Dz8KbBqxI5XYwfnvbugKVu33gWTev5qFtbhUH93Kg/rpT3JP+rsp+sSnnmdIdvDzMzV9SrOdHPPqe9TaOC9xq9mwBzWRU+xMxkBLF/vmhP8Do/6O5e/T6LxQqjZFx9jYsvTXW3QHO/Ji8ti66MOGPtRhXIyNBTfmSR2SnXZyqmmE5BV/o9/uRu9N8mEmpohB7y3www49GLGmdVVdvvwcxPyYLfl/QTLX62LZ9wbJFJZD60lWczcIuw/6kJcjJHNdHfNh/uwOqr3nJEh9bSWZn+dJ5u83AD/jog9zcWUQc+AyyTwf54t8BXEZMTmIProi/KNJZr1LYv9Q2D/5Any/3Z9k7r8zot9T8KdE7hxEHZG9Kr7tSeZ6OijrrUnmOMWQTH5Xby7+bJH+aYL4vFwreWok+y9FvKJJVnNmng993ljiM0tyUxl8mWaQzPVfSeIzV/hIVnN7GeJ3rEoQ/btB6rNGEHWxG/m4XyeI8/CQyI0kXyfAn72V2HsO8c/UUfy/JnXYA3yb7oj/A4TvAfIzdLjUw2PE49BEyf9Tyc8iwT/zYb4uEn9fok7H9Qli7mf2g78m+O+X8yM+bpHbkiy1rPYZAv0UqYLox3l+5Le48NG68qsF6qFtGz/8H0Uyxynghz0jSOZ9qoI/53Dxr7Yf5+EwyUcTwQ8XmfhUXYyU+HUimfcZE8Tc6uFHfsZLvvuJvVNRP7fJPrXPLJI5TiOxX5P5kr/xfsRzGda/TZX4rAV/9ll+nJNbguj7OX70W4zUw3w/6vdI0Op19f1JsXe+xOuC4HX8zkl9TPCr/KTjK5E+76ZK74z05VdnUyGZI5NkjuRx4c/Bo+w/jri/m4i87j8StLAqX/ukLnidOXZK3VMcFG6j5G0M7Oy4Dn2wfaTUyQqZC0OhP3CJ8A1CXOovkL7pJ37Olrh2J5m/nyV93kXiNEP6oL3keZr0eWuSVSwkb4396Jtpsn8d+PNpmsS1huw/E3lZVUHyNIdkvp8VR53UXgB7EwpLXlcjftUjxN59Ulcpxd+7UudPfOB/ILl5QTLz/pS5c0nmlCcmUsU9zsodzs89ImeLsXOb+BxQ+UivzzWerdsPoWf6FEANpzsMW9YV+P3es9fQBRb1N0KwyTTW5lJ9M59kJ54v+9oO7tWLBzFrFhVGr389CJvjCv9+D7xp6AKL+B8qjL58r7EWF+xaSrITT5QZtR1c458O4kwoUhi1muUQsAGHHTwHPxu6glXP4uB1WViLC3VVguREeOhn0nawLwHEtu+1AphtnYUra8HEdqg1W1c92x7COf+dZMVlYeF/FZL5+YbWnXiizGy+r5zZB64JJeDjjX3QXVEisR1YQ18vL4Hz/5qBFS6VD5GbnDe4nHgOn7JDaq/UAZwNa4qhZzscQE+HiyW2g/unnl6zdVX8dxVDf9Y3uMCN/uR1B571spr1MXM/ZseuCrgzHt2P2ktfMbEdWIP/5ypgTh42sJqrOhtWAXNujsHlxBNlNm0Hx3D6ftTwJVrj3j+xX96NHXZgDf4kqYhz7LiBBRfykFveq+cZXE48p4M+LfT8iJYzKGsAs/JjDHzqFEj8nsBuvtdrli7mbmmS2Z4IuRfnDuC+stLglr3U91EB5K1xDOZ8M5Kd/LSlmz5z9ftuqhBiObgM7irDYnG2ajsziZ1YQ5+OK4N52zdWY1E/+WJxru8gmfe7HML7yesy8Cs6hL7OEsBenUK4h5QMYF4VFpn9ENuU/40COL/Shyzb8D3J6v4ZQH1OBv5+IIC+CoUQJzfw27+EcE96UQbrZWOBD5GcyD/IHvoM4zjxnlNDOKunwPdDM8TXphKn9BIk5h8Za+uynDwWdfKsDPruZAj3peyws8kc5GByeZLVXuJHZfFjbsjaS9l3XOJAfirepLGIc5oA+CvHol9ul5H3yFj4vcu2XfEuIdlhL/PkoE8v9pv3mheGLckigI0P492lQQT8Tmv7fW2XraueIZIZkzEC/XQ4jPM4m+YSOQJy9e0WN/pxcxh9UDsC83mF8JePwFxYHMZ7TYBk2Iq4lCE5sT3KL6+eE9wrrU5bv6Gp2PY4jTP4qtuYE7au4j/iRt111boWFvtcdMOutga3E0+UOc15de8UcpPCgzvW21OISS6PYYetq+o4pQfn2Suta2HRj14P3iMeGtxOPFHmMuzo+/QUZv9PN3r18ynMELdph62r/Erige4HS1dj4X82D+b7K4PbgefvI8zfE6qdxJ3ZlQMzt91J2Jwrh2GHravinyIH3iNba10bq+rKnQP1U8vgToSHfm5tB8e0Szzuis08OHOWxCPHg8x42LqqXtt7cM/qr3VtrIp/Tw/qYYbB7cTzK4a2g2N5Mx49UdqDnsiYgFpqZdph6SLOVaFb+2m86FpY+N/Qg7mTIsHmduD56dN2qHeQBJwhIY/9mzPv+8i0w9ZVOmHZZ+afseNuIF6T5xhcTjxf5cw6fXka75alEMtgZDzeGdqZdti6av/KHtz3ksb/hlV10UT6J5PB5cQTZV4jL0M3nMbMfubG2XXlNHogk2mHrav8eu+We4alq7GYc8nhf8ldBrcTT5T5tB1ca2uOI7ZuL3o+9jhiW8pr2GHpog69Xpyv+7SujVXfF/fiHFlncDvwjMtv/j774jhm7+UciJ3rBGrvjdm3li78v4Y+LPLjuOhaWJyXL3PAjrcGtxNPlAWM381HtD+BO9oM3OlHHCSZn/fN9wNbV90bVtGaeh/TuhYW/scXhJ3jHdw2XtlRkD5F5D5c/VYccrbPizP8WRx+swtLXpKyEf+hSyqFNC/XyqM4zNatXsyOD3GYcSEHr3PNiSWVwpr3D3u6Uh3FnfOuyfsfuiRG0qcS82INMbrsxUzKeRS/8T8R3uSi5+BRz0xa18aq/rjtRfwzG9zsi/LH6btT1+GPc82J/QVQSwMEFAAAAAgA7kiGTMSd7OjHAAAApAEAAAcAAABOTVQuc2h4hc6vCsJQFAbw4/w/tUytIitbFMM0uoVZ1bBk0eQTiPgAhiWTYWFv4LIMsQgGmQYfwGwwCVbxXO4XBQ/38OMeDnyHyFDpd90eaSJ+dEmSpLo17FrUn+zrun18m+MwMezncBc6Hd2m/9Xm7nHPYQw/0pQDV/DMDoiUBttiR3KubLibRJmYd3iUPch5ToEu9OFVmtegBwN4lzkFcd+aPbELoqK4Z8ZGcq/4kqoWXELklxToQh8iv6xBDwZQ5JtEFRVacMp22VD8v1BLAQI/ABQAAAAIAO5IhkyO288lcAAAACICAAAHACQAAAAAAAAAIAAAAAAAAABOTVQuZGJmCgAgAAAAAAABABgAubrR63XN0wG198vrdc3TAbX3y+t1zdMBUEsBAj8AFAAAAAgA7kiGTC3vY0IHAQAAkAEAAAcAJAAAAAAAAAAgAAAAlQAAAE5NVC5wcmoKACAAAAAAAAEAGACJRs/rdc3TAYlGz+t1zdMBiUbP63XN0wFQSwECPwAUAAAACADuSIZMvIcyHmEBAABvAgAABwAkAAAAAAAAACAAAADBAQAATk1ULnFwagoAIAAAAAAAAQAYAIlGz+t1zdMBiUbP63XN0wGJRs/rdc3TAVBLAQI/ABQAAAAIAO5Ihkxyru7zYQoAAJQbAAAHACQAAAAAAAAAIAAAAEcDAABOTVQuc2hwCgAgAAAAAAABABgAubrR63XN0wG198vrdc3TAeo5+ep1zdMBUEsBAj8AFAAAAAgA7kiGTMSd7OjHAAAApAEAAAcAJAAAAAAAAAAgAAAAzQ0AAE5NVC5zaHgKACAAAAAAAAEAGAC5utHrdc3TAbX3y+t1zdMBtffL63XN0wFQSwUGAAAAAAUABQC9AQAAuQ4AAAAA[/img]
Po zmianie wygląda dokładnie tak samo:

Cytat:from osgeo import ogr
from itertools import takewhile
from qgis.core import *
#sciezka do pliku z tekstem
plik = 'D:/Testy/NMT-N34107Cb1_c.asc'
warstwa_liniowa = QgsVectorLayer("LineString?crs=epsg:2180", 'test', "memory")
warstwa_liniowa.startEditing()
True
feat = QgsFeature(warstwa_liniowa.pendingFields())
linie = []
with open(plik) as f:
       array = []
       for line in f:
           if line.startswith('Start'):
               data = takewhile(lambda x: not x.startswith("End"),f)
               linia = ogr.Geometry(ogr.wkbLineString)
               for wiersz in data:
                   wsp = wiersz.split(' ')
                   linia.AddPoint(float(wsp[1]),float(wsp[0]) )
               gotowa_geom = linia.ExportToWkt()
               geom = QgsGeometry.fromWkt(gotowa_geom)
               feat=QgsFeature(warstwa_liniowa.pendingFields())
               feat.setGeometry(geom)
               warstwa_liniowa.addFeatures([feat])
warstwa_liniowa.commitChanges()
  File "<input>", line 15
    warstwa_liniowa.commitChanges()
                  ^
SyntaxError: invalid syntax
warstwa_liniowa.updateExtents()
QgsMapLayerRegistry.instance().addMapLayer(warstwa_liniowa, True)
<qgis._core.QgsVectorLayer object at 0x0000000024F80158>

Natomiast wygląda na to, że Twój skrypt działa.

Nie wiem, czy ma to jakieś znaczenie, ale mam wersję 2.18.12
Wersja tutaj nie ma znaczenia, ja używam 2.18.17.

Nie odpalaj go z konsoli tylko zapisz sobie do pliku za pomocą wbudowanego edytora - tak jak na moim zrzucie.
Układ sobie zmieniłeś, wygląda na to, że wpasowuje się we właściwe miejsce na OSM

https://screenpresso.com/=NqMOb
Rewelacja! Działa!

Pierwszy kontakt z konsolą boli.

Dziękuję xmaziax - zapraszam na piwo.
W sumie jeszcze jedna ważna sprawa, o której zapomniałem wspomnieć - jak zrobić, żeby podczas tworzenia nie gubiła się współrzędna "z"?
Linie generuje się świetnie, współrzędne węzłów są poprawne, ale ucieka rzędna która jest mi potrzebna...
Na wyjściu otrzymujesz 1 element zbudowany z wielu punktów. Punkty składowe mają różne wysokości, więc wynikowej wysokości nie można przyjąć bezpośrednio z punktów.

Poniżej masz wersję skryptu z wyliczonymi średnimi wartościami dla linii.

Kod:
from osgeo import ogr
from itertools import takewhile
from qgis.core import *
from qgis.PyQt.QtCore import *
from numpy import mean as srednia


#sciezka do pliku z tekstem
plik = 'h:/python/asci/nmt-n34107cb1_o.asc'
obiekty = []
i = 1
with open(plik) as f:  
       for line in f:
           if line.startswith('Start'):
               data = takewhile(lambda x: not x.startswith("End"),f)
               linia = ogr.Geometry(ogr.wkbLineString)
               wys = []
               for wiersz in data:
                   wsp = wiersz.split(' ')
                   linia.AddPoint(float(wsp[1]),float(wsp[0]) )
                   wys.append(float(wsp[2]))
               wys_sr = format(srednia(wys),'.2f')
               gotowa_geom = linia.ExportToWkt()
               geom = QgsGeometry.fromWkt(gotowa_geom)
               ob = [i, wys_sr, geom]
               obiekty.append(ob)
           i += 1

warstwa_liniowa = QgsVectorLayer("LineString?crs=epsg:2180", 'test', "memory")
pr = warstwa_liniowa.dataProvider()
warstwa_liniowa.startEditing()
pr.addAttributes([QgsField('ID', QVariant.Int)])
pr.addAttributes([QgsField('WYS', QVariant.Double)])
warstwa_liniowa.commitChanges()

caps = warstwa_liniowa.dataProvider().capabilities()
if caps & QgsVectorDataProvider.AddFeatures:
   feat=QgsFeature(warstwa_liniowa.pendingFields())
   for  feature in obiekty:
       feat.setAttribute('ID', feature[0])
       feat.setAttribute('WYS', feature[1])
       feat.setGeometry(feature[2])
       warstwa_liniowa.addFeatures([feat])
       (res, outFeats) = warstwa_liniowa.dataProvider().addFeatures([feat])

warstwa_liniowa.commitChanges()
warstwa_liniowa.updateExtents()
QgsMapLayerRegistry.instance().addMapLayer(warstwa_liniowa, True)
       
       
Dziękuję bardzo!

Dla niektórych linii nie ma problemu, ale niektóre mają zmianę rzędnych nawet kilka-kilkanaście metrów. Będzie to używane do tworzenia modelu 3D, więc muszę znaleźć metodę jak wykorzystać współrzędną "z" i przypisać ją do każdego węzła osobno.
Stron: 1 2