TransWikia.com

Create Lines between points using points ID

Geographic Information Systems Asked by daviddafonseca on May 1, 2021

I am trying to recreate in QGIS a few diferent sets of constellations. I already have a layer with all the stars I need, with ID and coordinates (downloaded from here astronomyData.zip), and I also have a table with a line for each constellation, and the pairs of stars that need to be connected by ID (link here)

What I need is for some way to have the software look at each pair of IDs, find the coordinate for each star with the same ID the stars layer, and connect each pair with a line.

I feel like this should be easy, but I can’t figure out a way.

2 Answers

There doesn't seem to be a simple way to do this with the GUI, because of the (paired) data structure your data have.

However, I found the use case interesting and had attempted a script before the other answer published (although, I didn't published it waiting for other answers using the GUI).

This is the result with an azimuthal projection:

enter image description here

And this is the script I ended up with:

stars = QgsProject.instance().mapLayersByName("stars6")[0]
stars_dict = {str(feature["HIP"]):feature.geometry() for feature in stars.getFeatures()}

# Prepare constellation layer
layer = QgsVectorLayer("MultiLineString?crs={}".format(stars.crs().authid()), "constellations", "memory")
layer.dataProvider().addAttributes([QgsField("name", QVariant.String)])
layer.updateFields()
new_features = list()

constellations_file = open('/path/to/constellationship.fab', 'r') 
for line in constellations_file:
    constellation_name, hip_part = line.split("  ")
    constellation_name = constellation_name.split(" ")[0]
    hip_list = hip_part.strip().split(" ")

    # Build constellation lines
    mls = QgsMultiLineString()
    for start, end in zip(hip_list[::2], hip_list[1::2]):
        if not (start in stars_dict and end in stars_dict):
            print("Pair {}-{} not found in stars layer! Skipping...".format(start, end))
            continue
            
        g = QgsGeometry.fromPolyline([stars_dict[start].get(), stars_dict[end].get()])
        mls.addGeometry(g.get().clone())

    new_feature = QgsFeature()
    new_feature.setGeometry(mls)
    new_feature.setAttributes([constellation_name])
    new_features.append(new_feature)

constellations_file.close()
layer.dataProvider().addFeatures(new_features)
QgsProject.instance().addMapLayer(layer)

It creates a segment per each pair of stars; it doesn't create segments with length 0; and, the resulting layer has the corresponding constellation name.

Issues with the data

  • In constellationship.fab, the line CMa 17... doesn't have two spaces after number of points (CMa 17 33160 34045...). I had to adjust it manually (find the modified file here).
  • Star 33165 does not exist!

Usage

  1. Load the stars6.shp file into QGIS.
  2. Adjust the /path/to/constellationship.fab to reflect your local path. (Use the modified file mentioned above when talking about issues with the data.)
  3. Run the script from the QGIS Python console.

Correct answer by Germán Carrillo on May 1, 2021

I dont think this is 100 % correct but maybe you can fix it or understand whats causing it? Some constellations look a bit weird with line going across the map.

import csv

pointlyr = QgsProject.instance().mapLayersByName('StarsNamed2')[0]
pointdict = {f['Stars_Named_HIP']:QgsPoint(f.geometry().asPoint()) for f in pointlyr.getFeatures()} # Create a dictionary with star_id (?) as key and QgsPoint as value
orderfile = '/home/bera/Downloads/Linespoints/constellationship.fab'

#Create a line layer in memory
vl = QgsVectorLayer("LineString?crs='EPSG:4326'&index=yes", "constellations_temp", "memory")
provider = vl.dataProvider()

errorlist = []

with open(orderfile, newline='') as csvfile: #For each line in the text file find the matching point in the dictionary, store in a list and create a line
    spamreader = csv.reader(csvfile, delimiter=' ', quotechar='|')
    for row in spamreader:
        gLine = []
        for star in row[3:]:
            try:
                if float(star) in pointdict:
                    gLine.append(pointdict[float(star)])
            except:
                errorlist.append(row)
        gLine = QgsGeometry.fromPolyline(gLine)
        f = QgsFeature()
        f.setGeometry(gLine)
        provider.addFeature(f)
QgsProject.instance().addMapLayer(vl)
print(errorlist) #The text file isnt perfect, one error line

enter image description here

Answered by BERA on May 1, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP