uber / h3

Hexagonal hierarchical geospatial indexing system
https://h3geo.org
Apache License 2.0
4.83k stars 459 forks source link

Visualization problems with center to center connections. #759

Closed SkybuckFlying closed 1 year ago

SkybuckFlying commented 1 year ago

I am trying to draw/render lines from cell center to cell center... See picture below...

In case you wanna see crazy video, check here:

https://www.youtube.com/live/-EsixpO9mBk?feature=share

Anyway it would be great if this library had some support to do things like this more easily...

I am currently use code from c++ geodesy https://github.com/sslnjz/geodesy Converted to Delphi and such... and my own function ClosestPointDistanceToClosestPointOnSegment below to try and accomplish this, but all of this math seems to collapse under certain conditions, not sure if that is truely happening but seems so... or very maybe uber returning some weird data, or some computation is failing:

Plus some x,y,z to lat, long conversion code...

// engineer who provided it to you for assistance. // Define some earth constants const MAJA = (6378137.0); // semi major axis of ref ellipsoid const FLAT = (1.0/298.2572235630); // flattening coef of ref ellipsoid. // These are derived values from the above WGS-84 values const ESQR = (FLAT (2.0-FLAT)); // 1st eccentricity squared const OMES = (1.0 - ESQR); // 1 minus eccentricity squared const EFOR = (ESQR ESQR); // Sqr of the 1st eccentricity squared const ASQR = (MAJA * MAJA); // semi major axis squared := nlmaja**

procedure ConvertECEFToLTP ( const ParaEcef : TEcef; var ParaLtp : TLtp ); var a0,a1,a2,a3,a4,b0,b1,b2,b3 : double; b,c0,opqk,qk,qkc,qks,f,fprm : double; k : integer; ytemp, xtemp : double; begin // Convert from ParaEcef (XYZ) to LTP (ParaLtp.Lat/ParaLtp.Lon/ParaLtp.Alt)

//------------------------------------------- // b := (00 + 11) / semi_major_axis_squared // c := (zz) / semi_major_axis_squared //------------------------------------------- b := ( (ParaEcef.x ParaEcef.x) + (ParaEcef.y ParaEcef.y) ) / ASQR; c0 := (ParaEcef.z ParaEcef.z) / ASQR;

//------------------------------------------- // a0 := c one_minus_eccentricity_sqr // a1 := 2 a0 // a2 := a0 + b - first_eccentricity_to_fourth // a3 := -2.0 first_eccentricity_to_fourth // a4 := - first_eccentricity_to_fourth //------------------------------------------- a0 := OMES c0; a1 := 2.0 a0; a2 := a0 + b - EFOR; a3 := -2.0 EFOR; a4 := - EFOR;

//------------------------------------------- // b0 := 4 a0, b1 := 3 a1 // b2 := 2 a2, b3 := a3 //------------------------------------------- b0 := 4.0 a0; b1 := 3.0 a1; b2 := 2.0 a2; b3 := a3;

//------------------------------------------- // Compute First Eccentricity Squared //------------------------------------------- qk := ESQR;

// for (k := 1; k <:= 3; k++) for k := 1 to 3 do begin qks := qk qk; qkc := qk qks; f := (a0 qks qks) + (a1 qkc) + (a2 qks) + (a3 qk) + a4; fprm := (b0 qkc) + (b1 qks) + (b2 qk) + b3; qk := qk - (f / fprm); end;

//------------------------------------------- // Compute latitude, longitude, altitude //------------------------------------------- opqk := 1.0 + qk; if ( (ParaEcef.x = 0.0) and (ParaEcef.y = 0.0) ) then // on the earth's axis begin // We are sitting EXACTLY on the earth's axis // Probably at the center or on or near one of the poles ParaLtp.Lon := 0.0; // as good as any other value if (ParaEcef.z >= 0.0) then begin ParaLtp.Lat := PI / 2; // ParaLtp.Alt above north pole end else begin ParaLtp.Lat := -PI / 2; // ParaLtp.Alt above south pole end; end else begin ytemp := opqk ParaEcef.z; xtemp := sqrt ( ( ParaEcef.x ParaEcef.x ) + ( ParaEcef.y ParaEcef.y ) ); ParaLtp.Lat := arctan2( ytemp, xtemp ); ParaLtp.Lon := arctan2( ParaEcef.y, ParaEcef.x ); end; ParaLtp.Alt := ( 1.0 - OMES / ESQR qk ) MAJA sqrt ( ( b / (opqk * opqk) ) + c0 );

end; // ConvertECEFToLTP()

// written by Skybuck Flying ! ;) :) =D function ClosestPointDistanceToClosestPointOnSegment ( ParaPointX, ParaPointY, ParaPointZ : double; ParaStartX, ParaStartY, ParaStartZ : double; ParaEndX, ParaEndY, ParaEndZ : double; ParaRadius : double {= 6371e3} // (Mean) radius of earth (defaults to radius in metres). ) : double; var vPoint, vSegmentPoint1, vSegmentPoint2, vClosestPointOnSegment : TLatLonNvectorSpherical; begin vSegmentPoint1.FromXYZ2( ParaStartX, ParaStartY, ParaStartZ, ParaRadius ); vSegmentPoint2.FromXYZ2( ParaEndX, ParaEndY, ParaEndZ, ParaRadius ); vPoint.FromXYZ2( ParaPointX, ParaPointY, ParaPointZ, ParaRadius );

vClosestPointOnSegment := nearestPointOnSegment ( vPoint, vSegmentPoint1, vSegmentPoint2, ParaRadius );

result := distanceTo ( vPoint, vClosestPointOnSegment, ParaRadius );

end;

function ClosestPointDistanceToClosestPointOnSegmentV2 ( ParaPointX, ParaPointY, ParaPointZ : double; ParaStartX, ParaStartY, ParaStartZ : double; ParaEndX, ParaEndY, ParaEndZ : double; ParaRadius : double {= 6371e3} // (Mean) radius of earth (defaults to radius in metres). ) : double; var vPoint, vSegmentPoint1, vSegmentPoint2, vClosestPointOnSegment : TLatLonNvectorSpherical; vEcef : Tecef; vLtp : Tltp; begin // convert start x,y,z to lat lng vEcef.X := ParaStartX; vEcef.Y := ParaStartY; vEcef.Z := ParaStartZ;

ConvertECEFToLTP ( vEcef, vLtp );

vSegmentPoint1.mLat := vLtp.Lat; vSegmentPoint1.mLon := vLtp.Lon;

// convert end x,y,z to lat lang vEcef.X := ParaEndX; vEcef.Y := ParaEndY; vEcef.Z := ParaEndZ;

ConvertECEFToLTP ( vEcef, vLtp );

vSegmentPoint2.mLat := vLtp.Lat; vSegmentPoint2.mLon := vLtp.Lon;

// ParaRadius ??

// convert point x,y,z to lat lng vEcef.X := ParaPointX; vEcef.Y := ParaPointY; vEcef.Z := ParaPointZ;

ConvertECEFToLTP ( vEcef, vLtp );

vPoint.mLat := vLtp.Lat; vPoint.mLon := vLtp.Lon;

vClosestPointOnSegment := nearestPointOnSegment ( vPoint, vSegmentPoint1, vSegmentPoint2, ParaRadius );

result := distanceTo ( vPoint, vClosestPointOnSegment, ParaRadius );

end;

Any idea what could be wrong ? If you need more code let me know...

ProblemV1

isaacbrodsky commented 1 year ago

As I understand this issue is about your rendering code -- I don't think this is the place to get help on that.