Lisp to draw a line between GPS points

 
;draw real distances from GPS coordinates
;load points from csv files
(Defun C:gpsd ( / dataL csv_name ofil rd-pos n Pi1 P1 cntr P2n ang ED lat1 lat2 Dlat Dlong a c RD ent entdata)
;*******************load files***************
(setq dataL '());define an empty list
(setq csv_name "");define an empty name for file
(setq csv_name (getfiled "Select Point List .csv file" "" "csv" 2));get the file name 
(setq ofil (open CSV_name "r")); define "ofile" as open the file
(while (setq rd-pos (read-line ofil));read each line from file
(setq dataL (cons rd-pos dataL));add each row to list
);close while
(close ofil);close file
(setq dataL (reverse dataL));put list in right order
(setq n (length dataL));get the number of points
(setq Pi1 (nth 0 dataL));get initial GPS point
(setq P1 '(0 0));set initial Real point to "0,0"
(setq xi1 (car Pi1));get initial x from first point
(setq yi1 (cadr Pi1));get initial y from first point
(setq cntr 1)
(while (<= cntr n)
(setq P2n (nth cntr dataL));choose second GPS point by cntr
(setq Pxn (car P2n));get x from second point
(setq Pyn (cadr P2n));get y from second point
(setq ang (angle Pi1 P2n));get angle between points
(setq ang (* 180.0 (/ ang pi)));convert angles from radians to degrees
;calculate real distance
(setq ED 6371000);earth diameter in metres
(setq lat1 (* pi (/ xi1 180.0)));lat in rad
(setq lat2 (* pi (/ Pxn 180.0)))
(setq Dlat (* pi(/ (- Pxn xi1) 180.0)))
(setq Dlong (* pi(/ (- Pyn yi1) 180.0))) 
(setq a (+ (sin(/ 2 Dlat)^2) (* (* (cos lat1) (cos lat2)) (sin (/ 2 Dlong)^2))))
(setq c (* 2 (atan (sqrt a) (sqrt (- 1 a)))));calculate earth angle
(setq RD (* ED c));calculate real distance between GPs Points
(setq RD (* RD 1000));convert distance from metres to mm
(setq P2 (polar P1 ang RD));calculate Real P2 by distance and angle
(command ".line" P1 P2 "");draw a line as real distance between GPS points
(setq ent (entlast));select last entoty created
(setq entdata (entget ent));extract info from it
(setq P1 (cdr (assoc 11 entdata)));reset Real P1 as (X2,Y2) from last line
(setq Pi1 P2n);reset GPS Pi1 as next point
(setq xi1 (car Pi1));recalculate x from new point
(setq yi1 (cadr Pi1));recalculate y from new point
;reset counter
(setq cntr (+ 1 cntr)) 
);close while loop
);close lisp

Please read the CODE POSTING GUIDELINES and edit your post to include CODE TAGS.

When you open CSV file in notepad - not excel, you'll notice that every cell in a row of CSV is separated by "," character, so you must find Xcoord and Ycoord separately and put them in list of points as these coordinates are later retrieved with car or cadr functions... Also ^ operand in lisp is unknown - check (expt) function for exponentiation...

Not checked but these are my remarks... (highlighted) - you must test it and later improve...

;draw real distances from GPS coordinates
;load points from csv files
(Defun C:gpsd ( / Xcoord Ycoord dataL csv_name ofil rd-pos n Pi1 P1 cntr P2n ang ED lat1 lat2 Dlat Dlong a c RD ent entdata)
;*******************load files***************
(setq dataL '());define an empty list
(setq csv_name "");define an empty name for file
(setq csv_name (getfiled "Select Point List .csv file" "" "csv" 2));get the file name 
(setq ofil (open CSV_name "r")); define "ofile" as open the file
(while (setq rd-pos (read-line ofil));read each line from file
[highlight]
(setq Xcoord (substr rd-pos 1 (vl-string-position (ascii ",") rd-pos)))
(setq Ycoord (substr rd-pos (+ (vl-string-position (ascii ",") rd-pos) 2) (- (strlen rd-pos) (vl-string-position (ascii ",") rd-pos) 1)))
(setq dataL (cons (list (read Xcoord) (read Ycoord)) dataL));add each row to list
[/highlight]
);close while
(close ofil);close file
(setq dataL (reverse dataL));put list in right order
(setq n (length dataL));get the number of points
(setq Pi1 (nth 0 dataL));get initial GPS point
(setq P1 '(0 0));set initial Real point to "0,0"
(setq xi1 (car Pi1));get initial x from first point
(setq yi1 (cadr Pi1));get initial y from first point
(setq cntr 1)
(while (<= cntr n)
(setq P2n (nth cntr dataL));choose second GPS point by cntr
(setq Pxn (car P2n));get x from second point
(setq Pyn (cadr P2n));get y from second point
(setq ang (angle Pi1 P2n));get angle between points
(setq ang (* 180.0 (/ ang pi)));convert angles from radians to degrees
;calculate real distance
(setq ED 6371000);earth diameter in metres
(setq lat1 (* pi (/ xi1 180.0)));lat in rad
(setq lat2 (* pi (/ Pxn 180.0)))
(setq Dlat (* pi(/ (- Pxn xi1) 180.0)))
(setq Dlong (* pi(/ (- Pyn yi1) 180.0))) 
[highlight](setq a (+ (expt (sin (/ 2 Dlat)) 2) (* (* (cos lat1) (cos lat2)) (expt (sin (/ 2 Dlong)) 2))))[/highlight]
(setq c (* 2 (atan (sqrt a) (sqrt (- 1 a)))));calculate earth angle
(setq RD (* ED c));calculate real distance between GPs Points
(setq RD (* RD 1000));convert distance from metres to mm
(setq P2 (polar P1 ang RD));calculate Real P2 by distance and angle
(command ".line" P1 P2 "");draw a line as real distance between GPS points
(setq ent (entlast));select last entoty created
(setq entdata (entget ent));extract info from it
(setq P1 (cdr (assoc 11 entdata)));reset Real P1 as (X2,Y2) from last line
(setq Pi1 P2n);reset GPS Pi1 as next point
(setq xi1 (car Pi1));recalculate x from new point
(setq yi1 (cadr Pi1));recalculate y from new point
;reset counter
(setq cntr (+ 1 cntr)) 
);close while loop
);close lisp

M.R.

I have marked in red the trouble maker code - the lines you got from data file are strings, while you attempt to treat them as list:

...
(setq Pi1 (nth 0 dataL));get initial GPS point
(setq P1 '(0 0));set initial Real point to "0,0"
[color=red](setq xi1 (car Pi1));get initial x from first point
(setq yi1 (cadr Pi1));get initial y from first point[/color]
...

Thanks, I will try tomorrow Marko's version, I was thinking that I am reading the each row from csv as a list of coordinates for each point, obviously I am wrong.

Regards,

Cristian

I need your help again, this time I believe the angle gets a negative value and my program stops, I tried to force my calculation to bring me positive values using "abs", but still nothing

I have also change my code a bit, I am now calculating the distances relative to the first poing, get the real x&y, I store them into a list and at the end I want to draw a line to connect all the dots.

 
;draw real distances from GPS coordinates
;load points from csv files
(Defun C:gpx ( / Xcoord Ycoord dataL Rpt csv_name ofil rd-pos n Pi1 P1 cntr P2n ang ED lat1 lat2 Dlat Dlong a c RD ent entdata)
;*******************load files***************
(setq dataL '());define an empty list
(setq Rpt '());define an empty list for points
(setq csv_name "");define an empty name for file
(setq csv_name (getfiled "Select Point List .csv file" "" "csv" 2));get the file name 
(setq ofil (open CSV_name "r")); define "ofile" as open the file
   (while (setq rd-pos (read-line ofil));read each line from file
       (setq Xcoord (substr rd-pos 1 (vl-string-position (ascii ",") rd-pos)))
       (setq Ycoord (substr rd-pos (+ (vl-string-position (ascii ",") rd-pos) 2) (- (strlen rd-pos) (vl-string-position (ascii ",") rd-pos) 1)))
       (setq dataL (cons (list (read Xcoord) (read Ycoord)) dataL));add each row to list
   );close while
(close ofil);close file
(setq dataL (reverse dataL));put list in right order
(setq n (length dataL));get the number of points
(setq Pi1 (nth 0 dataL));get initial GPS point
(setq P1 '(0 0));set initial Real point to "0,0"
(setq xi1 (car Pi1));get initial x from first point
(setq yi1 (cadr Pi1));get initial y from first point
(setq Rpt (cons P1 Rpt));add P1 on list as 0,0
(setq cntr 1)
(while (<= cntr n)
   (setq P2n (nth cntr dataL));choose second GPS point by cntr
   (setq Pxn (car P2n));get x from second point
   (setq Pyn (cadr P2n));get y from second point
   (setq ang (angle Pi1 P2n));get angle between points
   (setq ang (* 180.0 (/ ang pi)));convert angles from radians to degrees
;calculate real distance
  (setq ED 6371000);earth diameter in metres
  (setq lat1 (abs(* pi (/ xi1 180.0))));lat in rad
  (setq lat2 (abs(* pi (/ Pxn 180.0))))
  (setq Dlat (abs (* pi(/ (- Pxn xi1) 180.0))));abs added 
  (setq Dlong (abs (* pi(/ (- Pyn yi1) 180.0))));abs added 
  (setq a (+ (expt (sin (/ 2 Dlat)) 2) (* (* (cos lat1) (cos lat2)) (expt (sin (/ 2 Dlong)) 2))))
  (setq c (* 2 (atan (sqrt a) (sqrt (- 1 a)))));calculate earth angle
  (setq RD (* ED c));calculate real distance between GPs Points
  (setq RD (* RD 1000));convert distance from metres to mm
  (setq P2 (polar P1 ang RD));calculate Real P2 by distance and angle
  (setq Rpt (cons P2 Rpt));add next point to the list  
;reset counter
 (setq cntr (+ 1 cntr)) 
);close while loop
(setq Rpt (reverse RPT));put the list in the right order
(command ".line" Rpt "");draw a line as real distance between GPS points
);close lisp

"; error: function undefined for argument: -0.0577494"

for points:

56.06613 3.4514

56.06671 3.45163

56.06694 3.44979

56.0663 3.44959

and

; error: bad argument type: 2D/3D point: nil

for this list of points:

56.045809 -3.467674

56.046754 -3.467043

56.047555 -3.466557

56.049877 -3.465697

what is wrong in my code?

Thanks,

Cristian

First, please pay attention to the number of items you try to walk on the list:

[s][color=red](while (<= cntr n)[/color][/s]
(while (< cntr n)

Second, the line cannot be drawn by using the list of points as argument:

[s][color=red](command ".line" Rpt "")[/color][/s]
(command "_.LINE")
(foreach point Rpt
(command "_non" point)
)
(command "")

Thanks Mircea,

I have changed my code with your recomandations but i still get this error: "; error: function undefined for argument: -0.00160928", I have attached my list of points as well as text file.

Thanks again,

Cristian

home road csv.2 .txt

There is an inconsistency between the code and example file – the separator used in the file is instead of comma considered in code; so you should use;

(ascii "\t")

instead of:

(ascii "\,")

The problem is from the fact that square root have meaning only for positive numbers; in at least one case the argument become negative.

(setq c (* 2 (atan (sqrt a) (sqrt [color=red](- 1 a)[/color]))))

You should validate again your formulas - I cannot check them since don't know from where come from.

The (at least first) trouble-maker pair of data is:

56.0534 -3.46219

The problem is from the fact that square root have meaning only for positive numbers; in at least one case the argument become negative.

(setq c (* 2 (atan (sqrt a) (sqrt [color=red](- 1 a)[/color]))))

You should validate again your formulas - I cannot check them since don't know from where come from.

 

The (at least first) trouble-maker pair of data is:

This is haversine formula, you may find more details in here.

Thanks,

Cristian

I've changed your home road.csv - added "," as separator, and my first code for drawing peace of Earth ground on sphere did the job, but second code gave wrong results, as there were many soilds that were unable to be unioned and after that intersect with sphere anvelope with 1m thick shell...

Here are my codes and your csv - modified by me...

Regards, M.R.

P.S. I've used Evgeniy Elpanov's triangulation for determining area code, but no use on your example, you must find some other way to determine area...

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;                                                                                 ;;;
;;; gpspronsph function takes GPS coordinates from CSV file,                        ;;;
;;; and project between points edges as arcs on sphere with radius of planet Earth. ;;;
;;;                                                                                 ;;;
;;; 25.07.2012. Marko Ribar, d.i.a. - first and final release                       ;;;
;;;                                                                                 ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun c:gpspronsph ( / file f csvrow Latt Long gpsdata ER ptx pty ptz pt ptlst k p1 p2 pm pmv1 pmsph )
 (vl-cmdf "_.UCS" "w")
 (if (ssget "_X")
   (progn 
     (alert "\nDWG must be empty - empty it and start over")
     (exit)
   )
 )
 (setq file (getfiled "Select file with GPS coordiantes" "" "csv" 2))
 (setq f (open file "r"))
 (while (setq csvrow (read-line f))
   (setq Latt (substr csvrow 1 (vl-string-position (ascii ",") csvrow)))
   (setq Long (substr csvrow (+ (vl-string-position (ascii ",") csvrow) 2) (- (strlen csvrow) (vl-string-position (ascii ",") csvrow) 1)))
   (setq gpsdata (cons (list (read Latt) (read Long)) gpsdata))
 )
 (setq gpsdata (reverse gpsdata))
 (setq ER 6371000.0) ; Radius of Earth in m 
 (foreach gps gpsdata
   (setq Latt (car gps))
   (setq Long (cadr gps))
   (setq ptx (* ER (cos (cvunit Latt "degrees" "radians")) (cos (cvunit Long "degrees" "radians"))))
   (setq pty (* ER (cos (cvunit Latt "degrees" "radians")) (sin (cvunit Long "degrees" "radians"))))
   (setq ptz (* ER (sin (cvunit Latt "degrees" "radians"))))
   (setq pt (list ptx pty ptz))
   (setq ptlst (cons pt ptlst))
 )
 (setq ptlst (reverse ptlst))
 (setq k -1)
 (repeat (length ptlst)
   (setq p1 (nth (setq k (1+ k)) ptlst))
   (if (/= k (- (length ptlst) 1)) 
     (setq p2 (nth (+ k 1) ptlst)) 
     (setq p2 (car ptlst))
   )
   (setq pm (mapcar '(lambda ( a b ) (/ (+ a b) 2.0)) p1 p2))
   (setq pmv1 (mapcar '(lambda ( a ) (/ a (distance '(0.0 0.0 0.0) pm))) pm))
   (setq pmsph (mapcar '(lambda ( a ) (* a ER)) pmv1))
   (vl-cmdf "_.UCS" "3" '(0.0 0.0 0.0) p1 p2)
   (vl-cmdf "_.ARC" (trans p1 0 1) "s" (trans pmsph 0 1) (trans p2 0 1))
   (vl-cmdf "_.UCS" "p")
 )
 (vl-cmdf "_.ZOOM" "_E")
 (princ)
)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;                                                                                 ;;;
;;; gpsareasph function takes GPS coordinates from CSV file,                        ;;;
;;; and project between points edges as arcs on sphere with radius of planet Earth, ;;;
;;; then it asks for offset distance to determine contour shape of peace of ground, ;;;
;;; for witch if finds approx. correct area value.                                  ;;;
;;; 25.07.2012. Marko Ribar, d.i.a. - first and final release                       ;;;
;;;                                                                                 ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(vl-load-com)
;;********************************************************;
;;                                                        ;
;; Written by  ElpanovEvgeniy                             ;
;; 17.10.2008                                             ;
;; edit 20.05.2011                                        ;
;; Program triangulate an irregular set of 3d points.     ;
;; Modified by ymg June 2011                              ;
;;********************************************************;
(defun triangulate (pl  /  a   b   c   i   i1   i2 
     bb  sl al  el  tl  l   ma  mi     
     ti  tr x1  x2  y1  y2  p    r  cp
    )
  
  (if pl
     (progn
(setq ti (car (_VL-TIMES))
       i  1
       i1 (/ (length pl) 100.)
       i2 0
       pl (vl-sort pl
      (function (lambda (a b) (< (car a) (car b))))
   ) 
       bb (list (apply 'mapcar (cons 'min pl))
                (apply 'mapcar (cons 'max pl))
                 )
              ;Replaced code to get minY and maxY with 3d Bounding Box Routine;
       ;A bit slower but clearer. minZ and maxZ kept for contouring    ;
       
       x1 (caar bb)   ;minX 
       x2 (caadr bb)  ;maxX 
       y1 (cadar bb)  ;minY 
       y2 (cadadr bb) ;maxY 
       
       
)
(setq cp (list (/ (+ x1 x2) 2.0) (/ (+ y1 y2) 2.0)); Midpoint of points cloud and center point of circumcircle through supertriangle.
        r (* (distance cp (list x1 y1)) 20) ;This could still be too small in certain case. No harm if we make it bigger.
       ma (+ (car cp) r);ma is maxX of supertriangle
       mi (- (car cp) r);mi is minX of supertriangle
       sl (list (list ma (cadr cp) 0);list of 3 points defining the supertriangle
  (list mi (+ (cadr cp) r) 0)
  (list mi (- (cadr cp) r) 0)
   ) 
   
       al  (list (cons x2 (cons cp (cons (* 20 r) sl))))
      ;al is a work list that contains active triangles each item is a list that contains:      ;
      ;     item 0: Xmax of points in triangle.                            ;
      ;     item 1: List 2d coordinates of center of circle circumscribing triangle.    ;
      ;     item 2: Radius of above circle.                         ;
      ;     item 3: List of 3 points defining the triangle                            ;
 
       ma (1- ma);Reducing ma to prepare for removal of triangles having a vertex
       mi (1+ mi);common with supertriangle. Done once triangulation is completed.
)               ;Increasing mi for same reason.

;Begin insertion of points

(repeat (length pl) 
    
    (setq  p (car pl);Get one point from point list
   pl (cdr pl);Remove point from point list
   el nil     ;Initialize edge list
    ) 
    (while al ;Active triangle list
       (setq tr (car al);Get one triangle from active triangle list.
      al (cdr al);Remove the triangle from the active list.
       ) 
       (cond
   ((< (car tr) (car p)) (setq tl (cons (cdddr tr) tl)));This triangle inactive. We store it's 3 vertex in tl (Final triangle list).
   ((< (distance p (cadr tr)) (caddr tr));p is inside the triangle.
    (setq tr (cdddr tr) ;Trim tr to vertex of triangle only.
    a (car tr)   ;  First point.
    b (cadr tr)  ;  Second point.
    c (caddr tr) ;  Third point.
          el (cons (list (+ (car a) (car b))   ;We create 3 new edges ab, bc and ca, 
           (+ (cadr a) (cadr b)) ;and store them in edge list.
    a
    b
     )
     (cons (list (+ (car b) (car c))
          (+ (cadr b) (cadr c))
          b
          c
    )
    (cons (list (+ (car c) (car a))
         (+ (cadr c) (cadr a))
         c
         a
          )
          el
    )
     ) 
      )
      
    )
   )
   (t (setq l (cons tr l)))   ;tr did not meet any cond so it is still active. 
       )         ;we store it a temporary list.
    );Go to next triangle of active list.
    
    (setq al l    ;Restore active triangle list from the temporary list.
    l nil  ;Re-initialize the temporary list to prepare for next insertion.
   
   el (vl-sort el ;Sort the edges list on X and Y
          (function (lambda (a b)
         (if (= (car a) (car b))
            (<= (cadr a) (cadr b))
            (< (car a) (car b))
         ) 
      ) 
          ) 
      ) 
    ) 
    ;Removes doubled edges, form new triangles, calculates circumcircles and add them to active list.
    (while el
       (if (and (= (caar el) (caadr el))
  (= (cadar el) (cadadr el))
    )
   (setq el (cddr el)) 
   (setq al (cons (getcircumcircle p (cddar el)) al)
  el (cdr el)
   ) 
       ) 
    ) 
    ;Spinner to show progress
    (if (and (< (setq i (1- i)) 1) (< i2 100))
       (progn
   (setvar
      "MODEMACRO"
      (strcat
  "     "
  (itoa (setq i2 (1+ i2)))
  " %    "
  (substr
     "||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||"
     1
     i2
  ) 
  (substr "..." 1 (- 100 i2))
      ) 
   ) 
   (setq i i1)
       ) 
    ) 
) ;Go to insert next point.
;We are done with the triangulation 

(foreach tr al (setq tl (cons (cdddr tr) tl)));What's left in active list is added to triangle list

;Purge triangle list of any triangle that has a common vertex with supertriangle.
(setq
    tl (vl-remove-if-not
   (function
      (lambda (a)
  (and (< mi (caadr a) ma) (< mi (caaddr a) ma))
      )
   )
   tl
       ) 
)
(foreach tr tl
           (vl-cmdf "_.UCS" "3" (car tr) (cadr tr) (caddr tr))
           (vl-cmdf "_.PLINE" (trans (car tr) 0 1) (trans (cadr tr) 0 1) (trans (caddr tr) 0 1) "c" "")
           (vl-cmdf "_.UCS" "p")
) 
     ) 
  ) 
  (setvar "MODEMACRO" "")
  (princ (strcat "\n "
   (rtos (/ (- (car (_VL-TIMES)) ti) 1000.) 2 4)
   " secs."
  )
  )
  (princ)
) 

;;*********************************************************;
;;          ;
;; Written by  ElpanovEvgeniy       ;
;; 17.10.2008         ;
;; Calculation of the centre of a circle and circle radius ;
;; for program triangulate       ;
;;          ;
;; Modified ymg june 2011 (renamed variables)     ;
;;*********************************************************;
(defun getcircumcircle (a el /  b c c2 cp r ang)
     
  (setq b (car el)
c (cadr el)
c2 (list (car c) (cadr c)) ;c2 is point c but in 2d
  ) 
  (if (not
  (zerop
     (setq ang (- (angle b c) (angle b a)))
  )
      )
     (progn (setq cp (polar c2
    (+ -1.570796326794896 (angle c a) ang)
    (setq r (/ (distance a c2) (sin ang) 2.0))
       )
     r (abs r)
     ) 
     (list (+ (car cp) r) cp r a b c)
     ) 
  ) 
)
(defun c:gpsareasph ( / osm plt file f csvrow Latt Long gpsdata ER ptx pty ptz pt ptlst sph1 sph2 sph ptt plptlst pl pln ss dst +dp -dp dd n ent entt sss ar )
 (vl-cmdf "_.UCS" "w")
 (if (ssget "_X")
   (progn 
     (alert "\nDWG must be empty - empty it and start over")
     (exit)
   )
 )
 (setq osm (getvar 'osmode))
 (setq plt (getvar 'plinetype))
 (setvar 'osmode 0)
 (setvar 'plinetype 2)
 (setq file (getfiled "Select file with GPS coordiantes" "" "csv" 2))
 (setq f (open file "r"))
 (while (setq csvrow (read-line f))
   (setq Latt (substr csvrow 1 (vl-string-position (ascii ",") csvrow)))
   (setq Long (substr csvrow (+ (vl-string-position (ascii ",") csvrow) 2) (- (strlen csvrow) (vl-string-position (ascii ",") csvrow) 1)))
   (setq gpsdata (cons (list (read Latt) (read Long)) gpsdata))
 )
 (setq gpsdata (reverse gpsdata))
 (setq ER 6371000.0) ; Radius of Earth in m 
 (foreach gps gpsdata
   (setq Latt (car gps))
   (setq Long (cadr gps))
   (setq ptx (* (- ER 1.0) (cos (cvunit Latt "degrees" "radians")) (cos (cvunit Long "degrees" "radians"))))
   (setq pty (* (- ER 1.0) (cos (cvunit Latt "degrees" "radians")) (sin (cvunit Long "degrees" "radians"))))
   (setq ptz (* (- ER 1.0) (sin (cvunit Latt "degrees" "radians"))))
   (setq pt (list ptx pty ptz))
   (setq ptlst (cons pt ptlst))
 )
 (setq ptlst (reverse ptlst))
 (vl-cmdf "_.SPHERE" '(0.0 0.0 0.0) ER)
 (setq sph1 (entlast))
 (vl-cmdf "_.SPHERE" '(0.0 0.0 0.0) (- ER 1.0))
 (setq sph2 (entlast))
 (vl-cmdf "_.SUBTRACT" sph1 "" sph2 "")
 (setq sph (entlast))
 (vl-cmdf "_.PLAN" "")
 (triangulate ptlst)
 (gc)
 (foreach pt ptlst
   (setq ptt (list (car pt) (cadr pt) 0.0))
   (setq plptlst (cons ptt plptlst))
 )
 (setq plptlst (reverse plptlst))
 (setq pl (entmakex (append (list '(0 . "LWPOLYLINE") '(100 . "AcDbEntity") '(100 . "AcDbPolyline")
                                     (cons 90 (length plptlst))
                                    '(70 . 1)
                            )
                            (mapcar '(lambda ( x ) (cons 10 x)) plptlst)
                            (list '(210 0.0 0.0 1.0))
                    )
          )
 )
 (initget 2)
 (setq dd (getreal "\nType offset distance for determining contour shape of GPS area - if not correct enter (-10.0) <10.0> : "))
 (if (null dd) (setq dd 10.0))
 (setq pln (vlax-vla-object->ename (car (vlax-safearray->list (vlax-variant-value (vla-offset (vlax-ename->vla-object pl) dd))))))
 (setq ptlst (mapcar 'cdr (vl-remove-if-not '(lambda ( a ) (eq (car a) 10)) (entget pln))))
 (entdel pl)
 (entdel pln)
 (setq ss (ssget "_CP" (cons (car ptlst) (reverse ptlst)) '((0 . "LWPOLYLINE"))))
 (vl-cmdf "_.REGION" ss "")
 (setq ss (ssget "_X" '((0 . "REGION"))))
 (setq dst 100000.0)
 (setq +dp (/ dst 2.0))
 (setq -dp (- +dp))
 (setq sss (ssadd))
 (repeat (setq n (sslength ss))
   (setq ent (ssname ss (setq n (1- n))))
   (setq entt (vlax-vla-object->ename (vla-copy (vlax-ename->vla-object ent))))
   (vl-cmdf "_.EXTRUDE" ent "" +dp)
   (setq ent (entlast))
   (vl-cmdf "_.EXTRUDE" entt "" -dp)
   (setq entt (entlast))
   (vl-cmdf "_.UNION" ent entt "")
   (ssadd (entlast) sss)
 )
 (vl-cmdf "_.UNION" sss "")
 (vl-cmdf "_.INTERSECT" sph (entlast) "")
 (vl-cmdf "_.ZOOM" "_E")
 (vl-cmdf "_.ERASE" (ssget "_X" '((0 . "LWPOLYLINE"))) "")
 (setq ar (vla-get-volume (vlax-ename->vla-object (entlast))))
 (prompt "\nArea between GPS coordinates is : ") (princ ar) (prompt " square meters")
 (setvar 'osmode osm)
 (setvar 'plinetype plt)
 (princ)
)

home road csv.mr .txt

gpspronsph.lsp

gpsareasph.lsp

Although not so precise, but I did get data ab area... It's around 75 ha...

So try, this - my second code :

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;                                                                                 ;;;
;;; gpsareasph function takes GPS coordinates from CSV file,                        ;;;
;;; and project between points edges as arcs on sphere with radius of planet Earth, ;;;
;;; then it asks for offset distance to determine contour shape of peace of ground, ;;;
;;; for witch if finds approx. correct area value.                                  ;;;
;;; 25.07.2012. Marko Ribar, d.i.a. - second and final release                      ;;;
;;;                                                                                 ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(vl-load-com)
(defun c:gpsareasph ( / osm file f csvrow Latt Long gpsdata ER ptx pty ptz pt ptlst sph1 sph2 sph ptt plptlst pl dst smallzpt ar )
 (vl-cmdf "_.UCS" "w")
 (if (ssget "_X")
   (progn 
     (alert "\nDWG must be empty - empty it and start over")
     (exit)
   )
 )
 (setq osm (getvar 'osmode))
 (setvar 'osmode 0)
 (setq file (getfiled "Select file with GPS coordiantes" "" "csv" 2))
 (setq f (open file "r"))
 (while (setq csvrow (read-line f))
   (setq Latt (substr csvrow 1 (vl-string-position (ascii ",") csvrow)))
   (setq Long (substr csvrow (+ (vl-string-position (ascii ",") csvrow) 2) (- (strlen csvrow) (vl-string-position (ascii ",") csvrow) 1)))
   (setq gpsdata (cons (list (read Latt) (read Long)) gpsdata))
 )
 (setq gpsdata (reverse gpsdata))
 (setq ER 6371000.0) ; Radius of Earth in m 
 (foreach gps gpsdata
   (setq Latt (car gps))
   (setq Long (cadr gps))
   (setq ptx (* (- ER 2.0) (cos (cvunit Latt "degrees" "radians")) (cos (cvunit Long "degrees" "radians"))))
   (setq pty (* (- ER 2.0) (cos (cvunit Latt "degrees" "radians")) (sin (cvunit Long "degrees" "radians"))))
   (setq ptz (* (- ER 2.0) (sin (cvunit Latt "degrees" "radians"))))
   (setq pt (list ptx pty ptz))
   (setq ptlst (cons pt ptlst))
 )
 (setq ptlst (reverse ptlst))
 (vl-cmdf "_.SPHERE" '(0.0 0.0 0.0) ER)
 (setq sph1 (entlast))
 (vl-cmdf "_.SPHERE" '(0.0 0.0 0.0) (- ER 1.0))
 (setq sph2 (entlast))
 (vl-cmdf "_.SUBTRACT" sph1 "" sph2 "")
 (setq sph (entlast))
 (vl-cmdf "_.PLAN" "")
 (foreach pt ptlst
   (setq ptt (list (car pt) (cadr pt) 0.0))
   (setq plptlst (cons ptt plptlst))
 )
 (setq plptlst (reverse plptlst))
 (setq pl (entmakex (append (list '(0 . "LWPOLYLINE") '(100 . "AcDbEntity") '(100 . "AcDbPolyline")
                                     (cons 90 (length plptlst))
                                    '(70 . 1)
                            )
                            (mapcar '(lambda ( x ) (cons 10 x)) plptlst)
                            (list '(210 0.0 0.0 1.0))
                    )
          )
 )
 (vl-cmdf "_.REGION" pl "")
 (setq dst 100000.0)
 (vl-cmdf "_.EXTRUDE" (entlast) "" dst)
 (setq smallzpt (car (vl-sort ptlst '(lambda ( a b ) (< (caddr a) (caddr b))))))
 (vl-cmdf "_.MOVE" (entlast) "" (list (car smallzpt) (cadr smallzpt) 0.0) smallzpt)
 (vl-cmdf "_.INTERSECT" sph (entlast) "")
 (vl-cmdf "_.ZOOM" "_E")
 (setq ar (vla-get-volume (vlax-ename->vla-object (entlast))))
 (prompt "\nArea between GPS coordinates is : ") (princ ar) (prompt " square meters")
 (setvar 'osmode osm)
 (princ)
)

M.R.

gpsareasph-new.lsp

This meight work, but in this case with these GPS coordinates on my comp it doesn't make region objects and so it can't complete operation to the end as predicted... But I think this is the best solution, for it is the most precise in calculating area... Still if GPS coordinates are near equator, maybe it can't make solid with SURFSCULPT - plane and faces don't interfere each other... This is only for A2012 and + (uses SURFSCULPT command) :

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;                                                                                 ;;;
;;; gpsareasph function takes GPS coordinates from CSV file,                        ;;;
;;; and project between points edges as arcs on sphere with radius of planet Earth, ;;;
;;; determines contour shape of peace of ground,                                    ;;;
;;; for witch if finds approx. correct area value.                                  ;;;
;;; 30.07.2012. Marko Ribar, d.i.a. - third and final release                       ;;;
;;;                                                                                 ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(vl-load-com)
(defun c:gpsareasph ( / osm file f csvrow Latt Long gpsdata ER ptx pty ptz pt ptlst sph1 sph2 sph k ss p1 p2 ar )
 (vl-cmdf "_.UCS" "w")
 (if (ssget "_X")
   (progn 
     (alert "\nDWG must be empty - empty it and start over")
     (exit)
   )
 )
 (setq osm (getvar 'osmode))
 (setvar 'osmode 0)
 (setq file (getfiled "Select file with GPS coordiantes" "" "csv" 2))
 (setq f (open file "r"))
 (while (setq csvrow (read-line f))
   (setq Latt (substr csvrow 1 (vl-string-position (ascii ",") csvrow)))
   (setq Long (substr csvrow (+ (vl-string-position (ascii ",") csvrow) 2) (- (strlen csvrow) (vl-string-position (ascii ",") csvrow) 1)))
   (setq gpsdata (cons (list (read Latt) (read Long)) gpsdata))
 )
 (setq gpsdata (reverse gpsdata))
 (setq ER 6371000.0) ; Radius of Earth in m 
 (foreach gps gpsdata
   (setq Latt (car gps))
   (setq Long (cadr gps))
   (setq ptx (* (* ER 5.0) (cos (cvunit Latt "degrees" "radians")) (cos (cvunit Long "degrees" "radians"))))
   (setq pty (* (* ER 5.0) (cos (cvunit Latt "degrees" "radians")) (sin (cvunit Long "degrees" "radians"))))
   (setq ptz (* (* ER 5.0) (sin (cvunit Latt "degrees" "radians"))))
   (setq pt (list ptx pty ptz))
   (setq ptlst (cons pt ptlst))
 )
 (setq ptlst (reverse ptlst))
 (vl-cmdf "_.SPHERE" '(0.0 0.0 0.0) ER)
 (setq sph1 (entlast))
 (vl-cmdf "_.SPHERE" '(0.0 0.0 0.0) (- ER 1.0))
 (setq sph2 (entlast))
 (vl-cmdf "_.SUBTRACT" sph1 "" sph2 "")
 (setq sph (entlast))
 (vl-cmdf "_.PLAN" "")
 (setq k -1)
 (setq ss (ssadd))
 (repeat (length ptlst)
   (setq p1 (nth (setq k (1+ k)) ptlst))
   (if (= k (- (length ptlst) 1)) (setq p2 (nth 0 ptlst)) (setq p2 (nth (1+ k) ptlst)))
   (vl-cmdf "_.UCS" "3" '(0.0 0.0 0.0) p1 p2)
   (vl-cmdf "_.PLINE" '(0.0 0.0 0.0) (trans p1 0 1) (trans p2 0 1) "c")
   (vl-cmdf "_.REGION" (entlast) "")
   (ssadd (entlast) ss)
   (vl-cmdf "_.UCS" "p")
 )
 (vl-cmdf "_.ZOOM" "_E")
 (vl-cmdf "_.RECTANGLE" (list (- (* ER 50.0)) (- (* ER 50.0)) (+ ER 1.0)) (list (+ (* ER 50.0)) (+ (* ER 50.0)) (+ ER 1.0)))
 (vl-cmdf "_.REGION" (entlast) "")
 (ssadd (entlast) ss)
 (vl-cmdf "_.SURFSCULPT" ss "")
 (vl-cmdf "_.INTERSECT" sph (entlast) "")
 (vl-cmdf "_.ZOOM" "_E")
 (setq ar (vla-get-volume (vlax-ename->vla-object (entlast))))
 (prompt "\nArea between GPS coordinates is : ") (princ ar) (prompt " square meters")
 (setvar 'osmode osm)
 (princ)
)

M.R.

gpsareasph-new-new.lsp

Hi Marko,

Thanks for your codes, I have tried them all, and I will use "gpsareasph-new.lsp‎" since I have just AutoCad 2009 and "gpsareasph-new-new.lsp‎" doesn't work, but I will keep it maybe in one day I will have 2012.

Thanks again,

Cristian

Although not so precise, but I did get data ab area... It's around 75 ha...

 

So try, this - my second code :

 

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; ;;;
;;; gpsareasph function takes GPS coordinates from CSV file, ;;;
;;; and project between points edges as arcs on sphere with radius of planet Earth, ;;;
;;; then it asks for offset distance to determine contour shape of peace of ground, ;;;
;;; for witch if finds approx. correct area value. ;;;
;;; 25.07.2012. Marko Ribar, d.i.a. - second and final release ;;;
;;; ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(vl-load-com)
(defun c:gpsareasph ( / osm file f csvrow Latt Long gpsdata ER ptx pty ptz pt ptlst sph1 sph2 sph ptt plptlst pl dst smallzpt ar )
(vl-cmdf "_.UCS" "w")
(if (ssget "_X")
(progn 
(alert "\nDWG must be empty - empty it and start over")
(exit)
)
)
(setq osm (getvar 'osmode))
(setvar 'osmode 0)
(setq file (getfiled "Select file with GPS coordiantes" "" "csv" 2))
(setq f (open file "r"))
(while (setq csvrow (read-line f))
(setq Latt (substr csvrow 1 (vl-string-position (ascii ",") csvrow)))
(setq Long (substr csvrow (+ (vl-string-position (ascii ",") csvrow) 2) (- (strlen csvrow) (vl-string-position (ascii ",") csvrow) 1)))
(setq gpsdata (cons (list (read Latt) (read Long)) gpsdata))
)
(setq gpsdata (reverse gpsdata))
(setq ER 6371000.0) ; Radius of Earth in m 
(foreach gps gpsdata
(setq Latt (car gps))
(setq Long (cadr gps))
(setq ptx (* (- ER 2.0) (cos (cvunit Latt "degrees" "radians")) (cos (cvunit Long "degrees" "radians"))))
(setq pty (* (- ER 2.0) (cos (cvunit Latt "degrees" "radians")) (sin (cvunit Long "degrees" "radians"))))
(setq ptz (* (- ER 2.0) (sin (cvunit Latt "degrees" "radians"))))
(setq pt (list ptx pty ptz))
(setq ptlst (cons pt ptlst))
)
(setq ptlst (reverse ptlst))
(vl-cmdf "_.SPHERE" '(0.0 0.0 0.0) ER)
(setq sph1 (entlast))
(vl-cmdf "_.SPHERE" '(0.0 0.0 0.0) (- ER 1.0))
(setq sph2 (entlast))
(vl-cmdf "_.SUBTRACT" sph1 "" sph2 "")
(setq sph (entlast))
(vl-cmdf "_.PLAN" "")
(foreach pt ptlst
(setq ptt (list (car pt) (cadr pt) 0.0))
(setq plptlst (cons ptt plptlst))
)
(setq plptlst (reverse plptlst))
(setq pl (entmakex (append (list '(0 . "LWPOLYLINE") '(100 . "AcDbEntity") '(100 . "AcDbPolyline")
(cons 90 (length plptlst))
'(70 . 1)
)
(mapcar '(lambda ( x ) (cons 10 x)) plptlst)
(list '(210 0.0 0.0 1.0))
)
)
)
(vl-cmdf "_.REGION" pl "")
(setq dst 100000.0)
(vl-cmdf "_.EXTRUDE" (entlast) "" dst)
(setq smallzpt (car (vl-sort ptlst '(lambda ( a b ) (< (caddr a) (caddr b))))))
(vl-cmdf "_.MOVE" (entlast) "" (list (car smallzpt) (cadr smallzpt) 0.0) smallzpt)
(vl-cmdf "_.INTERSECT" sph (entlast) "")
(vl-cmdf "_.ZOOM" "_E")
(setq ar (vla-get-volume (vlax-ename->vla-object (entlast))))
(prompt "\nArea between GPS coordinates is : ") (princ ar) (prompt " square meters")
(setvar 'osmode osm)
(princ)
)

M.R.

Thanks Marko,

I have tried all your codes, but this one works the best for me, since I have just Autocad 2009, and as you said, your last code does not work for me yet, maybe in one day when I will have 2012 version.

Thanks again,

Cristian

Thanks Cristian, I've modified further more 3rd code - better finding solid to intersect with on north hemisphere (SURFSCULPT is between regions and north dome surface)... Still with your GPS coordinates I can't convert PLINES into REGIONS, but with some other (larger angles) it does the job...

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;                                                                                 ;;;
;;; gpsareasph function takes GPS coordinates from CSV file,                        ;;;
;;; and project between points edges as arcs on sphere with radius of planet Earth, ;;;
;;; determines contour shape of peace of ground,                                    ;;;
;;; for witch if finds approx. correct area value.                                  ;;;
;;; 30.07.2012. Marko Ribar, d.i.a. - fourth and final release                      ;;;
;;; *** use this code for large distances between GPS coordinates ***               ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(vl-load-com)
(defun c:gpsareasph ( / osm file f csvrow Latt Long gpsdata ER ptx pty ptz pt ptlst sph1 sph2 sph k ss p1 p2 ar )
 (vl-cmdf "_.UCS" "w")
 (if (ssget "_X")
   (progn 
     (alert "\nDWG must be empty - empty it and start over")
     (exit)
   )
 )
 (setq osm (getvar 'osmode))
 (setvar 'osmode 0)
 (setq file (getfiled "Select file with GPS coordiantes" "" "csv" 2))
 (setq f (open file "r"))
 (while (setq csvrow (read-line f))
   (setq Latt (substr csvrow 1 (vl-string-position (ascii ",") csvrow)))
   (setq Long (substr csvrow (+ (vl-string-position (ascii ",") csvrow) 2) (- (strlen csvrow) (vl-string-position (ascii ",") csvrow) 1)))
   (setq gpsdata (cons (list (read Latt) (read Long)) gpsdata))
 )
 (setq gpsdata (reverse gpsdata))
 (setq ER 6371000.0) ; Radius of Earth in m 
 (foreach gps gpsdata
   (setq Latt (car gps))
   (setq Long (cadr gps))
   (setq ptx (* (* ER 2.0) (cos (cvunit Latt "degrees" "radians")) (cos (cvunit Long "degrees" "radians"))))
   (setq pty (* (* ER 2.0) (cos (cvunit Latt "degrees" "radians")) (sin (cvunit Long "degrees" "radians"))))
   (setq ptz (* (* ER 2.0) (sin (cvunit Latt "degrees" "radians"))))
   (setq pt (list ptx pty ptz))
   (setq ptlst (cons pt ptlst))
 )
 (setq ptlst (reverse ptlst))
 (vl-cmdf "_.SPHERE" '(0.0 0.0 0.0) ER)
 (setq sph1 (entlast))
 (vl-cmdf "_.SPHERE" '(0.0 0.0 0.0) (- ER 1e-4))
 (setq sph2 (entlast))
 (vl-cmdf "_.SUBTRACT" sph1 "" sph2 "")
 (setq sph (entlast))
 (vl-cmdf "_.PLAN" "")
 (setq k -1)
 (setq ss (ssadd))
 (repeat (length ptlst)
   (setq p1 (nth (setq k (1+ k)) ptlst))
   (if (= k (- (length ptlst) 1)) (setq p2 (nth 0 ptlst)) (setq p2 (nth (1+ k) ptlst)))
   (vl-cmdf "_.UCS" "3" '(0.0 0.0 0.0) p1 p2)
   (vl-cmdf "_.PLINE" '(0.0 0.0 0.0) (trans p1 0 1) (trans p2 0 1) "c")
   (vl-cmdf "_.REGION" (entlast) "")
   (ssadd (entlast) ss)
   (vl-cmdf "_.UCS" "_P")
 )
 (vl-cmdf "_.ZOOM" "_E")
 (vl-cmdf "_.SPHERE" '(0.0 0.0 0.0) (* ER 1.25))
 (vl-cmdf "_.CONVTOSURFACE" (entlast) "")
 (vl-cmdf "_.SLICE" (entlast) "" "XY" "" "")
 (entdel (entlast))
 (ssadd (entlast) ss)
 (vl-cmdf "_.SURFSCULPT" ss "")
 (vl-cmdf "_.INTERSECT" sph (entlast) "")
 (vl-cmdf "_.ZOOM" "_E")
 (vl-cmdf "_.AREA" "o" (entlast) "")
 (setq ar (getvar 'area))
 (setq ar (/ ar 2.0)
 (prompt "\nArea between GPS coordinates is : ") (princ ar) (prompt " square meters")
 (setvar 'osmode osm)
 (princ)
)

M.R.

So maybe when you get A2012, you'll look into this way...

gpsareasph-new-new-new.lsp

*** PROBLEM SOLVED ***

This is exactly what you should get :

Area between GPS coordinates is : 741875.0 square meters

Here is my last and the most efficient code - of course only for one Earth hemisphere and works I think for A2011 and +... Tested on A2012...

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;                                                                                 ;;;
;;; gpsareasph function takes GPS coordinates from CSV file,                        ;;;
;;; and project between points edges as arcs on sphere with radius of planet Earth, ;;;
;;; determines contour shape of peace of ground,                                    ;;;
;;; for witch if finds approx. correct area value.                                  ;;;
;;; 31.07.2012. Marko Ribar, d.i.a. - fifth and final release                       ;;;
;;; *** improved possibility of making REGION objects ***                           ;;;
;;; *** use this code for small distances between GPS coordinates ***               ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(vl-load-com)
(defun c:gpsareasph ( / osm file f csvrow Latt Long gpsdata ER pt1x pt1y pt1z pt1 pt1lst pt2x pt2y pt2z pt2 pt2lst sph1 sph2 sph k ss p1 p2 p3 p4 ar )
 (vl-cmdf "_.UCS" "w")
 (if (ssget "_X")
   (progn 
     (alert "\nDWG must be empty - empty it and start over")
     (exit)
   )
 )
 (setq osm (getvar 'osmode))
 (setvar 'osmode 0)
 (setq file (getfiled "Select file with GPS coordiantes" "" "csv" 2))
 (setq f (open file "r"))
 (while (setq csvrow (read-line f))
   (setq Latt (substr csvrow 1 (vl-string-position (ascii ",") csvrow)))
   (setq Long (substr csvrow (+ (vl-string-position (ascii ",") csvrow) 2) (- (strlen csvrow) (vl-string-position (ascii ",") csvrow) 1)))
   (setq gpsdata (cons (list (read Latt) (read Long)) gpsdata))
 )
 (setq gpsdata (reverse gpsdata))
 (setq ER 6371000.0) ; Radius of Earth in m 
 (foreach gps gpsdata
   (setq Latt (car gps))
   (setq Long (cadr gps))
   (setq pt1x (* (* ER 1.02) (cos (cvunit Latt "degrees" "radians")) (cos (cvunit Long "degrees" "radians"))))
   (setq pt1y (* (* ER 1.02) (cos (cvunit Latt "degrees" "radians")) (sin (cvunit Long "degrees" "radians"))))
   (setq pt1z (* (* ER 1.02) (sin (cvunit Latt "degrees" "radians"))))
   (setq pt1 (list pt1x pt1y pt1z))
   (setq pt1lst (cons pt1 pt1lst))
   (setq pt2x (* (* ER 0.98) (cos (cvunit Latt "degrees" "radians")) (cos (cvunit Long "degrees" "radians"))))
   (setq pt2y (* (* ER 0.98) (cos (cvunit Latt "degrees" "radians")) (sin (cvunit Long "degrees" "radians"))))
   (setq pt2z (* (* ER 0.98) (sin (cvunit Latt "degrees" "radians"))))
   (setq pt2 (list pt2x pt2y pt2z))
   (setq pt2lst (cons pt2 pt2lst))
 )
 (setq pt1lst (reverse pt1lst))
 (setq pt2lst (reverse pt2lst))
 (vl-cmdf "_.SPHERE" '(0.0 0.0 0.0) ER)
 (setq sph1 (entlast))
 (vl-cmdf "_.SPHERE" '(0.0 0.0 0.0) (- ER 1e-4))
 (setq sph2 (entlast))
 (vl-cmdf "_.SUBTRACT" sph1 "" sph2 "")
 (setq sph (entlast))
 (vl-cmdf "_.PLAN" "")
 (setq k -1)
 (setq ss (ssadd))
 (repeat (length pt1lst)
   (setq p1 (nth (setq k (1+ k)) pt1lst))
   (if (= k (- (length pt1lst) 1)) (setq p2 (nth 0 pt1lst)) (setq p2 (nth (1+ k) pt1lst)))
   (setq p3 (nth k pt2lst))
   (if (= k (- (length pt2lst) 1)) (setq p4 (nth 0 pt2lst)) (setq p4 (nth (1+ k) pt2lst)))
   (vl-cmdf "_.UCS" "3" p1 p2 p3)
   (vl-cmdf "_.PLINE" (trans p1 0 1) (trans p2 0 1) (trans p4 0 1) (trans p3 0 1) "c")
   (vl-cmdf "_.REGION" (entlast) "")
   (ssadd (entlast) ss)
   (vl-cmdf "_.UCS" "_P")
 )
 (vl-cmdf "_.SPHERE" '(0.0 0.0 0.0) (* ER 1.01))
 (vl-cmdf "_.CONVTOSURFACE" (entlast) "")
 (vl-cmdf "_.SLICE" (entlast) "" "XY" "" "")
 (entdel (entlast))
 (ssadd (entlast) ss)
 (vl-cmdf "_.SPHERE" '(0.0 0.0 0.0) (* ER 0.99))
 (vl-cmdf "_.CONVTOSURFACE" (entlast) "")
 (vl-cmdf "_.SLICE" (entlast) "" "XY" "" "")
 (entdel (entlast))
 (ssadd (entlast) ss)
 (vl-cmdf "_.SURFSCULPT" ss "")
 (vl-cmdf "_.INTERSECT" sph (entlast) "")
 (vl-cmdf "_.ZOOM" "_E")
 (vl-cmdf "_.AREA" "o" (entlast) "")
 (setq ar (getvar 'area))
 (setq ar (/ ar 2.0))
 (prompt "\nArea between GPS coordinates is : ") (princ ar) (prompt " square meters")
 (setvar 'osmode osm)
 (princ)
)

Many regards,

M.R.

gpsareasph-new-new-new-new.lsp

Thank you very much, your lisp's are briliant, i am writing lisp's as well, but mine are just at beginer level, I am still learning, but your lisp's are way above me, for my own curiosity, why did you had to use a "z" coordinate? is this how the formula works? or you really need a "z"?

I have checked your solid, and the calculation is perfect, from my phone I can get only the perimeter which is: 5888.34m.

Best whishes,

Cristian

Cristian, thank you for your comments, however I have to point that lsp for drawing arcs is OK... You can always use it to get lengths - distances... But lsp for getting area is difficult to predict weather angles between 2 near GPS coordinates are going to be large or small, so I suggest A2011 or + and using both of mine codes - fourth for large angles, and fifth for small as in your case... I think that combine them into sinlge code isn't soo good solution, but there maybe cases where you have both large angles and also small - in this case you have to find middle solution and only prey that ACAD will convert plines to regions...

M.R.