Lisp to draw a line between GPS points

[Cristi_an_24]

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 (- 1 a)))))

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

[marko_ribar]

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)
)

[marko_ribar]

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.

[marko_ribar]

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.

[Cristi_an_24]

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

[Cristi_an_24]

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

[marko_ribar]

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...

[marko_ribar]

*** 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.

[Cristi_an_24]

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


Screenshot_2012-07-31-08-06-02.jpg

[marko_ribar]

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.