Surface Area of regions on an ellipsoid Earth

(defun spherecutarea (latdeg sidelength deglen)
  "Calculate the area of a cut in a sphere at the latitude LATDEG
with the given SIDELENGTH and the length of one degree at the
Equator DEGLEN."
  (* deglen sidelength ; longitude 
     deglen sidelength (cos (* 2 float-pi (/ latdeg 360.0))))) ; latitude

(defun spherearea (latdeg sidelength)
  "Area of a segment of a sphere at LATDEG with the given
SIDELENGTH."
  (let* ((R 6371.0) ; km
         (deglen (/ (* 2 float-pi R) 360.0))) ; km^2
    (spherecutarea latdeg sidelength deglen)))

(defun spheresegmentarea (latdeg sidelength)
  "Calculate the area of a rectangular segment on a sphere at
latitude LATDEG with the given SIDELENGTH."
  (* 24728.6234228 sidelength sidelength 
     (cos (* float-pi (/ latdeg 180.0)))))

(defmacro turntofloatsingle (var)
  (list 'setq var (list 'float var)))

(defmacro turntofloat (&rest vars)
  "Turn a list of items to floats."
  (cons 'progn (mapcar 
                (lambda (var) 
                  (list 'turntofloatsingle var))
                vars)))

(defun ellipsoiddistance (a f L1 L2 φ1 φ2)
  "Calculate the distance of two arbitrary points on an ellipsoid.

  Parameters: Equator radius A, oblateness F and for point 1 and
  2 respectively the longitudes L1 and L2 and the latitudes φ1
  and φ2."
  ; ensure that we work on floats
  (turntofloat a f φ1 φ2 L1 L2)
  ; the first simplifications don’t depend on each other, 
  ; so we just use let to bind them
  (let ((F (/ (+ φ1 φ2) 2))
        (G (/ (- φ1 φ2) 2))
        (λ (/ (- L1 L2) 2)))
    (message (format "F %f G %f λ %f a %f f %f L1 %f L2 %f φ1 %f φ2 %f" 
                     F G λ a f L1 L2 φ1 φ2))
    ; the second don’t depend on each other either
    (let ((S (+ (* (expt (sin G) 2)
                   (expt (cos λ) 2))
                (* (expt (cos F) 2)
                   (expt (sin λ) 2))))
          (C (+ (* (expt (cos G) 2)
                   (expt (cos λ) 2))        
                (* (expt (sin F) 2)
                   (expt (sin λ) 2)))))
      ; now we have a few consecutive definitions, so we use let*
      ; which allows references to previous elements in the same let*.
      (let* ((ω (atan (sqrt (/ S C))))
             (R (/    (sqrt (* S C)) ω)))
        (let ((D (* 2 ω a))
              (H1 (/ (- (* 3 R)) (* 2 C)))
              (H2 (/ (+ (* 3 R)) (* 2 C))))
          ; All prepared. Now we just fit all this together. This is
          ; the last line, so the function returns the value.
          (* D (- 
                (+ 1 (* f H1 (expt (sin F) 2) (expt (cos G) 2)))
                (* f H2 (expt (cos F) 2) (expt (sin G) 2)))))))))

(defun ellipsoidrectanglearea (a f longitude latitude dlon dlat)
  (let ((L1 longitude)
        (L2 (+ longitude dlon))
        (φ1 latitude)
        (φ2 (+ latitude dlat)))
    (let ((lenlower (ellipsoiddistance a f L1 L2 φ1 φ1))
          (lenupper (ellipsoiddistance a f L1 L2 φ2 φ2))
          (lenwestern (ellipsoiddistance a f L1 L1 φ1 φ2))
          (leneastern (ellipsoiddistance a f L2 L2 φ1 φ2)))
      (if (not (= lenwestern leneastern))
          (error "Western and Eastern length are not equal. 
This violates the laws of geometry. We die. Western: %f Eastern: %f" 
                 lenwestern leneastern))
      (let ((horizontalmean (/ (+ lenlower lenupper) 2)))
        ; now just return length times width
        (* horizontalmean lenwestern)))))

<<ellipsoid-helpers>>
<<ellipsoid-distance>>
<<ellipsoid-rectanglearea>>

(defun ellipsoidrectangleareafromdeg (latdeg sidelength)
  "Calculate the rectangle area from the latitude LATDEG and the
SIDELENGTH given as degrees."
  (message (format "latdeg %f sidelength %f" latdeg sidelength))
  (let ((londeg 15) ; irrelevant due to symmetry
        (dlondeg sidelength)
        (dlatdeg sidelength)
        (a 6378.139)
        (f (/ 1 298.25642)))
    (let ((lon (* 2 float-pi (/ londeg 360.0))) ; 2π / 360
          (dlon (* 2 float-pi (/ dlondeg 360.0)))
          (lat (* 2 float-pi (/ latdeg 360.0))) 
          (dlat (* 2 float-pi (/ dlatdeg 360.0))))
      (ellipsoidrectanglearea a f lon lat dlon dlat))))

(defun ellipsoidrectangleareafromdegnumericalintegration (latdeg sidelength steps)
  "Calculate the rectangle area from the latidute LATDEG and the
  SIDELENGTH given as degrees by adding them in STEPS smaller steps per sidelength."
  (let ((area 0)
        (smallerside (/ (float sidelength)
                        (float steps))))
    (loop for i from 0 to (1- steps) by 1 do
          (message (format "i %f" i))
          (let ((smallerlat (+ latdeg (* smallerside i))))
            ; add stepts times the area since the longitudinal
            ; calculation does not change, so we only need to
            ; calculate in once.
            (setq area (+ area (* steps 
                                  (ellipsoidrectangleareafromdeg
                                   smallerlat smallerside))))))

    area))
; no return value
nil

# add map lines
m.drawcoastlines()
m.drawparallels(np.arange(-90.,120.,30.), 
                labels=[False,True,True,False])
m.drawmeridians(np.arange(0.,420.,60.), 
                labels=[True,False,False,True])
m.drawmapboundary(fill_color='aqua')

import numpy as np
import pylab as pl
import mpl_toolkits.basemap as bm
import netCDF4 as nc
def singlehemispherelats2map(northernlats):
    """Turn the northern lats (0-90) into a map (180,360)."""
    # duplicate the northernlats
    lats = np.zeros((180, ))
    lats[0:90] = northernlats[:0:-1,1]
    lats[90:] = northernlats[1:,1]
    # and blow them up into a map
    lons = np.ones((360, ))
    lats = np.matrix(lats)
    lons = np.matrix(lons)
    mapscaling = lons.transpose() * lats
    mapscaling = mapscaling.transpose()
    return mapscaling

# first read the file
with open("transcomellipticlat90-sum1000.dat") as f:
    northernlats = np.genfromtxt(f, delimiter=" ")
mapscaling = singlehemispherelats2map(northernlats)  
with open("transcomellipticlat90-sum1000vsdirect.dat") as f:
    northernlats = np.genfromtxt(f, delimiter=" ")
mapscalingdiff = singlehemispherelats2map(northernlats)
with open("transcomellipticlat90-direct.dat") as f:
    northernlats = np.genfromtxt(f, delimiter=" ")
mapscalingdirect = singlehemispherelats2map(northernlats)
with open("transcomellipticlat90-sphere.dat") as f:
    northernlats = np.genfromtxt(f, delimiter=" ")
mapscalingsphere = singlehemispherelats2map(northernlats)
with open("transcomellipticlat90-sum1000vssphere.dat") as f:
    northernlats = np.genfromtxt(f, delimiter=" ")
mapscalingdiffsphere = singlehemispherelats2map(northernlats)

# several different plots:
<<plotareamapperpixel>>
<<plotareamapperpixeldirect>>
<<plotareamapperpixelerror>>
<<plotareamapperpixelrelerror>>
<<plotareamapperpixelsphereerror>>

# plot it for representation
m = bm.Basemap()
m.imshow(mapscaling)
bar = pl.colorbar()
bar.set_label("area per pixel [$km^2$]")
<<addplotstyle>>
pl.title("Surface Area 1x1 [$km^2$]")
pl.savefig("eartharea_1x1.png")
pl.close()
print """\n\n#+caption:Area when summing 1000x1000 smaller areas
[[./eartharea_1x1.png]]"""

m = bm.Basemap()
m.imshow(mapscaling)
bar = pl.colorbar()
bar.set_label("area per pixel [$km^2$]")
# summon map style! :)
<<addplotstyle>>
pl.title("Surface Area 1x1, no numerical integration [$km^2$]")
pl.savefig("earthareadirect_1x1.png")
pl.close()
print "\n\n#+caption:Area when using just one square\n[[./earthareadirect_1x1.png]]"

m = bm.Basemap()
m.imshow(mapscalingdiff)
<<addplotstyle>>
bar = pl.colorbar()
bar.set_label("area per pixel [$km^2$]")
pl.title("Surface Area 1x1 difference: sum 1000 vs direct [$km^2$]")
pl.savefig("eartharea1000vs1_1x1.png")  # save as a clean netCDF4 file
pl.close()
print "\n\n#+caption:Difference between summing smaller squares", 
print "and just using one square\n[[./eartharea1000vs1_1x1.png]]"

m = bm.Basemap()
m.imshow(np.log(np.abs(mapscalingdiff/mapscaling)))
<<addplotstyle>>
bar = pl.colorbar()
bar.set_label("relative error per pixel, logarithmic")
pl.title("Surface Area 1x1 diff relative: sum 1000 vs direct")
pl.savefig("eartharea1000vs1rel_1x1.png")  # save as a clean netCDF4 file
pl.close()
print """\n\n#+caption:Relative Area Error by not integrating (logscale)
[[./eartharea1000vs1rel_1x1.png]]"""

m = bm.Basemap()
m.imshow(np.log(np.abs(mapscalingdiffsphere/mapscaling)))
<<addplotstyle>>
bar = pl.colorbar()
bar.set_label("relative error per pixel, logarithmic")
pl.title("Surface Area 1x1 diff relative: sum 1000 vs sphere")
pl.savefig("eartharea1000vssphererel_1x1.png")
pl.close()
print """\n\n#+caption:Relative Error from Sphere (logscale)
[[./eartharea1000vssphererel_1x1.png]]"""

<<readcsvareafiles>>  
<<plotareamaps>>
D = nc.Dataset("eartharea.nc", "w")
D.comment = "Created with tm5tools/ct2pyshell/transcomareas.org"
D.createDimension("longitude", 360)
D.createDimension("latitude", 180)
area = D.createVariable("1x1", "f8", ("latitude", "longitude"))
area.units = "km^2"
area.comment = "from 180W to 180E and from 90S to 90N"
area[:] = mapscaling
area = D.createVariable("1x1_1000vs1", "f8", ("latitude", "longitude"))
area.units = "km^2"
area.comment = ("Difference between the direct calculation of the "
"area and summing up 1000x1000 smaller areas."
"from 180W to 180E and from 90S to 90N")
area[:] = mapscalingdiff
area = D.createVariable("1x1_direct", "f8", ("latitude", "longitude"))
area.units = "km^2"
area.comment = ("Area calculated without numerical intergration (bigger errors!). "
"from 180W to 180E and from 90S to 90N")
area[:] = mapscalingdirect
area = D.createVariable("1x1_sphere", "f8", ("latitude", "longitude"))
area.units = "km^2"
area.comment = ("Area calculated on a simple sphere. "
"from 180W to 180E and from 90S to 90N")
area[:] = mapscalingsphere

(defmacro turntofloatsingle (var)
  (list 'setq var (list 'float var)))

<<turntofloat-single>>
(defmacro turntofloatbackticks (&rest vars)
  "Turn a list of items to floats using backtick notation."
  `(progn ,@(mapcar 
             (lambda (var) 
               `(turntofloatsingle ,var)) 
             vars)))

<<turntofloat-single>>
(defmacro turntofloat (&rest vars)
  "Turn a list of items to floats (without using backticks)."
  ; cons turns this into a call of progn on the list returned by
  ; mapcar
  (cons 'progn (mapcar 
                (lambda (var) 
                  (list 'turntofloatsingle var))
                vars)))

<<turntofloat-single>>
(defmacro turntofloatcollect (&rest vars)
  "Turn a list of items to floats, using the collect directive of loop."
  ; execute progn on the list returned by the loop
  (cons 'progn 
        ; loop ... collect returns a list of all the loop results.
        (loop for var in vars collect 
              (list 'turntofloatsingle var))))

<<turntofloat-single>>

; build the list explicitely to make it easier for me to understand
; what the macro does
(defmacro turntofloatexplicit (&rest vars)
  "Turn a list of items to floats (using explicit list building
instead of mapcar)."
  ; prepare an empty list of function calls
  (let ((funclist '()))
    ; for each variable add a call to the single-item macro
    (loop for var in vars do
          ; (list 'turntofloatsingle var) creates the call to
          ; turntofloatsingle with the variable which is referenced by
          ; var. Push puts that at the beginning of the funclist.
          (push (list 'turntofloatsingle var) funclist))
    ; to ensure the right order of operations, we reverse the funclist
    (setq funclist (reverse funclist))
    ; cons turns this into a call of progn on the list. We need progn,
    ; because the funclist contains multiple function calls.
    (cons 'progn funclist)))

<<turntofloat-single>>
; Common Lisp Macro to turn the place to a float in one step.
(defmacro turntofloatinline (&rest places)
  "Turn a list of items to floats using an inline function call."
  `(progn ,@(mapcar 
             (lambda (place) 
               `(callf float ,place)) places)))

<<turntofloat-collect>>
(setq a 1 b 3.8 c 2)
(turntofloatcollect a b c)
(message (number-to-string c))