0

For a project I have the lovely joy to create an equally spaced grid, to be placed on a map later in QGIS. To this end, I created 4 reference coordinates, defining the top left corner as well as the x and y spacing. Based on this grid, my Python programming yielded this extended grid:

Not exactly what is expected

As you can see it has several issues:

  1. Both in X and Y direction it looks wobbly, the points are not perfectly aligned
  2. The distance between the points is not the required distance

Problem 1 is, to my opinion, caused by the extension approach and losing significance in the digits. Problem 2 is caused due to a bad definition of coordinate sets.

However, I choose this approach as I did not find a way to define a rotated grid based on coordinates.

Preferably, I have following inputs to my script:

  • Starting coordinate (i.e. top left)
  • horizontal and vertical spacing
  • number of rows / columns

What should be adopted to mitigate before mentioned issues or, in the ideal world, how could a solution direction look like based on the preferred inputs?

from string import ascii_lowercase
from geojson import Point, Feature, FeatureCollection, dump

def getExtraPolatedPoint(p1,p2, ratio):
    'Creates a line extrapoled in p1->p2 direction'
    b = (p1[0]+ratio*(p2[0]-p1[0]), p1[1]+ratio*(p2[1]-p1[1]) )
    return b

L = [letter1+letter2 for letter1 in ascii_lowercase for letter2 in ascii_lowercase]

with open(r"4points.geojson") as json_file:
    res = json.load(json_file)
    
features = []
feature_collection = FeatureCollection(features)

tl = None

for a in res["features"]:
    b = a["geometry"]["coordinates"]
    if tl is not None:
        if tl[1]<b[1]:tl=b
        if tr[0]<b[0]:tr=b
        if br[1]>b[1]:br=b
        if bl[0]>b[0]:bl=b
    else:
        tl = tr = br = bl = b

refP1 = tl
refP2x = tr
refP2y = bl
refP2xy = br

for i in range(16):
    newP1  = getExtraPolatedPoint(refP1,refP2x,i)
    newP2xy = getExtraPolatedPoint(refP2y,refP2xy,i)
    
    for j in range(75):
        gridPoint =  getExtraPolatedPoint(newP1,newP2xy,j)
        point = Point(gridPoint)
        features.append(Feature(geometry=point, properties={"ID":L[i]+'{:02}'.format(j+1) }))

Update: It is a geo coordinate system, I.e. lat/lon coordinates. Rotation comes from the fact that we need to fit the lines on an object on the map, hence the above rotation.

Update 2: So, the pyproj was the helping hand. The missing term which connects the dots (no pun intended) is azimuth. Using this package I was able to create the maximum coordinates using the fwd function. Having that in place, the npts function created the necessary points in the spacing.

Sometimes thing are obvious, if you know what you are looking for...

significant
  • 153
  • 7
  • in what form are coordinates and spacing defined? Where does the rotation come in? – Marat May 10 '23 at 22:55
  • Thanks for asking, we look at lat/lon coordinates and rotation is due to the necessity to fit the points to points on an object on earth. Hence, in Qgis the result wil be used as an overlay. – significant May 11 '23 at 04:41
  • How is spacing defined - e.g., is it meters or degrees? What is the center of rotation - e.g., is it the start coordinate, center of all points, or something else? What is the angle of rotation (it is not in parameters)? – Marat May 11 '23 at 13:40
  • spacing is defined in meters (parameterized), rotation is due to the fact that underlaying object on the satilite image is on an angle. This angle is unknown, but if parameterised one can in some loops come to an estimated guess. – significant May 11 '23 at 18:47
  • sorry, maybe I'm missing something obvious - but how are you going to rotate these points without knowing the angle? Whose responsibility is to estimate this angle? – Marat May 11 '23 at 19:15
  • @Marat see update in the opening post – significant May 15 '23 at 19:49

1 Answers1

0

Working with distances in geographic units should be avoided unless the values are very large, this is most likely causing your issues.

My suggestion is to (1) project your start coordinates, (2) define the end ones according to the distance, number of total points per path and rotation (3) create a line between the start and end points, (4) use shapely interpolate to define the intermediate points, (5) convert the points back to geographic coordinates. Any remaining wobbly looks are probably due to the conversion between the coordinates / float precision

Below is a sample code for one path:

# imports
from shapely.geometry import LineString
import geopandas as gpd
import numpy as np
import pyproj

# User input
start_lat = 48.8566  # start latitude in degrees 
start_lon = 2.3522  # start longitude in degrees
dist = 10 # distance between points
n_pts = 100  # total number of points for each "line"
rotat = 260  # the rotation required in degrees
proj_epsg = 32631  # the EPSG reference for a projected coordinate system

# pyproj transformer from WGS84 to the projected system
geo2proj = pyproj.Transformer.from_crs(crs_from=4326, crs_to=proj_epsg, always_xy=True)

# project the start point
start_x, start_y = geo2proj.transform(start_lon, start_lat)

# calculate the end points considering the rotation and the distance
end_x = start_x + dist * n_pts * np.cos(np.deg2rad(rotat))
end_y = start_y + dist * n_pts * np.sin(np.deg2rad(rotat))

# create a path between the start and the end points
rotated_path = LineString([(start_x, start_y), (end_x, end_y)])

# create a geodataframe to store the outputs
gdf = gpd.GeoDataFrame(geometry=[], crs=proj_epsg)

# using Shapely interpolate insures the points will be on the path
for n in range(0, n_pts):
    gdf.loc[n, 'geometry'] = rotated_path.interpolate(n * dist)

# projected coordinates to geographic
gdf.to_crs(4326, inplace=True)
Isadora Rezende
  • 51
  • 1
  • 2