Creating a kernel density estimate map in R

In this post I'm going to create a kernel density estimate map in R from a file with latitude/longitude coordinates. WARNING: depending on your application the following gives incorrect results because a non-spherical kernel density estimator is used with spherical data (big thanks too Brian Rowlingson for pointing that out). I didn't find any spherical kernel density estimation implementation in R but if you know any please let me know. While exploring this issue I also wrote a similar implementation as the below code for the generation of a kernel density estimate with the spatial statistics package splancs instead of KernSmooth (splancs_kernel_density.R).

To create this map we'll use the bkde2D function in the KernSmooth package to generate the kernel density estimates and the raster package for outputting a raster file. If you don't have the KernSmooth and raster packages then they can as usually be installed with:

```install.packages("KernSmooth")
install.packages("raster")
```

In order to create the 2D binned kernel density map I first loaded the KernSmooth and raster package.

```library("KernSmooth")
library("raster")
```

Then we read the input file with the coordinates of the points that we want to generate the kernel density estimate off. This input file is an export from the Ocean Biographic Information System (OBIS) and represents distribution records from the seaweed Alaria esculenta.

```input <- "Alaria_esculenta.csv"
output <- "Alaria_esculenta_density.asc"

# get the coordinates
coordinates <- records[,2:3]
```

Next we compute the 2D binned kernel density estimate. In order to limit the size of the generated output files I set all values smaller then 0.00001 to 0.

```# compute the 2D binned kernel density estimate
est <- bkde2D(coordinates,
bandwidth=c(3,3),
gridsize=c(4320,2160),
range.x=list(c(-180,180),c(-90,90)))
est\$fhat[est\$fhat<0.00001] <- 0 ## ignore very small values
```

Finally we create a raster file from the output of bkde2D, inspect it visually and export it as an ascii file.

```# create raster
est.raster = raster(list(x=est\$x1,y=est\$x2,z=est\$fhat))
projection(est.raster) <- CRS("+init=epsg:4326")
xmin(est.raster) <- -180
xmax(est.raster) <- 180
ymin(est.raster) <- -90
ymax(est.raster) <- 90
# visually inspect the raster output
plot(est.raster)
# write an ascii file
writeRaster(est.raster,output,"ascii")
```

And the output looks like this. As you can see most records are located near 50°N and 0° which corresponds to the distribution records map generated by OBIS.
The output ascii raster file can be used as an input file for OccurrenceThinner which uses the kernel density estimate grid for filtering out records from densely sampled regions in order to avoid overfitting while building species distribution models (SDM).

Code and data for this post are in my blog github repository.

Azimuthal equidistant projection: An implementation in F#, Python and Julia

Recently I needed to calculate distance from one point to another set of points in order to find the nearest point and its distance to the origin point. I opted to project all points to the azimuthal equidistant projection. The most important property of this projection is that all distances from the center of the projection to all other points represents the true distance (assuming a spheric earth) to all other projected points. To calculate the distance you just have to calculate the euclidean distance and multiply it with the average earth radius (6371000.0 meters).

The formulas for the azimuthal equidistant projection can be found at mathworld.wolfram.com. This formula can then be directly translated into F# like this:

```open System
module AzimuthalEquidistantProjection =

let inline degToRad d = 0.0174532925199433 * d; // (1.0/180.0 * Math.PI) * d

let project centerlon centerlat lon lat =
// http://mathworld.wolfram.com/AzimuthalEquidistantProjection.html
let t1 = degToRad centerlat // latitude center of projection
let l0 = degToRad centerlon // longitude center of projection
let c = acos ((sin t1) * (sin t) + (cos t1) * (cos t) * (cos (l-l0)))
let k = c / (sin c)
let x = k * (cos t) * (sin (l-l0))
let y = k * ((cos t1) * (sin t) - (sin t1) * (cos t) * (cos (l-l0)))
(x, y)
```

I also implemented the azimuthal equidistant projection by using units of measure. I create one set of measures to distinguish between degrees and radians and another one to distinguish between x and y.

```open System
module AzimuthalEquidistantProjectionWithMeasures =
[<Measure>] type deg // degrees
[<Measure>] type x
[<Measure>] type y
type lon = float<x deg>
type lat = float<y deg>

let cos (d:float<'u rad>) = Math.Cos(float d)
let sin (d:float<'u rad>) = Math.Sin(float d)

let project (centerlon:lon) (centerlat:lat) (lon:lon) (lat:lat) =
// http://mathworld.wolfram.com/AzimuthalEquidistantProjection.html
let t1 = degToRad centerlat // latitude center of projection
let l0 = degToRad centerlon // longitude center of projection
let c = acos ((sin t1) * (sin t) + (cos t1) * (cos t) * (cos (l-l0)))
let k = c / (sin c)
let x:float<x> = k * (cos t) * (sin (l-l0))
|> LanguagePrimitives.FloatWithMeasure
let y:float<y> = k * ((cos t1) * (sin t) - (sin t1) * (cos t) * (cos (l-l0)))
|> LanguagePrimitives.FloatWithMeasure
(x, y)
```

In the past I've used the python OGR bindings and the ESRI Projection Engine for my map projection needs but this time I needed a pure python implementation so I translated the above code and optimized it a bit by precalculating some values that we'll need when we project the points during the initialization of my projection class. A similarly optimized version in F# is on fssnip.net.

```from math import cos, sin, acos, radians

class Point(object):
def __init__(self,x,y):
self.x = x
self.y = y

class AzimuthalEquidistantProjection(object):
"""
http://mathworld.wolfram.com/AzimuthalEquidistantProjection.html
"""
def __init__(self, center):
self.center = center
self.t1 = radians(center.y) ## latitude center of projection
self.l0 = radians(center.x) ## longitude center of projection
self.cost1 = cos(self.t1)
self.sint1 = sin(self.t1)

def project(self, point):
costcosll0 = cos(t) * cos(l-self.l0)
sint = sin(t)

c = acos ((self.sint1) * (sint) + (self.cost1) * costcosll0)
k = c / sin(c)
x = k * cos(t) * sin(l-self.l0)
y = k * (self.cost1 * sint - self.sint1 * costcosll0)
return Point(x, y)

import unittest
class Test_AzimuthalEquidistantProjection(unittest.TestCase):
def test_project(self):
p = AzimuthalEquidistantProjection(Point(0.0,0.0))
r = p.project(Point(1.0,1.0))
self.assertAlmostEqual(0.01745152022, r.x)
self.assertAlmostEqual(0.01745417858, r.y)

p = AzimuthalEquidistantProjection(Point(1.0,2.0))
r = p.project(Point(3.0,4.0))
self.assertAlmostEqual(0.03482860733, r.x)
self.assertAlmostEqual(0.03494898734, r.y)

p = AzimuthalEquidistantProjection(Point(-10.0001, 80.0001))
r = p.project(Point(7.935, 63.302))
self.assertAlmostEqual(0.1405127567, r.x)
self.assertAlmostEqual(-0.263406547, r.y)

if __name__ == '__main__':
unittest.main()
```

I've also implemented the straightforward projection code in Julia. Julia is a promising language built for technical computing. It's still only on version 0.2 but a lot of packages have already been created for the language.

```function degToRad(d)
d*0.0174532925199433
end

function project(centerlon, centerlat, lon, lat)