Geospatial Analytics

geospatial-analytics-ai999x

Map projections and coordinates, geodesic distance and area, spatial statistics, remote-sensing indices, solar position and terrain, and urban and network analysis.

Geospatial Question Answering Token Classification Text Classification Sentence Similarity Feature Extraction

Hugging Face ↗GitHub ↗See it on the capability orb

99capabilities

14topics

99formulas typeset

8worked examples below

1.0.0version

About this model

What Geospatial Analytics does

Geospatial Analytics (geospatial-analytics-ai999x) is an offline Python package for geospatial analytics. It covers map projections and coordinates, geodesic distance and area, spatial statistics, remote-sensing indices, solar position and terrain, and urban and network analysis.

It has 99 capabilities — 99 computational procedures — in 14 topics, the largest being Map projections, Spatial statistics, Ellipsoidal geodesy, Terrain analysis, Coordinate systems and conversions and Point pattern analysis. Ask a question in words with ask(), or call a capability by name with typed inputs with solve(). Every answer carries its value and unit, the steps that produced it and the formula in LaTeX.

When a question is not determined — a value is missing, an input is out of range, or the request is outside the package — it says so and names what it needs, instead of guessing. It runs on an ordinary CPU without a network connection (Python >=3.10), and the same question always gets the same answer.

Install

pip install geospatial-analytics-ai999x
hf download AI999X/geospatial-analytics-ai999x --local-dir geospatial-analytics-ai999x

What it covers

Cartography2Ground distance from a map distance and a representative fraction · Map scale (representative fraction) from a map and a ground distance Coordinate notation5Decimal degrees to degrees, minutes and seconds · Decode a geohash: cell centre and size · Degrees, minutes and seconds to decimal degrees …Coordinate systems and conversions7Auxiliary latitudes: geocentric, reduced, conformal, authalic and rectifying · Geocentric (ECEF) to geodetic coordinates · Geodetic to geocentric (ECEF) coordinates …Ellipsoidal geodesy9Arc length along a parallel of latitude · Area and perimeter of a geodesic polygon on the ellipsoid · Area of a latitude-longitude cell on the ellipsoid …Facility location3Closest facility for every node of a network · Maximal covering location problem (MCLP) on a network · p-median facility location on a network Geostatistics5Empirical (Matheron) semivariogram in distance bins · Fit an exponential or Gaussian variogram model (nugget, partial sill, range) · Inverse distance weighted (IDW) interpolation at a point …Map projections22Albers equal-area conic projection: coordinates and scale factors · Azimuthal equidistant projection (ellipsoidal, geodesic) · Equirectangular (equidistant cylindrical / plate carree) projection …Network accessibility6Gravity-based (Hansen) accessibility on a road network · Network connectivity indices: cyclomatic number, alpha, beta and gamma · Network detour index (circuity) between two nodes …Planar vector analysis3Discrete Frechet distance between two polylines · Douglas-Peucker line simplification · Shape compactness indices of a polygon (Polsby-Popper, Schwartzberg, convex hull, Reock)Point pattern analysis7Kernel density estimate of a point pattern at a location · Mean centre and standard distance of a point pattern · Nearest-neighbour index (Clark-Evans ratio) of a point pattern …Remote sensing6At-sensor brightness temperature from thermal radiance · Ground sample distance and swath of a nadir-looking camera · Thematic map accuracy: overall accuracy, kappa, producer's and user's accuracy …Spatial statistics10Bivariate Moran's I (spatial cross-correlation) · Distance-band spatial weights · Geary's C spatial autocorrelation with normal and randomization inference …Terrain analysis8D8 flow accumulation of a small DEM · D8 flow direction of every cell of a small DEM · Hillshade (analytical relief shading) of a 3x3 elevation window …Urban geography6Huff model: probabilities of shopping at each store from an origin (network distances) · Rank-size rule (Zipf) fit of a settlement system · Residential segregation (evenness): dissimilarity, Gini and Theil entropy indices …

Try it

Worked examples

Each question below was put to geospatial-analytics-ai999x and the answer shown is the one it returned, with the steps it took.

Ellipsoidal geodesym.ask(…)

Geodesic destination: lat1 = 52.2, lon1 = 0.12, azimuth = 45, distance = 100 km, ellipsoid = WGS84

lat2 = 52.83083647

lon2 = 1.169145951

azimuth2 = 45.83253655

Show the steps

Geodesic destination point on the ellipsoid (direct problem)
\frac{s}{b} = \int_{\sigma_1}^{\sigma_2}\sqrt{1 + k^2\sin^2\sigma}\,d\sigma,\quad \lambda_{12} = \omega_{12} - f\sin\alpha_0\int_{\sigma_1}^{\sigma_2}\frac{2-f}{1 + (1-f)\sqrt{1+k^2\sin^2\sigma}}\,d\sigma
with lat1 = 52.2, lon1 = 0.12, azimuth = 45, distance = 100 km, ellipsoid = WGS84
→ lat2 = 52.83083647, lon2 = 1.169145951, azimuth2 = 45.83253655

Capability: Geodesic destination point on the ellipsoid (direct problem)

Map projectionsm.ask(…)

UTM forward: lat = 51.5, lon = -0.13, ellipsoid = WGS84

easting = 699195.8759 m

northing = 5709335.019 m

k = 1.000087189

gamma = 2.246819165

zone = 30

hemisphere = N

Show the steps

Latitude and longitude to UTM coordinates (with scale factor and convergence)
E = 500\,000 + 0.9996\,A\eta,\quad N = N_0 + 0.9996\,A\xi
with lat = 51.5, lon = -0.13, ellipsoid = WGS84
→ easting = 699195.8759 m, northing = 5709335.019 m, k = 1.000087189 m, gamma = 2.246819165 m, zone = 30 m, hemisphere = N m

Capability: Latitude and longitude to UTM coordinates (with scale factor and convergence)

Ellipsoidal geodesym.ask(…)

Geodesic polygon area: lats = [48.0, 49.5, 47.2], lons = [2.0, 4.5, 6.1], ellipsoid = GRS80

area = 3.36697348e+10 m^2

perimeter = 850957.6213 m

Show the steps

Area and perimeter of a geodesic polygon on the ellipsoid
S = -\oint \frac{a^2}{2}\left[q(\varphi) - q(\varphi_c)\right]d\lambda,\quad q(\varphi) = (1-e^2)\left[\frac{\sin\varphi}{1-e^2\sin^2\varphi} + \frac{\operatorname{artanh}(e\sin\varphi)}{e}\right]
with lats = [48, 49.5, 47.2], lons = [2, 4.5, 6.1], ellipsoid = GRS80
→ area = 3.36697348e+10 m^2, perimeter = 850957.6213 m

Capability: Area and perimeter of a geodesic polygon on the ellipsoid

Remote sensingm.ask(…)

Vegetation indices: red = 0.08, nir = 0.45, blue = 0.04, L = 0.5

NDVI = 0.6981132075

EVI = 0.5674846626

SAVI = 0.5388349515

OSAVI = 0.5362318841

MSAVI = 0.5468871126

EVI2 = 0.5633373934

Show the steps

Vegetation indices: NDVI, EVI, SAVI, OSAVI, MSAVI and EVI2 from band reflectances
NDVI = \frac{\rho_{N} - \rho_R}{\rho_N + \rho_R},\ EVI = 2.5\frac{\rho_N - \rho_R}{\rho_N + 6\rho_R - 7.5\rho_B + 1},\ SAVI = (1+L)\frac{\rho_N - \rho_R}{\rho_N + \rho_R + L}
with red = 0.08, nir = 0.45, blue = 0.04, L = 0.5
→ NDVI = 0.6981132075, EVI = 0.5674846626, SAVI = 0.5388349515, OSAVI = 0.5362318841, MSAVI = 0.5468871126, EVI2 = 0.5633373934

Capability: Vegetation indices: NDVI, EVI, SAVI, OSAVI, MSAVI and EVI2 from band reflectances

Spatial statisticsm.solve(…)

m.solve("geo_morans_i", y=[8, 6, 6, 3, 2, 1], W=[[0, 1, 0, 0, 0, 0], [1, 0, 1, 0, 0, 0], [0, 1, 0, 1, 0, 0], [0, 0, 1, 0, 1, 0], [0, 0, 0, 1, 0, 1], [0, 0, 0, 0, 1, 0]], transform="row")

I = 0.65625

EI = -0.2

VI_norm = 0.140952381

z_norm = 2.280680117

p_norm = 0.02256738258

VI_rand = 0.1817926233

Show the steps

Global Moran's I spatial autocorrelation with normal and randomization inference
I = \frac{n}{S_0}\frac{\sum_i\sum_j w_{ij}z_iz_j}{\sum_i z_i^2},\quad E[I] = -\frac{1}{n-1},\quad V_N[I] = \frac{n^2S_1 - nS_2 + 3S_0^2}{(n^2-1)S_0^2} - E[I]^2
with y = [8, 6, 6, 3, 2, 1], W = [0 1 0 0 0 0; 1 0 1 0 0 0; 0 1 0 1 0 0; 0 0 1 0 1 0; 0 0 0 1 0 1; 0 0 0 0 1 0], transform = row
→ I = 0.65625, EI = -0.2, VI_norm = 0.140952381, z_norm = 2.280680117, p_norm = 0.02256738258, VI_rand = 0.1817926233, z_rand = 2.008225402, p_rand = 0.04461934483

Capability: Global Moran's I spatial autocorrelation with normal and randomization inference

Cartographym.ask(…)

Map scale from distances: map_distance = 5 cm, ground_distance = 2.5 km

scale = 50000

rf = 2e-05

Show the steps

Map scale (representative fraction) from a map and a ground distance
S = \frac{D}{d},\quad RF = \frac{1}{S}
with map_distance = 5 cm, ground_distance = 2.5 km
→ scale = 50000, rf = 2e-05

Capability: Map scale (representative fraction) from a map and a ground distance

Coordinate notationm.ask(…)

Decimal degrees to DMS: decimal = 51.477222, axis = latitude

degrees = 51

minutes = 28

seconds = 37.9992

hemisphere = N

Show the steps

Decimal degrees to degrees, minutes and seconds
d = \lfloor|x|\rfloor,\quad m = \lfloor60(|x| - d)\rfloor,\quad s = 3600\left(|x| - d - \tfrac{m}{60}\right)
with decimal = 51.477222, axis = latitude
→ degrees = 51, minutes = 28, seconds = 37.9992, hemisphere = N

Capability: Decimal degrees to degrees, minutes and seconds

Coordinate systems and conversionsm.ask(…)

Auxiliary latitudes: lat = 45, ellipsoid = WGS84

geocentric = 44.80757678

reduced = 44.90378785

conformal = 44.80768406

authalic = 44.87170287

rectifying = 44.85568199

Show the steps

Auxiliary latitudes: geocentric, reduced, conformal, authalic and rectifying
\tan\psi = (1-e^2)\tan\varphi,\ \tan\beta = (1-f)\tan\varphi,\ \chi = \operatorname{gd}\psi_{iso},\ \xi = \arcsin\frac{q}{q_p},\ \mu = \frac{\pi}{2}\frac{m(\varphi)}{m(\pi/2)}
with lat = 45, ellipsoid = WGS84
→ geocentric = 44.80757678, reduced = 44.90378785, conformal = 44.80768406, authalic = 44.87170287, rectifying = 44.85568199

Capability: Auxiliary latitudes: geocentric, reduced, conformal, authalic and rectifying

Use it

From Python or the command line

Every snippet was run as shown; the output under it is what it printed.

import geospatial_analytics_ai999x as m

a = m.ask("Geodesic destination: lat1 = 52.2, lon1 = 0.12, azimuth = 45, distance = 100 km, ellipsoid = WGS84")
print(a.outcome)
print(a.values)
for step in a.steps:
    print(step["name"])

Output

SOLVED
{'lat2': 52.8308364668, 'lon2': 1.16914595053, 'azimuth2': 45.832536554}
Geodesic destination point on the ellipsoid (direct problem)

import geospatial_analytics_ai999x as m

a = m.solve("geo_morans_i", y=[8, 6, 6, 3, 2, 1], W=[[0, 1, 0, 0, 0, 0], [1, 0, 1, 0, 0, 0], [0, 1, 0, 1, 0, 0], [0, 0, 1, 0, 1, 0], [0, 0, 0, 1, 0, 1], [0, 0, 0, 0, 1, 0]], transform="row")
print(a.outcome)
print(a.values)

Output

SOLVED
{'I': 0.65625, 'EI': -0.2, 'VI_norm': 0.140952380952, 'z_norm': 2.28068011671, 'p_norm': 0.0225673825813, 'VI_rand': 0.181792623299, 'z_rand': 2.00822540212, 'p_rand': 0.0446193448302}

import geospatial_analytics_ai999x as m

# an input is missing: the answer says what it needs instead of guessing
a = m.solve("geo_morans_i")
print(a.outcome)
print(a.message)

Output

NEEDS_INPUT
geo_morans_i needs y, W, transform (values at the n locations; spatial weights matrix (n x n, zero diagonal, row i lists the weights of i's neighbours); row-standardize the weights (row) or use them as given (none)). Example: solve('geo_morans_i', y=[8, 6, 6, 3, 2, 1], W=[[0, 1, 0, 0, 0, 0], [1, 0, 1, 0, 0, 0], [0, 1, 0, 1, 0, 0], [0, 0, 1, 0, 1, 0], [0, 0, 0, 1, 0, 1], [0, 0, 0, 0, 1, 0]], transform='row').

geospatial-analytics-ai999x "Geodesic destination: lat1 = 52.2, lon1 = 0.12, azimuth = 45, distance = 100 km, ellipsoid = WGS84"

Output

lat2 = 52.8308, lon2 = 1.16915, azimuth2 = 45.8325
1. Geodesic destination point on the ellipsoid (direct problem): Bessel's auxiliary sphere: solve b*I1(sigma) = s by Newton steps with the arc-length integral by Gauss-Legendre quadrature, then the longitude integral  with lat1 = 52.2, lon1 = 0.12, azimuth = 45, distance = 100 km, ellipsoid = WGS84  ->  lat2 = 52.8308, lon2 = 1.16915, azimuth2 = 45.8325
Assumptions: start latitude strictly between -90 and 90; 0 < distance <= 15000 km
…

Figures

Computed, not drawn by hand

Every value in these figures was computed by the package, each point a run of solve(); in a notebook, plot() and plot_relation() draw a capability's curves the same way.

Great circle against rhumb line, London to New York

\frac{s}{b} = \int_{\sigma_1}^{\sigma_2}\sqrt{1 + k^2\sin^2\sigma}\,d\sigma,\quad \lambda_{12} = \omega_{12} - f\sin\alpha_0\int_{\sigma_1}^{\sigma_2}\frac{2-f}{1 + (1-f)\sqrt{1+k^2\sin^2\sigma}}\,d\sigma

London to New York on a Mercator map (WGS84): the geodesic leaves at azimuth 288.373° and is 5585.81 km long; the rhumb line keeps bearing 258.088° and is 5810.37 km. Each path is 40 destination points computed by the package and projected by it; the rhumb line is straight on this map. Marked: the geodesic midpoint, latitude 52.385°, longitude -41.2786°.

World graticule in the Mollweide equal-area projection

2\theta + \sin2\theta = \pi\sin\varphi,\quad x = \frac{2\sqrt2}{\pi}R\,\Delta\lambda\cos\theta,\quad y = \sqrt2R\sin\theta

The world graticule in the Mollweide equal-area projection (sphere of radius 6371 km): every parallel and meridian is a chain of points projected by the package (576 points); marked: latitude 45°, longitude 90° → x = 7261.19 km, y = 5334.27 km.

The world from London: azimuth and distance

x = FE + s_{12}\sin\alpha_1,\quad y = FN + s_{12}\cos\alpha_1

Eight cities seen from London as the package computes them on the WGS84 ellipsoid — each arrow points along the initial azimuth, clockwise from north, with length the geodesic distance (an azimuthal equidistant view): New York 5585.81 km at 288.373°; Tokyo 9583.99 km at 31.6564°; Cape Town 9635.07 km at 164.613°; Rio de Janeiro 9251.66 km at 219.388°; Moscow 2509.19 km at 64.2336°; Mumbai 7199.57 km at 90.6929°; Reykjavík 1895.1 km at 326.304°; Nairobi 6805.98 km at 136.642°.

Mean centre and standard distance of a point pattern

SD = \sqrt{\frac{\sum_i (x_i - \bar x)^2 + (y_i - \bar y)^2}{n}}

A pattern of 11 points with the mean centre and standard distance computed by the package: marked: x = 52.4, y = 49.5982; the circle has radius 40.2158, the standard distance (root-mean-square distance of the points from the centre).

Line simplification at two tolerances

d_{max} = \max_k \operatorname{dist}(P_k, \overline{P_iP_j}) > \varepsilon \Rightarrow \text{keep } P_k

A 12-vertex line simplified by the Douglas–Peucker algorithm at two tolerances, the kept vertices as returned by the package: tolerance 0.4 keeps 8 vertices (0, 2, 4, 6, 7, 9, 10, 11); tolerance 1 keeps 4 vertices (0, 2, 4, 11).

Elevation grid and its flow accumulation

A(c) = 1 + \sum_{u\to c}A(u)

A 4 × 4 elevation grid (left, metres, cell size 10 m) and the D8 flow accumulation the package computes from it (right): each cell shows how many cells drain through it, itself included; the maximum, 16, is the outlet in the lowest corner. Each panel is shaded by its own maximum.

Formulas

The mathematics inside

12 of its 99 formulas; every one is in the capability browser below.

Geodesic destination point on the ellipsoid (direct problem)

Latitude and longitude to UTM coordinates (with scale factor and convergence)

E = 500\,000 + 0.9996\,A\eta,\quad N = N_0 + 0.9996\,A\xi

Area and perimeter of a geodesic polygon on the ellipsoid

S = -\oint \frac{a^2}{2}\left[q(\varphi) - q(\varphi_c)\right]d\lambda,\quad q(\varphi) = (1-e^2)\left[\frac{\sin\varphi}{1-e^2\sin^2\varphi} + \frac{\operatorname{artanh}(e\sin\varphi)}{e}\right]

Vegetation indices: NDVI, EVI, SAVI, OSAVI, MSAVI and EVI2 from band reflectances

NDVI = \frac{\rho_{N} - \rho_R}{\rho_N + \rho_R},\ EVI = 2.5\frac{\rho_N - \rho_R}{\rho_N + 6\rho_R - 7.5\rho_B + 1},\ SAVI = (1+L)\frac{\rho_N - \rho_R}{\rho_N + \rho_R + L}

Global Moran's I spatial autocorrelation with normal and randomization inference

I = \frac{n}{S_0}\frac{\sum_i\sum_j w_{ij}z_iz_j}{\sum_i z_i^2},\quad E[I] = -\frac{1}{n-1},\quad V_N[I] = \frac{n^2S_1 - nS_2 + 3S_0^2}{(n^2-1)S_0^2} - E[I]^2

Map scale (representative fraction) from a map and a ground distance

S = \frac{D}{d},\quad RF = \frac{1}{S}

Decimal degrees to degrees, minutes and seconds

d = \lfloor|x|\rfloor,\quad m = \lfloor60(|x| - d)\rfloor,\quad s = 3600\left(|x| - d - \tfrac{m}{60}\right)

Auxiliary latitudes: geocentric, reduced, conformal, authalic and rectifying

\tan\psi = (1-e^2)\tan\varphi,\ \tan\beta = (1-f)\tan\varphi,\ \chi = \operatorname{gd}\psi_{iso},\ \xi = \arcsin\frac{q}{q_p},\ \mu = \frac{\pi}{2}\frac{m(\varphi)}{m(\pi/2)}

Ground distance from a map distance and a representative fraction

D = d \times S

Decode a geohash: cell centre and size

[\varphi_{min}, \varphi_{max}] \times [\lambda_{min}, \lambda_{max}] \ \text{by alternate bisection}

Geocentric (ECEF) to geodetic coordinates

\varphi = \operatorname{atan2}\!\left(Z + e^2N\sin\varphi,\ \sqrt{X^2+Y^2}\right),\quad h = p\cos\varphi + Z\sin\varphi - \frac{a^2}{N}

Arc length along a parallel of latitude

s = N(\varphi)\cos\varphi\,|\Delta\lambda|,\quad N = \frac{a}{\sqrt{1-e^2\sin^2\varphi}}

Capabilities

Every capability, with an example

Search by name, formula, input or topic. Each capability opens to show what it takes, what it returns, when it applies, and a question it answered.

99 capabilities

Cartography 2

Ground distance from a map distance and a representative fractiongeo_map_scale_ground_distance

The real-world distance represented by a distance measured on a map of scale 1 : scale.

D = d \times S

Inputs

map_distance distance measured on the map [m]
scale scale denominator (1 : scale)

Returns

ground_distance ground distance [m]
ground_km ground distance in kilometres

Applies when

positive map distance and scale

Example

m.ask("Map scale ground distance: map_distance = 4.2 cm, scale = 50000")

→ ground_distance = 2100 m, ground_km = 2.1

Map scale (representative fraction) from a map and a ground distancegeo_map_scale_from_distances

The scale denominator and representative fraction of a map on which a known ground distance measures a given map distance.

S = \frac{D}{d},\quad RF = \frac{1}{S}

Inputs

map_distance distance on the map [m]
ground_distance corresponding ground distance [m]

Returns

scale scale denominator S (scale 1 : S)
rf representative fraction 1/S

Applies when

ground distance > map distance > 0

Example

m.ask("Map scale from distances: map_distance = 5 cm, ground_distance = 2.5 km")

→ scale = 50000, rf = 2e-05

Coordinate notation 5

Decimal degrees to degrees, minutes and secondsgeo_decimal_to_dms

A signed latitude or longitude in decimal degrees written as whole degrees, whole minutes, seconds and a hemisphere letter.

d = \lfloor|x|\rfloor,\quad m = \lfloor60(|x| - d)\rfloor,\quad s = 3600\left(|x| - d - \tfrac{m}{60}\right)

Inputs

decimal signed angle (decimal degrees)
axis whether the angle is a latitude (N/S) or a longitude (E/W) choice

Returns

degrees whole degrees
minutes whole minutes
seconds seconds
hemisphere N, S, E or W

Applies when

|latitude| <= 90, |longitude| <= 180

Example

m.ask("Decimal degrees to DMS: decimal = 51.477222, axis = latitude")

→ degrees = 51, minutes = 28, seconds = 37.9992, hemisphere = N

Decode a geohash: cell centre and sizegeo_geohash_decode

The centre latitude and longitude of a geohash cell and its half-height and half-width in degrees.

[\varphi_{min}, \varphi_{max}] \times [\lambda_{min}, \lambda_{max}] \ \text{by alternate bisection}

Inputs

geohash geohash string (characters 0-9, b-z without a, i, l, o) text

Returns

lat_center latitude of the cell centre
lon_center longitude of the cell centre
lat_error half the cell height (degrees)
lon_error half the cell width (degrees)

Applies when

1 to 12 valid characters

Example

m.ask("Geohash decode: geohash = u4pruydqqvj")

→ lat_center = 57.6491, lon_center = 10.4074, lat_error = 6.70552e-07, lon_error = 6.70552e-07

Degrees, minutes and seconds to decimal degreesgeo_dms_to_decimal

A latitude or longitude written in degrees, minutes and seconds with a hemisphere letter, converted to signed decimal degrees and radians.

\varphi = \pm\left(d + \frac{m}{60} + \frac{s}{3600}\right)

Inputs

degrees whole degrees (>= 0) integer
minutes whole minutes (0-59) integer
seconds seconds (0 <= s < 60)
hemisphere hemisphere letter (S and W are negative) choice

Returns

decimal signed decimal degrees
radians signed angle in radians

Applies when

0 <= minutes <= 59
0 <= seconds < 60
at most 90 (N/S) or 180 (E/W) degrees

Example

m.ask("DMS to decimal degrees: degrees = 51, minutes = 28, seconds = 38, hemisphere = N")

→ decimal = 51.4772, radians = 0.898447

Geohash of a point and its cellgeo_geohash_encode

The geohash of a latitude and longitude at a given number of characters, with the centre and half-sizes of the geohash cell.

\text{bits} = \lambda_1\varphi_1\lambda_2\varphi_2\ldots,\quad \text{character} = \text{base32}(5\ \text{bits})

Inputs

lat latitude (decimal degrees)
lon longitude (decimal degrees)
precision number of geohash characters (1-12) integer

Returns

geohash geohash
lat_center latitude of the cell centre
lon_center longitude of the cell centre
lat_error half the cell height (degrees)
lon_error half the cell width (degrees)

Applies when

1 to 12 characters

Example

m.ask("Geohash encode: lat = 57.64911, lon = 10.40744, precision = 11")

→ geohash = u4pruydqqvj, lat_center = 57.6491, lon_center = 10.4074, lat_error = 6.70552e-07, lon_error = 6.70552e-07

MGRS grid reference of a point (1 m)geo_mgrs_encode

The Military Grid Reference System (MGRS / USNG) reference of a WGS84 latitude and longitude to 1 m: zone, latitude band, 100 km square letters and 5-digit easting and northing.

\text{MGRS} = z\,B\,C_{100k}R_{100k}\,\lfloor E\rfloor \bmod 10^5\,\lfloor N\rfloor \bmod 10^5

Inputs

lat latitude (WGS84) (decimal degrees)
lon longitude (WGS84) (decimal degrees)

Returns

mgrs MGRS reference (no spaces)

Applies when

80 S <= latitude <= 84 N (polar regions use UPS)

Example

m.ask("MGRS grid reference: lat = 50, lon = 10")

→ mgrs = 32UNA7166639109

Coordinate systems and conversions 7

Auxiliary latitudes: geocentric, reduced, conformal, authalic and rectifyinggeo_auxiliary_latitudes

The five auxiliary latitudes of a geodetic latitude on the WGS84 or GRS80 ellipsoid (used by the sphere-based forms of map projections and geodesics).

\tan\psi = (1-e^2)\tan\varphi,\ \tan\beta = (1-f)\tan\varphi,\ \chi = \operatorname{gd}\psi_{iso},\ \xi = \arcsin\frac{q}{q_p},\ \mu = \frac{\pi}{2}\frac{m(\varphi)}{m(\pi/2)}

Inputs

lat geodetic latitude (decimal degrees)
ellipsoid reference ellipsoid choice

Returns

geocentric geocentric latitude (degrees)
reduced reduced (parametric) latitude (degrees)
conformal conformal latitude (degrees)
authalic authalic latitude (degrees)
rectifying rectifying latitude (degrees)

Applies when

latitude strictly between -90 and 90

Example

m.ask("Auxiliary latitudes: lat = 45, ellipsoid = WGS84")

→ geocentric = 44.8076, reduced = 44.9038, conformal = 44.8077, authalic = 44.8717, rectifying = 44.8557

Geocentric (ECEF) to geodetic coordinatesgeo_ecef_to_geodetic

Geodetic latitude, longitude and ellipsoidal height of a point given by its Earth-centred, Earth-fixed coordinates X, Y, Z on the WGS84 or GRS80 ellipsoid.

\varphi = \operatorname{atan2}\!\left(Z + e^2N\sin\varphi,\ \sqrt{X^2+Y^2}\right),\quad h = p\cos\varphi + Z\sin\varphi - \frac{a^2}{N}

Inputs

X geocentric X [m]
Y geocentric Y [m]
Z geocentric Z [m]
ellipsoid reference ellipsoid choice

Returns

lat geodetic latitude (degrees)
lon longitude (degrees)
h ellipsoidal height [m]

Applies when

6200 km <= distance from the centre <= 100000 km
not on the polar axis

Example

m.ask("ECEF to geodetic: X = 4517590.879 m, Y = 796576.168 m, Z = 4487348.409 m, ellipsoid = WGS84")

→ lat = 44.56, lon = 10, h = 49468.7 m

Geodetic to geocentric (ECEF) coordinatesgeo_geodetic_to_ecef

Earth-centred, Earth-fixed Cartesian coordinates X, Y, Z of a point given by geodetic latitude, longitude and ellipsoidal height on the WGS84 or GRS80 ellipsoid.

X = (N+h)\cos\varphi\cos\lambda,\ Y = (N+h)\cos\varphi\sin\lambda,\ Z = \left(N(1-e^2)+h\right)\sin\varphi

Inputs

lat geodetic latitude (decimal degrees)
lon longitude (decimal degrees)
h ellipsoidal height [m]
ellipsoid reference ellipsoid choice

Returns

X X (towards latitude 0, longitude 0) [m]
Y Y (towards longitude 90 E) [m]
Z Z (towards the north pole) [m]

Applies when

latitude between -90 and 90
height between -100 km and 100000 km

Example

m.ask("Geodetic to ECEF: lat = 45, lon = 10, h = 100 m, ellipsoid = WGS84")

→ X = 4.44903e+06 m, Y = 784484 m, Z = 4.48742e+06 m

Geodetic to local east-north-up (ENU) coordinatesgeo_geodetic_to_enu

Local tangent-plane coordinates east, north and up of a point relative to an origin (both given by latitude, longitude and ellipsoidal height), and the straight-line (slant) range.

\begin{bmatrix}e\\n\\u\end{bmatrix} = \begin{bmatrix}-\sin\lambda_0 & \cos\lambda_0 & 0\\ -\sin\varphi_0\cos\lambda_0 & -\sin\varphi_0\sin\lambda_0 & \cos\varphi_0\\ \cos\varphi_0\cos\lambda_0 & \cos\varphi_0\sin\lambda_0 & \sin\varphi_0\end{bmatrix}\begin{bmatrix}\Delta X\\ \Delta Y\\ \Delta Z\end{bmatrix}

Inputs

lat0 latitude of the local origin (decimal degrees)
lon0 longitude of the local origin (decimal degrees)
h0 ellipsoidal height of the origin [m]
lat latitude of the point (decimal degrees)
lon longitude of the point (decimal degrees)
h ellipsoidal height of the point [m]
ellipsoid reference ellipsoid choice

Returns

e east coordinate [m]
n north coordinate [m]
u up coordinate [m]
slant_range straight-line distance from the origin [m]

Applies when

origin not at a pole

Example

m.ask("Geodetic to ENU: lat0 = 46.017, lon0 = 7.75, h0 = 1673 m, lat = 45.976, lon = 7.658, h = 4531 m, ellipsoid = WGS84")

→ e = -7134.76 m, n = -4556.32 m, u = 2852.39 m, slant_range = 8933.14 m

Local east-north-up (ENU) to geodetic coordinatesgeo_enu_to_geodetic

Latitude, longitude and ellipsoidal height of a point given by its local east, north and up coordinates relative to an origin on the WGS84 or GRS80 ellipsoid.

\mathbf{X} = \mathbf{X}_0 + R^{T}(\varphi_0,\lambda_0)\,[e\ n\ u]^{T}

Inputs

lat0 latitude of the local origin (decimal degrees)
lon0 longitude of the local origin (decimal degrees)
h0 ellipsoidal height of the origin [m]
e east coordinate [m]
n north coordinate [m]
u up coordinate [m]
ellipsoid reference ellipsoid choice

Returns

lat latitude (degrees)
lon longitude (degrees)
h ellipsoidal height [m]

Applies when

origin not at a pole
the point not deep inside the Earth

Example

m.ask("ENU to geodetic: lat0 = 46.017, lon0 = 7.75, h0 = 1673 m, e = -7123.1 m, n = -4557.3 m, u = 2848.2 m, ellipsoid = WGS84")

→ lat = 45.976, lon = 7.65815, h = 4526.8 m

Look angles: azimuth, elevation and range from an observer to a targetgeo_look_angles

The azimuth (clockwise from true north), elevation angle above the local horizontal plane and slant range from an observer to a target, both given by latitude, longitude and ellipsoidal height (a satellite, aircraft or summit).

A = \operatorname{atan2}(e, n),\quad \varepsilon = \operatorname{atan2}\!\left(u, \sqrt{e^2+n^2}\right)

Inputs

lat0 observer latitude (decimal degrees)
lon0 observer longitude (decimal degrees)
h0 observer ellipsoidal height [m]
lat target latitude (decimal degrees)
lon target longitude (decimal degrees)
h target ellipsoidal height [m]
ellipsoid reference ellipsoid choice

Returns

azimuth azimuth (degrees clockwise from north)
elevation elevation angle (degrees, negative below the horizon)
slant_range straight-line distance [m]

Applies when

observer not at a pole
target distinct from the observer

Example

m.ask("Look angles: lat0 = 40, lon0 = -105, h0 = 1600 m, lat = 0, lon = -100, h = 35786000 m, ellipsoid = WGS84")

→ azimuth = 172.243, elevation = 43.4643, slant_range = 3.75137e+07 m

UTM zone, latitude band and central meridiangeo_utm_zone

The Universal Transverse Mercator zone number (with the Norway and Svalbard exceptions), the MGRS latitude band letter, the hemisphere and the zone's central meridian for a latitude and longitude.

z = \left\lfloor\frac{\lambda + 180^\circ}{6^\circ}\right\rfloor + 1,\quad \lambda_0 = 6z - 183^\circ

Inputs

lat latitude (decimal degrees)
lon longitude (decimal degrees)

Returns

zone UTM zone number (1-60)
band latitude band letter (C-X)
hemisphere N or S
central_meridian central meridian of the zone (degrees)

Applies when

80 S <= latitude <= 84 N

Example

m.ask("UTM zone: lat = 51.5, lon = -0.13")

→ zone = 30, band = U, hemisphere = N, central_meridian = -3

Ellipsoidal geodesy 9

Arc length along a parallel of latitudegeo_parallel_arc

The distance along a parallel of latitude between two longitudes (the shorter way round) on the WGS84 or GRS80 ellipsoid, and the radius of the parallel.

s = N(\varphi)\cos\varphi\,|\Delta\lambda|,\quad N = \frac{a}{\sqrt{1-e^2\sin^2\varphi}}

Inputs

lat latitude of the parallel (decimal degrees)
lon1 first longitude (decimal degrees)
lon2 second longitude (decimal degrees)
ellipsoid reference ellipsoid choice

Returns

distance arc length along the parallel [m]
radius radius of the parallel N cos(phi) [m]

Applies when

latitude strictly between -90 and 90
different longitudes

Example

m.ask("Parallel arc length: lat = 60, lon1 = 0, lon2 = 1, ellipsoid = WGS84")

→ distance = 55800 m, radius = 3.1971e+06 m

Area and perimeter of a geodesic polygon on the ellipsoidgeo_geodesic_polygon_area

The area enclosed by a polygon whose edges are geodesics on the WGS84 or GRS80 ellipsoid, and its perimeter (vertices given as latitude and longitude lists, in either orientation).

S = -\oint \frac{a^2}{2}\left[q(\varphi) - q(\varphi_c)\right]d\lambda,\quad q(\varphi) = (1-e^2)\left[\frac{\sin\varphi}{1-e^2\sin^2\varphi} + \frac{\operatorname{artanh}(e\sin\varphi)}{e}\right]

Inputs

lats latitudes of the vertices in order (decimal degrees) vector
lons longitudes of the vertices in the same order (decimal degrees) vector
ellipsoid reference ellipsoid choice

Returns

area enclosed area on the ellipsoid [m^2]
perimeter perimeter (sum of the geodesic edge lengths) [m]

Applies when

a simple polygon not enclosing a pole
vertices within 88 degrees of the equator
edges at most 90 degrees long

Example

m.ask("Geodesic polygon area: lats = [0, 0, 1, 1], lons = [0, 1, 1, 0], ellipsoid = WGS84")

→ area = 1.23088e+10 m^2, perimeter = 443771 m

Area of a latitude-longitude cell on the ellipsoidgeo_graticule_cell_area

The area of the region between two parallels and two meridians (a graticule cell or grid square) on the WGS84 or GRS80 ellipsoid, and its share of the whole ellipsoid's surface.

A = \frac{a^2}{2}\left[q(\varphi_2) - q(\varphi_1)\right](\lambda_2 - \lambda_1)

Inputs

lat1 southern latitude (decimal degrees)
lat2 northern latitude (decimal degrees)
lon1 western longitude (decimal degrees)
lon2 eastern longitude (east of lon1) (decimal degrees)
ellipsoid reference ellipsoid choice

Returns

area area of the cell [m^2]
fraction_of_ellipsoid area as a fraction of the ellipsoid's total surface

Applies when

lat1 < lat2
0 < lon2 - lon1 <= 360

Example

m.ask("Graticule cell area: lat1 = 0, lat2 = 1, lon1 = 0, lon2 = 1, ellipsoid = WGS84")

→ area = 1.23085e+10 m^2, fraction_of_ellipsoid = 2.41311e-05

Destination along a rhumb line (loxodrome) on the ellipsoidgeo_rhumb_line_destination

The point reached by travelling a given distance on a constant bearing (rhumb line) from a start point on the WGS84 or GRS80 ellipsoid.

m_2 = m_1 + s\cos\alpha,\quad \Delta\lambda = \tan\alpha\,(\psi_2 - \psi_1)

Inputs

lat1 latitude of the start point (decimal degrees)
lon1 longitude of the start point (decimal degrees)
bearing constant bearing, clockwise from north (decimal degrees)
distance distance travelled [m]
ellipsoid reference ellipsoid choice

Returns

lat2 latitude of the destination (degrees)
lon2 longitude of the destination (degrees)

Applies when

the line stays more than 1 degree from the poles
distance > 0

Example

m.ask("Rhumb line destination: lat1 = 50, lon1 = -5, bearing = 30, distance = 500 km, ellipsoid = WGS84")

→ lat2 = 53.8917, lon2 = -1.36121

Geodesic destination point on the ellipsoid (direct problem)geo_geodesic_destination

The point reached by travelling a given distance along a geodesic from a start point at a given initial azimuth on the WGS84 or GRS80 ellipsoid, and the azimuth of the geodesic on arrival.

Inputs

lat1 latitude of the start point (decimal degrees)
lon1 longitude of the start point (decimal degrees)
azimuth initial azimuth, clockwise from north (decimal degrees)
distance geodesic distance travelled [m]
ellipsoid reference ellipsoid choice

Returns

lat2 latitude of the destination (degrees)
lon2 longitude of the destination (degrees)
azimuth2 forward azimuth at the destination (degrees clockwise from north)

Applies when

start latitude strictly between -90 and 90
0 < distance <= 15000 km

Example

m.ask("Geodesic destination: lat1 = 52.2, lon1 = 0.12, azimuth = 45, distance = 100 km, ellipsoid = WGS84")

→ lat2 = 52.8308, lon2 = 1.16915, azimuth2 = 45.8325

Meridian arc length between two latitudesgeo_meridian_arc

The distance along a meridian of the WGS84 or GRS80 ellipsoid between two latitudes, and the meridian distance of each latitude from the equator.

m(\varphi) = \frac{a}{1+n}\left[\left(1+\tfrac{n^2}{4}+\tfrac{n^4}{64}\right)\varphi - \tfrac32\left(n-\tfrac{n^3}{8}\right)\sin2\varphi + \tfrac{15}{16}\left(n^2-\tfrac{n^4}{4}\right)\sin4\varphi - \tfrac{35}{48}n^3\sin6\varphi + \tfrac{315}{512}n^4\sin8\varphi\right]

Inputs

lat1 first latitude (decimal degrees)
lat2 second latitude (decimal degrees)
ellipsoid reference ellipsoid choice

Returns

distance meridian arc length between the latitudes [m]
m1 meridian distance from the equator to lat1 (negative south) [m]
m2 meridian distance from the equator to lat2 (negative south) [m]

Applies when

latitudes between -90 and 90 and different

Example

m.ask("Meridian arc length: lat1 = 0, lat2 = 90, ellipsoid = WGS84")

→ distance = 1.0002e+07 m, m1 = 0 m, m2 = 1.0002e+07 m

Midpoint of a geodesic on the ellipsoidgeo_geodesic_midpoint

The point halfway along the shortest path (geodesic) between two points on the WGS84 or GRS80 ellipsoid, with the geodesic length.

P_{mid} = \mathrm{direct}\left(P_1,\ \alpha_1,\ s_{12}/2\right)

Inputs

lat1 latitude of the first point (decimal degrees)
lon1 longitude of the first point (decimal degrees)
lat2 latitude of the second point (decimal degrees)
lon2 longitude of the second point (decimal degrees)
ellipsoid reference ellipsoid choice

Returns

lat_mid latitude of the midpoint (degrees)
lon_mid longitude of the midpoint (degrees)
distance length of the geodesic between the two points [m]

Applies when

distinct points at most 150 degrees apart
latitudes strictly between -90 and 90

Example

m.ask("Geodesic midpoint: lat1 = 51.5, lon1 = -0.13, lat2 = 40.71, lon2 = -74.01, ellipsoid = WGS84")

→ lat_mid = 52.385, lon_mid = -41.2786, distance = 5.58581e+06 m

Radii of curvature of the ellipsoid at a latitudegeo_radii_of_curvature

The meridional radius of curvature M, the prime-vertical radius of curvature N and the Gaussian mean radius sqrt(MN) of the WGS84 or GRS80 ellipsoid at a latitude.

M = \frac{a(1-e^2)}{(1-e^2\sin^2\varphi)^{3/2}},\quad N = \frac{a}{\sqrt{1-e^2\sin^2\varphi}},\quad R_G = \sqrt{MN}

Inputs

lat latitude (decimal degrees)
ellipsoid reference ellipsoid choice

Returns

M meridional radius of curvature [m]
N prime-vertical radius of curvature [m]
R_gauss Gaussian mean radius sqrt(M N) [m]

Applies when

latitude between -90 and 90

Example

m.ask("Radii of curvature: lat = 45, ellipsoid = WGS84")

→ M = 6.36738e+06 m, N = 6.38884e+06 m, R_gauss = 6.3781e+06 m

Rhumb-line (loxodrome) distance and bearing on the ellipsoidgeo_rhumb_line_inverse

The constant bearing and the length of the rhumb line (loxodrome) between two points on the WGS84 or GRS80 ellipsoid, the shorter way round in longitude.

\tan\alpha = \frac{\Delta\lambda}{\Delta\psi},\quad s = \frac{\Delta m}{\cos\alpha},\quad \psi = \operatorname{arsinh}(\tan\varphi) - e\operatorname{artanh}(e\sin\varphi)

Inputs

lat1 latitude of the start point (decimal degrees)
lon1 longitude of the start point (decimal degrees)
lat2 latitude of the end point (decimal degrees)
lon2 longitude of the end point (decimal degrees)
ellipsoid reference ellipsoid choice

Returns

distance rhumb-line length [m]
bearing constant bearing (degrees clockwise from north)

Applies when

latitudes within 89 degrees of the equator
distinct points

Example

m.ask("Rhumb line distance and bearing: lat1 = 50, lon1 = -5, lat2 = 58, lon2 = 3, ellipsoid = WGS84")

→ distance = 1.03244e+06 m, bearing = 30.4064

Facility location 3

Closest facility for every node of a networkgeo_closest_facility

For every node of a road network, the nearest facility by shortest-path distance and that distance, with the mean distance to the nearest facility.

f^*(i) = \arg\min_{f\in F} d(i, f)

Inputs

edges undirected network as weighted edges between nodes numbered 0 .. n-1 edges
facilities facility nodes vector

Returns

nearest nearest facility of each node
distance distance to it
mean_distance mean distance to the nearest facility

Applies when

every node reaches a facility
no node equally near two facilities

Example

m.ask("Closest facility: edges = [0-1 2; 1-2 3; 0-2 6; 2-3 1; 3-4 4], facilities = [0, 3]")

→ nearest = [0.0, 0.0, 3.0, 3.0, 3.0], distance = [0.0, 2.0, 1.0, 0.0, 4.0], mean_distance = 1.4

Maximal covering location problem (MCLP) on a networkgeo_maximal_covering_location

The largest demand that p facilities, chosen from candidate nodes, can cover within a coverage distance on a road network, and the covered share.

\max_{S\subseteq C,|S| = p}\ \sum_i w_i\,\mathbf 1\left[\min_{j\in S}d_{ij}\le s\right]

Inputs

edges undirected network as weighted edges between nodes numbered 0 .. n-1 edges
demand demand at nodes 0 .. n-1 vector
candidates candidate facility nodes vector
p number of facilities integer
coverage coverage distance (network distance <= coverage counts)

Returns

covered maximum covered demand
fraction covered share of the total demand

Applies when

at most 20000 candidate combinations

Example

m.ask("Maximal covering location: edges = [0-1 2; 1-2 3; 0-2 6; 2-3 1; 3-4 4; 4-5 2; 1-5 7], demand = [10, 20, 15, 30, 5, 25], candidates = [0, 1, 2, 3, 4, 5], p = 2, coverage = 3")

→ covered = 95, fraction = 0.904762

p-median facility location on a networkgeo_p_median

The p facility sites, chosen from candidate nodes, that minimize the total demand-weighted shortest-path distance from every node to its nearest open facility.

\min_{S\subseteq C, |S| = p}\ \sum_i w_i\min_{j\in S} d_{ij}

Inputs

edges undirected network as weighted edges between nodes numbered 0 .. n-1 edges
demand demand at nodes 0 .. n-1 vector
candidates candidate facility nodes vector
p number of facilities to open integer

Returns

objective minimum total weighted distance
facilities chosen facility nodes (ascending)
mean_distance demand-weighted mean distance

Applies when

connected network
at most 20000 candidate combinations
a unique optimal set

Example

m.ask("p-median: edges = [0-1 2; 1-2 3; 0-2 6; 2-3 1; 3-4 4; 4-5 2; 1-5 7], demand = [10, 20, 15, 30, 5, 25], candidates = [0, 1, 2, 3, 4, 5], p = 2")

→ objective = 150, facilities = [2.0, 5.0], mean_distance = 1.42857

Geostatistics 5

Empirical (Matheron) semivariogram in distance binsgeo_empirical_semivariogram

The classical (Matheron) empirical semivariogram of sampled values: half the mean squared difference of all pairs whose separation falls in each distance bin, and the pair counts.

\hat\gamma(h_k) = \frac{1}{2N_k}\sum_{(i,j)\in k}(z_i - z_j)^2

Inputs

x x coordinates of the samples vector
y y coordinates of the samples vector
z sample values vector
bin_edges increasing distance bin edges (a bin includes its lower edge) vector

Returns

gamma semivariance of each bin
counts number of pairs in each bin

Applies when

every bin contains at least one pair
bins [e_k, e_k+1)

Example

m.ask("Empirical semivariogram: x = [0, 1, 2, 3, 4], y = [0, 0, 0, 0, 0], z = [1, 3, 2, 5, 4], bin_edges = [0, 1.5, 2.5, 4.5]")

→ gamma = [1.875, 1.5, 4.33333333333], counts = [4.0, 3.0, 3.0]

Fit an exponential or Gaussian variogram model (nugget, partial sill, range)geo_variogram_model_fit

The nugget, partial sill, range and sill of an exponential or Gaussian variogram model fitted by ordinary least squares to an empirical semivariogram (lags and semivariances).

\gamma(h) = c_0 + c\left(1 - e^{-3h/a}\right)\ \text{or}\ c_0 + c\left(1 - e^{-(7h/4a)^2}\right),\quad \min_{c_0, c, a}\sum_k\left(\hat\gamma_k - \gamma(h_k)\right)^2

Inputs

h lag distances vector
gamma semivariances at the lags vector
model variogram model (pykrige parameterization, practical range) choice

Returns

nugget nugget c0
psill partial sill c
range practical range a
sill sill c0 + c

Applies when

the least-squares range has an interior minimum (between min(h)/20 and 20 max(h))

Example

m.ask("Variogram model fit: h = [1, 2, 4, 6, 8, 12, 16, 20], gamma = [1.6, 2.9, 4.6, 5.6, 6.3, 6.8, 7, 7.05], model = exponential")

→ nugget = -0.0628821, psill = 7.16435, range = 11.3285, sill = 7.10147

Inverse distance weighted (IDW) interpolation at a pointgeo_idw_interpolation

The inverse-distance-weighted (Shepard) estimate of a variable at a location from sample points, with a stated power (2 is the usual choice).

\hat z(\mathbf p) = \frac{\sum_i d_i^{-p}z_i}{\sum_i d_i^{-p}}

Inputs

x x coordinates of the samples vector
y y coordinates of the samples vector
z sample values vector
px x of the location
py y of the location
power distance power p (> 0)

Returns

estimate interpolated value

Applies when

power > 0 (at most 10)
a location on a sample returns that sample's value

Example

m.ask("IDW interpolation: x = [0, 10, 0, 10], y = [0, 0, 10, 10], z = [1, 2, 3, 4], px = 2, py = 3, power = 2")

→ estimate = 1.65898

Ordinary kriging estimate and kriging variance at a pointgeo_ordinary_kriging

The ordinary kriging estimate of a variable at a location from samples, with a given variogram model (nugget, sill, range), and the kriging variance.

\begin{bmatrix}\Gamma & \mathbf 1\\ \mathbf 1^T & 0\end{bmatrix}\begin{bmatrix}\boldsymbol\lambda\\ \mu\end{bmatrix} = \begin{bmatrix}\boldsymbol\gamma_0\\ 1\end{bmatrix},\quad \hat z = \boldsymbol\lambda^T\mathbf z,\quad \sigma^2_{OK} = \boldsymbol\lambda^T\boldsymbol\gamma_0 + \mu

Inputs

x x coordinates of the samples vector
y y coordinates of the samples vector
z sample values vector
px x of the location
py y of the location
model variogram model (pykrige parameterization) choice
nugget nugget
sill sill (nugget + partial sill)
range range (practical range for exponential and gaussian)

Returns

estimate kriging estimate
variance kriging variance

Applies when

distinct sample locations
sill > nugget >= 0, range > 0
the variogram is 0 at zero distance (exact interpolation)

Example

m.ask("Ordinary kriging: x = [0, 10, 0, 10], y = [0, 0, 10, 10], z = [1, 2, 3, 4], px = 4, py = 6, model = spherical, nugget = 0, sill = 2, range = 20")

→ estimate = 2.6048, variance = 0.877115

Polynomial trend surface (first or second order) by least squaresgeo_trend_surface

The coefficients and R^2 of a first-order (plane) or second-order (quadratic) trend surface fitted to sampled values by least squares.

z = b_0 + b_1x + b_2y\ \left[+\, b_3x^2 + b_4xy + b_5y^2\right]

Inputs

x x coordinates vector
y y coordinates vector
z values vector
order order of the surface (1 or 2) integer

Returns

coefficients b0, b1 (x), b2 (y)[, b3 (x^2), b4 (xy), b5 (y^2)]
r_squared coefficient of determination

Applies when

more samples than coefficients
values not all equal

Example

m.ask("Trend surface: x = [0, 1, 2, 0, 1, 2], y = [0, 0, 0, 1, 1, 1], z = [1, 3, 5.2, 2, 4.1, 6], order = 1")

→ coefficients = [1.01666666667, 2.05, 0.966666666667], r_squared = 0.99872

Map projections 22

Albers equal-area conic projection: coordinates and scale factorsgeo_albers_equal_area_forward

Easting and northing of a point in the ellipsoidal Albers equal-area conic projection with two standard parallels, with the scale factors along the parallel (k) and the meridian (h = 1/k).

n = \frac{m_1^2 - m_2^2}{q_2 - q_1},\quad \rho = \frac{a\sqrt{C - nq}}{n},\quad x = FE + \rho\sin\theta,\ y = FN + \rho_0 - \rho\cos\theta,\ k = \frac{\rho n}{am} = \frac1h

Inputs

lat latitude (decimal degrees)
lon longitude (decimal degrees)
lat1 first standard parallel (decimal degrees)
lat2 second standard parallel (decimal degrees)
lat0 latitude of origin (decimal degrees)
lon0 central meridian (decimal degrees)
ellipsoid reference ellipsoid choice
FE false easting (default 0) [m] optional
FN false northing (default 0) [m] optional

Returns

x easting [m]
y northing [m]
k scale factor along the parallel
h scale factor along the meridian

Applies when

parallels within 89 degrees of the equator and not symmetric about it

Example

m.ask("Albers equal-area forward: lat = 40, lon = -100, lat1 = 29.5, lat2 = 45.5, lat0 = 23, lon0 = -96, ellipsoid = GRS80")

→ x = -338391 m, y = 1.8941e+06 m, k = 0.990969, h = 1.00911

Azimuthal equidistant projection (ellipsoidal, geodesic)geo_azimuthal_equidistant_forward

Easting and northing of a point in the ellipsoidal azimuthal equidistant projection about a centre: the geodesic distance from the centre plotted at the geodesic azimuth (used for range rings and airline maps).

x = FE + s_{12}\sin\alpha_1,\quad y = FN + s_{12}\cos\alpha_1

Inputs

lat latitude (decimal degrees)
lon longitude (decimal degrees)
lat0 latitude of the centre (decimal degrees)
lon0 longitude of the centre (decimal degrees)
ellipsoid reference ellipsoid choice
FE false easting (default 0) [m] optional
FN false northing (default 0) [m] optional

Returns

x easting [m]
y northing [m]
distance geodesic distance from the centre [m]
azimuth geodesic azimuth from the centre (degrees)

Applies when

point distinct from the centre and at most 150 degrees from it
centre within 89 degrees of the equator

Example

m.ask("Azimuthal equidistant forward: lat = 40.71, lon = -74.01, lat0 = 51.5, lon0 = -0.13, ellipsoid = WGS84")

→ x = -5.30106e+06 m, y = 1.76067e+06 m, distance = 5.58581e+06 m, azimuth = 288.373

Equirectangular (equidistant cylindrical / plate carree) projectiongeo_equirectangular_forward

x and y of a point in the equirectangular (equidistant cylindrical) projection of a sphere of radius R with a standard parallel (lat_ts = 0 is the plate carree).

x = R(\lambda-\lambda_0)\cos\varphi_{ts},\quad y = R(\varphi-\varphi_0)

Inputs

lat latitude (decimal degrees)
lon longitude (decimal degrees)
lat_ts standard parallel (true scale) (decimal degrees)
lat0 latitude of origin (decimal degrees)
lon0 central meridian (decimal degrees)
R radius of the sphere [m]

Returns

x x [m]
y y [m]

Applies when

sphere of radius R > 0

Example

m.ask("Equirectangular forward: lat = 45, lon = 10, lat_ts = 0, lat0 = 0, lon0 = 0, R = 6371000 m")

→ x = 1.11195e+06 m, y = 5.00377e+06 m

Gnomonic projection (sphere)geo_gnomonic_forward

x and y of a point in the gnomonic projection (great circles become straight lines) of a sphere of radius R centred on a point.

x = \frac{R\cos\varphi\sin\Delta\lambda}{\cos c},\quad y = \frac{R\left(\cos\varphi_0\sin\varphi - \sin\varphi_0\cos\varphi\cos\Delta\lambda\right)}{\cos c}

Inputs

lat latitude (decimal degrees)
lon longitude (decimal degrees)
lat0 latitude of the centre (decimal degrees)
lon0 longitude of the centre (decimal degrees)
R radius of the sphere [m]

Returns

x x [m]
y y [m]

Applies when

the point within about 84 degrees of the centre (cos c >= 0.1)

Example

m.ask("Gnomonic forward: lat = 40, lon = 10, lat0 = 45, lon0 = 0, R = 6371000 m")

→ x = 857808 m, y = -508966 m

Inverse Albers equal-area conic projectiongeo_albers_equal_area_inverse

Latitude and longitude of a point given by its easting and northing in the ellipsoidal Albers equal-area conic projection with two standard parallels.

q = \frac{C - \rho^2n^2/a^2}{n},\quad \varphi:\ q(\varphi) = q,\quad \lambda = \lambda_0 + \frac{\theta}{n}

Inputs

x easting [m]
y northing [m]
lat1 first standard parallel (decimal degrees)
lat2 second standard parallel (decimal degrees)
lat0 latitude of origin (decimal degrees)
lon0 central meridian (decimal degrees)
ellipsoid reference ellipsoid choice
FE false easting (default 0) [m] optional
FN false northing (default 0) [m] optional

Returns

lat latitude (degrees)
lon longitude (degrees)

Applies when

parallels within 89 degrees of the equator and not symmetric about it
the point inside the projected area

Example

m.ask("Albers equal-area inverse: x = -338796.0 m, y = 1952260.0 m, lat1 = 29.5, lat2 = 45.5, lat0 = 23, lon0 = -96, ellipsoid = GRS80")

→ lat = 40.5185, lon = -100.034

Inverse Lambert conformal conic projection (two standard parallels)geo_lambert_conformal_conic_inverse

Latitude and longitude of a point given by its easting and northing in the ellipsoidal Lambert conformal conic projection with two standard parallels.

\rho = \operatorname{sgn}(n)\sqrt{x'^2 + (\rho_0 - y')^2},\quad t = \left(\frac{\rho}{aF}\right)^{1/n},\quad \lambda = \lambda_0 + \frac{\theta}{n}

Inputs

x easting [m]
y northing [m]
lat1 first standard parallel (decimal degrees)
lat2 second standard parallel (decimal degrees)
lat0 latitude of origin (decimal degrees)
lon0 central meridian (decimal degrees)
ellipsoid reference ellipsoid choice
FE false easting (default 0) [m] optional
FN false northing (default 0) [m] optional

Returns

lat latitude (degrees)
lon longitude (degrees)

Applies when

parallels within 89 degrees of the equator and not symmetric about it

Example

m.ask("Lambert conformal conic inverse: x = -339545.0 m, y = 1937380.0 m, lat1 = 33, lat2 = 45, lat0 = 23, lon0 = -96, ellipsoid = GRS80")

→ lat = 40.2966, lon = -100.016

Inverse Mercator projection (ellipsoidal)geo_mercator_inverse

Latitude and longitude of a point given by its x and y in the ellipsoidal Mercator projection with a central meridian and standard parallel.

\psi = \frac{y}{k_0a},\quad \tau' = \sinh\psi,\quad \tau = \tan\varphi:\ \tau\sqrt{1+\sigma^2} - \sigma\sqrt{1+\tau^2} = \tau'

Inputs

x easting [m]
y northing [m]
lon0 central meridian (decimal degrees)
lat_ts latitude of true scale (decimal degrees)
ellipsoid reference ellipsoid choice

Returns

lat latitude (degrees)
lon longitude (degrees)

Applies when

|lat_ts| <= 85

Example

m.ask("Mercator inverse: x = -14471.5 m, y = 6680752.0 m, lon0 = 0, lat_ts = 0, ellipsoid = WGS84")

→ lat = 51.5224, lon = -0.13

Lambert azimuthal equal-area projection (ellipsoidal)geo_lambert_azimuthal_equal_area_forward

Easting and northing of a point in the ellipsoidal oblique Lambert azimuthal equal-area projection about a centre (as in ETRS89-LAEA, EPSG:3035).

B = R_q\sqrt{\frac{2}{1 + \sin\beta_1\sin\beta + \cos\beta_1\cos\beta\cos\Delta\lambda}},\ x = BD\cos\beta\sin\Delta\lambda,\ y = \frac{B}{D}\left(\cos\beta_1\sin\beta - \sin\beta_1\cos\beta\cos\Delta\lambda\right)

Inputs

lat latitude (decimal degrees)
lon longitude (decimal degrees)
lat0 latitude of the centre (decimal degrees)
lon0 longitude of the centre (decimal degrees)
ellipsoid reference ellipsoid choice
FE false easting (default 0) [m] optional
FN false northing (default 0) [m] optional

Returns

x easting [m]
y northing [m]

Applies when

centre within 89 degrees of the equator
point within 120 degrees of the centre

Example

m.ask("Lambert azimuthal equal-area forward: lat = 50, lon = 5, lat0 = 52, lon0 = 10, ellipsoid = GRS80, FE = 4321000 m, FN = 3210000 m")

→ x = 3.9628e+06 m, y = 2.99972e+06 m

Lambert conformal conic projection (two standard parallels): coordinates, scale factor and convergencegeo_lambert_conformal_conic_forward

Easting and northing of a point in the ellipsoidal Lambert conformal conic projection with two standard parallels (EPSG method 9802), with the point scale factor and the grid convergence.

n = \frac{\ln m_1 - \ln m_2}{\ln t_1 - \ln t_2},\ \rho = aFt^n,\ x = FE + \rho\sin\theta,\ y = FN + \rho_0 - \rho\cos\theta,\ k = \frac{\rho n}{am}

Inputs

lat latitude (decimal degrees)
lon longitude (decimal degrees)
lat1 first standard parallel (decimal degrees)
lat2 second standard parallel (decimal degrees)
lat0 latitude of origin (decimal degrees)
lon0 central meridian (decimal degrees)
ellipsoid reference ellipsoid choice
FE false easting (default 0) [m] optional
FN false northing (default 0) [m] optional

Returns

x easting [m]
y northing [m]
k point scale factor
gamma grid convergence (degrees)

Applies when

parallels within 89 degrees of the equator and not symmetric about it

Example

m.ask("Lambert conformal conic forward: lat = 40, lon = -100, lat1 = 33, lat2 = 45, lat0 = 23, lon0 = -96, ellipsoid = GRS80")

→ x = -339644 m, y = 1.9046e+06 m, k = 0.994666, gamma = -2.52199

Lambert cylindrical equal-area projection (sphere)geo_cylindrical_equal_area_forward

x and y of a point in the cylindrical equal-area projection of a sphere of radius R with a standard parallel (0 Lambert, 30 Behrmann, 45 Gall-Peters).

x = R(\lambda-\lambda_0)\cos\varphi_s,\quad y = \frac{R\sin\varphi}{\cos\varphi_s}

Inputs

lat latitude (decimal degrees)
lon longitude (decimal degrees)
lat_ts standard parallel (decimal degrees)
lon0 central meridian (decimal degrees)
R radius of the sphere [m]

Returns

x x [m]
y y [m]

Applies when

standard parallel strictly between -90 and 90

Example

m.ask("Cylindrical equal-area forward: lat = 45, lon = 10, lat_ts = 45, lon0 = 0, R = 6371000 m")

→ x = 786267 m, y = 6.371e+06 m

Latitude and longitude to UTM coordinates (with scale factor and convergence)geo_utm_forward

UTM zone, hemisphere, easting and northing of a point on the WGS84 or GRS80 ellipsoid, with the point scale factor and grid convergence.

E = 500\,000 + 0.9996\,A\eta,\quad N = N_0 + 0.9996\,A\xi

Inputs

lat latitude (decimal degrees)
lon longitude (decimal degrees)
ellipsoid reference ellipsoid choice

Returns

easting UTM easting [m]
northing UTM northing [m]
zone UTM zone number
hemisphere N or S
k point scale factor
gamma grid convergence (degrees)

Applies when

80 S <= latitude <= 84 N

Example

m.ask("UTM forward: lat = 51.5, lon = -0.13, ellipsoid = WGS84")

→ easting = 699196 m, northing = 5.70934e+06 m, k = 1.00009, gamma = 2.24682, zone = 30, hemisphere = N

Latitude and longitude to Web Mercator (EPSG:3857)geo_web_mercator_forward

Web Mercator (pseudo-Mercator, EPSG:3857) x and y of a WGS84 latitude and longitude, as used by web maps.

x = R\lambda,\quad y = R\ln\tan\left(\frac{\pi}{4} + \frac{\varphi}{2}\right),\quad R = 6378137\ \mathrm{m}

Inputs

lat latitude (WGS84) (decimal degrees)
lon longitude (WGS84) (decimal degrees)

Returns

x Web Mercator x [m]
y Web Mercator y [m]

Applies when

|latitude| <= 85.06 degrees (the square web-map extent)

Example

m.ask("Web Mercator forward: lat = 51.5, lon = -0.13")

→ x = -14471.5 m, y = 6.71022e+06 m

Mercator projection (ellipsoidal): coordinates and scale factorgeo_mercator_forward

x and y of a point in the ellipsoidal Mercator projection with a central meridian and standard parallel (lat_ts = 0 is World Mercator, EPSG:3395), with the point scale factor.

x = k_0a(\lambda-\lambda_0),\quad y = k_0a\left[\operatorname{arsinh}(\tan\varphi) - e\operatorname{artanh}(e\sin\varphi)\right],\quad k = \frac{k_0\sqrt{1-e^2\sin^2\varphi}}{\cos\varphi}

Inputs

lat latitude (decimal degrees)
lon longitude (decimal degrees)
lon0 central meridian (decimal degrees)
lat_ts latitude of true scale (standard parallel) (decimal degrees)
ellipsoid reference ellipsoid choice

Returns

x easting [m]
y northing [m]
k point scale factor

Applies when

|latitude| <= 89.5
|lat_ts| <= 85

Example

m.ask("Mercator forward: lat = 51.5, lon = -0.13, lon0 = 0, lat_ts = 0, ellipsoid = WGS84")

→ x = -14471.5 m, y = 6.67676e+06 m, k = 1.60309

Mollweide equal-area projection (sphere)geo_mollweide_forward

x and y of a point in the Mollweide (homalographic) equal-area projection of a sphere of radius R.

2\theta + \sin2\theta = \pi\sin\varphi,\quad x = \frac{2\sqrt2}{\pi}R\,\Delta\lambda\cos\theta,\quad y = \sqrt2R\sin\theta

Inputs

lat latitude (decimal degrees)
lon longitude (decimal degrees)
lon0 central meridian (decimal degrees)
R radius of the sphere [m]

Returns

x x [m]
y y [m]

Applies when

sphere of radius R > 0

Example

m.ask("Mollweide forward: lat = 45, lon = 90, lon0 = 0, R = 6371000 m")

→ x = 7.26119e+06 m, y = 5.33427e+06 m

Orthographic projection (sphere)geo_orthographic_forward

x and y of a point in the orthographic projection (the view of a globe from infinitely far away) of a sphere of radius R centred on a point.

x = R\cos\varphi\sin\Delta\lambda,\quad y = R\left(\cos\varphi_0\sin\varphi - \sin\varphi_0\cos\varphi\cos\Delta\lambda\right)

Inputs

lat latitude (decimal degrees)
lon longitude (decimal degrees)
lat0 latitude of the centre (decimal degrees)
lon0 longitude of the centre (decimal degrees)
R radius of the sphere [m]

Returns

x x [m]
y y [m]

Applies when

the point on the visible hemisphere (cos c >= 0)

Example

m.ask("Orthographic forward: lat = 40, lon = 10, lat0 = 45, lon0 = 0, R = 6371000 m")

→ x = 847485 m, y = -502841 m

Polar stereographic projection (standard parallel): coordinates and scale factorgeo_polar_stereographic_forward

Easting and northing of a point in the ellipsoidal polar stereographic projection defined by a latitude of true scale (EPSG variant B, as in EPSG:3413 and EPSG:3031), with the point scale factor.

\rho = a\,m_c\frac{t}{t_c},\quad t = e^{-\psi},\quad x = FE + \rho\sin\Delta\lambda,\ y = FN \mp \rho\cos\Delta\lambda,\quad k = \frac{\rho}{am}

Inputs

lat latitude (decimal degrees)
lon longitude (decimal degrees)
pole pole at the projection centre choice
lat_ts latitude of true scale (decimal degrees)
lon0 longitude of the downward (north pole) or upward (south pole) y axis (decimal degrees)
ellipsoid reference ellipsoid choice
FE false easting (default 0) [m] optional
FN false northing (default 0) [m] optional

Returns

x easting [m]
y northing [m]
k point scale factor

Applies when

the point in the pole's hemisphere (not at the pole)
30 <= |lat_ts| <= 89.99 in the same hemisphere

Example

m.ask("Polar stereographic forward: lat = 75, lon = -40, pole = north, lat_ts = 70, lon0 = -45, ellipsoid = WGS84")

→ x = 142402 m, y = -1.62766e+06 m, k = 0.986664

Sinusoidal equal-area projection (ellipsoidal)geo_sinusoidal_forward

x and y of a point in the ellipsoidal sinusoidal (Sanson-Flamsteed) equal-area projection, as used for MODIS land grids.

x = \frac{a\cos\varphi}{\sqrt{1-e^2\sin^2\varphi}}(\lambda - \lambda_0),\quad y = m(\varphi)

Inputs

lat latitude (decimal degrees)
lon longitude (decimal degrees)
lon0 central meridian (decimal degrees)
ellipsoid reference ellipsoid choice

Returns

x easting [m]
y northing [m]

Applies when

|latitude| <= 89.9

Example

m.ask("Sinusoidal forward: lat = 45, lon = 10, lon0 = 0, ellipsoid = WGS84")

→ x = 788468 m, y = 4.98494e+06 m

Tissot's indicatrix: scale factors and distortion of a projection at a pointgeo_tissot_indicatrix

The meridian and parallel scale factors h and k, the angle between meridian and parallel, the semi-axes a and b of Tissot's indicatrix, the maximum angular distortion and the areal scale of a world projection of the sphere at a point.

a' = \sqrt{h^2+k^2+2hk\sin\theta'},\ b' = \sqrt{h^2+k^2-2hk\sin\theta'},\ a = \tfrac{a'+b'}{2},\ b = \tfrac{a'-b'}{2},\ \omega = 2\arcsin\frac{a-b}{a+b},\ s = hk\sin\theta'

Inputs

projection projection (spherical form, central meridian 0) choice
lat latitude (decimal degrees)
lon longitude measured from the central meridian (decimal degrees)

Returns

h scale factor along the meridian
k scale factor along the parallel
theta_prime angle between the projected meridian and parallel (degrees, at most 90)
a maximum scale factor (semi-major axis)
b minimum scale factor (semi-minor axis)
omega maximum angular distortion (degrees)
s areal scale factor

Applies when

|latitude| <= 80

Example

m.ask("Tissot indicatrix: projection = mercator, lat = 60, lon = 0")

→ h = 2, k = 2, theta_prime = 90, a = 2, b = 2, omega = 0, s = 4

Transverse Mercator projection: coordinates, scale factor and grid convergencegeo_transverse_mercator_forward

Easting and northing of a point in a transverse Mercator projection with a given origin, central-meridian scale factor and false origin on the WGS84 or GRS80 ellipsoid, with the point scale factor and the grid convergence.

\xi = \xi' + \sum_{j=1}^{6}\alpha_j\sin2j\xi'\cosh2j\eta',\quad \eta = \eta' + \sum_{j=1}^{6}\alpha_j\cos2j\xi'\sinh2j\eta',\quad x = FE + k_0A\eta,\ y = FN + k_0A(\xi - \xi_0)

Inputs

lat latitude (decimal degrees)
lon longitude (decimal degrees)
lat0 latitude of origin (decimal degrees)
lon0 central meridian (decimal degrees)
k0 scale factor on the central meridian
ellipsoid reference ellipsoid choice
FE false easting (default 0) [m] optional
FN false northing (default 0) [m] optional

Returns

x easting [m]
y northing [m]
k point scale factor
gamma grid convergence (degrees, positive east of the central meridian in the north)

Applies when

within 10 degrees of longitude of the central meridian
latitude within 89 degrees of the equator

Example

m.ask("Transverse Mercator forward: lat = 52.5, lon = 13.4, lat0 = 0, lon0 = 15, k0 = 0.9996, ellipsoid = GRS80, FE = 500000 m, FN = 0 m")

→ x = 391391 m, y = 5.81786e+06 m, k = 0.999745, gamma = -1.26949

UTM coordinates to latitude and longitudegeo_utm_inverse

Latitude and longitude on the WGS84 or GRS80 ellipsoid of a point given by its UTM easting, northing, zone and hemisphere.

\xi' = \xi - \sum_{j=1}^{6}\beta_j\sin2j\xi\cosh2j\eta,\quad \eta' = \eta - \sum_{j=1}^{6}\beta_j\cos2j\xi\sinh2j\eta

Inputs

easting UTM easting [m]
northing UTM northing [m]
zone UTM zone number integer
hemisphere hemisphere of the zone (N: false northing 0; S: false northing 10000 km) choice
ellipsoid reference ellipsoid choice

Returns

lat latitude (degrees)
lon longitude (degrees)

Applies when

100 km <= easting <= 900 km
0 <= northing <= 10000 km

Example

m.ask("UTM inverse: easting = 699316.0 m, northing = 5710164.0 m, zone = 30, hemisphere = N, ellipsoid = WGS84")

→ lat = 51.5074, lon = -0.127803

Web Mercator (EPSG:3857) to latitude and longitudegeo_web_mercator_inverse

WGS84 latitude and longitude of a point given by its Web Mercator (pseudo-Mercator) x and y.

\varphi = \arctan\sinh\frac{y}{R},\quad \lambda = \frac{x}{R}

Inputs

x Web Mercator x [m]
y Web Mercator y [m]

Returns

lat latitude (degrees)
lon longitude (degrees)

Applies when

inside the web-map square

Example

m.ask("Web Mercator inverse: x = -14471.5 m, y = 6710219.1 m")

→ lat = 51.5, lon = -0.13

Web-map tile (slippy map / XYZ) and quadkey containing a pointgeo_web_map_tile

The x and y indices of the web-map (XYZ / slippy map) tile containing a WGS84 point at a zoom level, and the tile's quadkey.

x_t = \left\lfloor 2^z\,\frac{\lambda + 180^\circ}{360^\circ}\right\rfloor,\quad y_t = \left\lfloor 2^z\left(\frac12 - \frac{\operatorname{arsinh}(\tan\varphi)}{2\pi}\right)\right\rfloor

Inputs

lat latitude (decimal degrees)
lon longitude (decimal degrees)
zoom zoom level (0-24) integer

Returns

tile_x tile column (from 180 W)
tile_y tile row (from the north edge)
quadkey quadkey of the tile

Applies when

|latitude| <= 85.0511
zoom 0-24

Example

m.ask("Web map tile: lat = 51.5, lon = -0.13, zoom = 12")

→ tile_x = 2046, tile_y = 1362, quadkey = 031313131130

Network accessibility 6

Gravity-based (Hansen) accessibility on a road networkgeo_gravity_accessibility

The Hansen gravity accessibility of every node of a network: the sum of the opportunities (jobs, services) at the other nodes discounted by a power or negative-exponential function of the shortest-path distance.

A_i = \sum_{j\ne i} O_j\,f(d_{ij}),\quad f(d) = d^{-\beta}\ \text{or}\ e^{-\beta d}

Inputs

edges undirected network as weighted edges between nodes numbered 0 .. n-1 edges
opportunities opportunities at nodes 0 .. n-1 vector
beta distance-decay parameter
impedance impedance function d^-beta or exp(-beta d) choice

Returns

accessibility accessibility of each node

Applies when

positive edge lengths
a node's own opportunities are excluded; unreachable nodes contribute 0

Example

m.ask("Gravity accessibility: edges = [0-1 2; 1-2 3; 0-2 6; 2-3 1], opportunities = [100, 50, 80, 20], beta = 1, impedance = power")

→ accessibility = [44.3333333333, 81.6666666667, 56.6666666667, 109.166666667]

Network connectivity indices: cyclomatic number, alpha, beta and gammageo_network_connectivity_indices

The graph-theoretic connectivity measures of a (planar) transport network: vertices, edges, components, cyclomatic number and Kansky's alpha, beta and gamma indices.

\mu = e - v + p,\quad \alpha = \frac{e - v + p}{2v - 5},\quad \beta = \frac{e}{v},\quad \gamma = \frac{e}{3(v-2)}

Inputs

edges network edges between nodes numbered from 0 (lengths, if given, are ignored) edges

Returns

vertices number of vertices v
edges number of distinct edges e
components connected components p
cyclomatic cyclomatic number e - v + p
alpha (e - v + p)/(2v - 5)
beta e / v
gamma e / (3(v - 2))

Applies when

at least 3 vertices
planar formulas for alpha and gamma

Example

m.ask("Network connectivity indices: edges = [0-1; 1-2; 2-3; 3-0; 0-2]")

→ vertices = 4, edges = 5, components = 1, cyclomatic = 2, alpha = 0.666667, beta = 1.25, gamma = 0.833333

Network detour index (circuity) between two nodesgeo_detour_index

The shortest network distance between two nodes of a street network (edge lengths from the node coordinates), the straight-line distance and their ratio, the detour index or circuity.

DI = \frac{d_{net}(o, t)}{\lVert\mathbf x_o - \mathbf x_t\rVert}

Inputs

coords x and y of nodes 0 .. n-1 (n x 2, projected coordinates) matrix
edges street segments between nodes (lengths are computed from the coordinates) edges
origin origin node
destination destination node

Returns

network_distance shortest network distance
euclidean_distance straight-line distance
detour_index network / straight-line distance (>= 1)

Applies when

straight edges between the node coordinates
destination reachable

Example

m.ask("Detour index: coords = [[0, 0], [3, 0], [3, 4], [0, 4]], edges = [0-1; 1-2; 2-3], origin = 0, destination = 2")

→ network_distance = 7, euclidean_distance = 5, detour_index = 1.4

Network service area (isochrone) of a facilitygeo_service_area

The nodes of a road network within a travel-cost cutoff of a facility (the service area or isochrone), their number and the largest distance reached.

SA(s, c) = \{v : d(s, v) \le c\}

Inputs

edges undirected network as weighted edges between nodes numbered 0 .. n-1 edges
source facility node
cutoff maximum network distance or travel time

Returns

reachable nodes within the cutoff (ascending)
count number of nodes reached (the facility included)
farthest largest network distance within the cutoff

Applies when

positive edge costs
distance <= cutoff counts as reached

Example

m.ask("Service area: edges = [0-1 2; 1-2 3; 0-2 6; 2-3 1; 3-4 4], source = 0, cutoff = 6")

→ reachable = [0.0, 1.0, 2.0, 3.0], count = 4, farthest = 6

Shimbel accessibility index of every nodegeo_shimbel_index

The Shimbel (total shortest-path distance) index of every node of a connected network, the network total (dispersion) and the most accessible node (smallest index).

A_i = \sum_j d_{ij},\quad D(G) = \sum_i A_i

Inputs

edges undirected network as weighted edges between nodes numbered 0 .. n-1 edges

Returns

shimbel sum of shortest-path distances from each node
total network dispersion (sum over all nodes)
most_accessible node with the smallest Shimbel index (lowest number on ties)

Applies when

connected network with nodes 0 .. n-1

Example

m.ask("Shimbel index: edges = [0-1 2; 1-2 3; 0-2 6; 2-3 1]")

→ shimbel = [13.0, 9.0, 9.0, 11.0], total = 42, most_accessible = 1

Two-step floating catchment area (2SFCA) accessibilitygeo_two_step_floating_catchment

The 2SFCA spatial accessibility of every node to services (doctors, schools): supply-to-population ratios within each site's catchment, summed over the sites within reach of each node.

R_j = \frac{S_j}{\sum_{k: d_{kj}\le d_0}P_k},\quad A_i = \sum_{j: d_{ij}\le d_0}R_j

Inputs

edges undirected network as weighted edges between nodes numbered 0 .. n-1 edges
population population at nodes 0 .. n-1 vector
supply_nodes nodes with a service site vector
capacity capacity of each site (e.g. physicians) vector
threshold catchment threshold (network distance)

Returns

accessibility accessibility of each node (supply per person)
ratios supply-to-population ratio of each site

Applies when

a catchment includes nodes at network distance <= threshold
a site with no population in reach has ratio 0

Example

m.ask("Two-step floating catchment: edges = [0-1 2; 1-2 3; 0-2 6; 2-3 1], population = [1000, 500, 800, 200], supply_nodes = [1, 3], capacity = [5, 2], threshold = 3")

→ accessibility = [0.00217391304348, 0.00217391304348, 0.00417391304348, 0.002], ratios = [0.00217391304348, 0.002]

Planar vector analysis 3

Discrete Frechet distance between two polylinesgeo_discrete_frechet_distance

The discrete Frechet ('dog-walking') distance between two polylines (GPS tracks, routes) measured over their vertices.

c_{i,j} = \max\left(d(P_i, Q_j),\ \min(c_{i-1,j},\ c_{i-1,j-1},\ c_{i,j-1})\right)

Inputs

x1 x coordinates of the first line vector
y1 y coordinates of the first line vector
x2 x coordinates of the second line vector
y2 y coordinates of the second line vector

Returns

distance discrete Frechet distance (coordinate units)

Applies when

at least 2 vertices per line

Example

m.ask("Frechet distance: x1 = [0, 1, 2], y1 = [0, 0, 0], x2 = [0, 2], y2 = [1, 1]")

→ distance = 1.41421

Douglas-Peucker line simplificationgeo_douglas_peucker_simplification

The vertices of a polyline kept by the Ramer-Douglas-Peucker simplification at a distance tolerance (indices from 0).

d_{max} = \max_k \operatorname{dist}(P_k, \overline{P_iP_j}) > \varepsilon \Rightarrow \text{keep } P_k

Inputs

x x coordinates of the vertices in order vector
y y coordinates of the vertices vector
tolerance distance tolerance (same unit as the coordinates)

Returns

kept indices of the kept vertices (from 0, in order)
n_kept number of kept vertices

Applies when

distinct vertices
tolerance > 0

Example

m.ask("Douglas-Peucker simplification: x = [0, 1, 2, 3, 4, 5, 6, 7], y = [0, 0.1, -0.1, 5, 6, 7, 8.1, 9], tolerance = 1")

→ kept = [0.0, 2.0, 3.0, 7.0], n_kept = 4

Shape compactness indices of a polygon (Polsby-Popper, Schwartzberg, convex hull, Reock)geo_polygon_compactness

Area, perimeter and the standard compactness measures of a polygon in projected coordinates (used for districts and land parcels): Polsby-Popper, Schwartzberg, convex-hull ratio and Reock.

PP = \frac{4\pi A}{P^2},\quad S = \frac{P}{2\sqrt{\pi A}},\quad C_H = \frac{A}{A_{hull}},\quad R = \frac{A}{\pi r_{mec}^2}

Inputs

x x coordinates of the vertices in order (not repeating the first) vector
y y coordinates of the vertices vector

Returns

area area (square coordinate units)
perimeter perimeter (coordinate units)
polsby_popper 4 pi A / P^2
schwartzberg P / (2 sqrt(pi A))
convex_hull_ratio area / convex hull area
reock area / area of the minimum enclosing circle

Applies when

a simple polygon with distinct vertices

Example

m.ask("Polygon compactness: x = [0, 4, 4, 0], y = [0, 0, 3, 3]")

→ area = 12, perimeter = 14, polsby_popper = 0.76937, schwartzberg = 1.14007, convex_hull_ratio = 1, reock = 0.611155

Point pattern analysis 7

Kernel density estimate of a point pattern at a locationgeo_kernel_density

The two-dimensional kernel density (per unit area, integrating to 1) and the intensity (events per unit area) of a point pattern at a location, with a stated kernel and bandwidth.

\hat f(\mathbf p) = \frac1n\sum_i K_h(\lVert\mathbf p - \mathbf x_i\rVert),\quad K^{gauss}_h = \frac{e^{-u^2/2}}{2\pi h^2},\ K^{epa}_h = \frac{2(1-u^2)}{\pi h^2},\ K^{unif}_h = \frac{1}{\pi h^2}

Inputs

x x coordinates of the events vector
y y coordinates of the events vector
px x of the location
py y of the location
bandwidth bandwidth h (coordinate units)
kernel kernel (gaussian with sd h; epanechnikov and uniform with radius h) choice

Returns

density density (per square coordinate unit, integrates to 1)
intensity intensity = n x density (events per square unit)

Applies when

bandwidth > 0

Example

m.ask("Kernel density: x = [0, 1, 2, 5], y = [0, 1, 0, 5], px = 1, py = 0.5, bandwidth = 2, kernel = gaussian")

→ density = 0.0267648, intensity = 0.107059

Mean centre and standard distance of a point patterngeo_mean_center_standard_distance

The mean centre of a set of points and the standard distance (root-mean-square distance from the mean centre).

SD = \sqrt{\frac{\sum_i (x_i - \bar x)^2 + (y_i - \bar y)^2}{n}}

Inputs

x x coordinates vector
y y coordinates vector

Returns

mean_x x of the mean centre
mean_y y of the mean centre
std_distance standard distance

Applies when

population form (divide by n)

Example

m.ask("Mean center and standard distance: x = [0, 4, 4, 0], y = [0, 0, 3, 3]")

→ mean_x = 2, mean_y = 1.5, std_distance = 2.5

Nearest-neighbour index (Clark-Evans ratio) of a point patterngeo_nearest_neighbour_index

The mean observed nearest-neighbour distance of a point pattern, the mean expected under complete spatial randomness in a study area, their ratio R (R < 1 clustered, R > 1 dispersed), the z-score and two-sided p-value.

R = \frac{\bar d_{obs}}{0.5/\sqrt{n/A}},\quad z = \frac{\bar d_{obs} - \bar d_{exp}}{0.26136/\sqrt{n^2/A}}

Inputs

x x coordinates of the points vector
y y coordinates of the points vector
area area of the study region (square coordinate units)

Returns

d_obs mean observed nearest-neighbour distance
d_exp expected mean distance 0.5/sqrt(n/A)
R nearest-neighbour index d_obs/d_exp
z z-score
p_value two-sided p-value

Applies when

distinct points
no edge correction

Example

m.ask("Nearest neighbour index: x = [1, 2, 8, 9, 5], y = [1, 2, 8, 9, 5], area = 100")

→ d_obs = 1.9799, d_exp = 2.23607, R = 0.885438, z = -0.490069, p_value = 0.624085

Quadrat analysis: chi-square test of complete spatial randomnessgeo_quadrat_test

Counts of points in an nx x ny grid of quadrats over the bounding rectangle of a point pattern, the chi-square statistic against complete spatial randomness, its p-value and the variance-to-mean ratio.

\chi^2 = \sum_{c}\frac{(O_c - E)^2}{E},\quad E = \frac{n}{n_xn_y},\quad VMR = \frac{\chi^2}{n_xn_y - 1}

Inputs

x x coordinates of the points vector
y y coordinates of the points vector
nx number of quadrat columns integer
ny number of quadrat rows integer

Returns

chi2 chi-square statistic
df degrees of freedom (cells - 1)
p_value upper-tail p-value
vmr variance-to-mean ratio of the counts (chi2/df)
counts count of each quadrat, row by row from the south-west

Applies when

the window is the bounding rectangle of the points
cells include their west and south edges; the last column and row include the east and north edges

Example

m.ask("Quadrat test: x = [1, 2, 3, 8, 9, 2, 1, 7, 8, 9], y = [1, 1, 2, 8, 9, 8, 9, 1, 2, 1], nx = 2, ny = 2")

→ chi2 = 0.4, df = 3, p_value = 0.940242, vmr = 0.133333, counts = [3.0, 3.0, 2.0, 2.0]

Ripley's K and L functions of a point patterngeo_ripley_k

Ripley's K function and its variance-stabilized L function at given distances for a point pattern, without edge correction, in the bounding rectangle of the points.

\hat K(d) = \frac{A}{n^2}\sum_{i\ne j}\mathbf{1}(d_{ij} < d),\quad L(d) = \sqrt{K(d)/\pi}

Inputs

x x coordinates of the points vector
y y coordinates of the points vector
d distances at which to evaluate K vector

Returns

K K(d) at each distance
L L(d) = sqrt(K/pi) at each distance
area area of the bounding rectangle

Applies when

window = bounding rectangle of the points
pairs closer than d (strictly) are counted; no edge correction

Example

m.ask("Ripley's K: x = [0, 1, 2, 8, 9, 10], y = [0, 1, 0, 9, 10, 9], d = [1.5, 3, 12]")

→ K = [22.2222222222, 33.3333333333, 44.4444444444], L = [2.65961520268, 3.25735007935, 3.76126389032], area = 100

Standard deviational ellipse of a point patterngeo_standard_deviational_ellipse

The semi-axes of the standard deviational ellipse of a point pattern (CrimeStat correction sqrt(2) sqrt(n/(n-2))) and the orientation of its major axis.

s_{x,y} = \sqrt{\frac{2n}{n-2}\lambda_{1,2}},\quad \lambda_{1,2} = \frac{s_{xx}+s_{yy}}{2} \pm \sqrt{\frac{(s_{xx}-s_{yy})^2}{4} + s_{xy}^2}

Inputs

x x coordinates vector
y y coordinates vector

Returns

sx semi-major axis
sy semi-minor axis
orientation bearing of the major axis, degrees clockwise from north (y axis), in [0, 180)

Applies when

n >= 3, points neither collinear nor circularly symmetric

Example

m.ask("Standard deviational ellipse: x = [66.22, 22.52, 31.01, 9.47, 30.78, 75.21, 79.26, 8.23, 98.73, 89.78, 65.19], y = [32.54, 22.39, 81.21, 31.02, 60.1, 58.93, 7.68, 39.93, 77.17, 42.53, 92.08]")

→ sx = 50.1303, sy = 37.9522, orientation = 70.1437

Weighted mean centre and weighted standard distancegeo_weighted_mean_center

The weighted mean centre of a set of points (e.g. population-weighted centroid) and the weighted standard distance about it.

\bar x_w = \frac{\sum w_ix_i}{\sum w_i},\quad SD_w = \sqrt{\frac{\sum w_i\left[(x_i - \bar x_w)^2 + (y_i - \bar y_w)^2\right]}{\sum w_i}}

Inputs

x x coordinates vector
y y coordinates vector
w non-negative weights vector

Returns

wmean_x x of the weighted mean centre
wmean_y y of the weighted mean centre
weighted_std_distance weighted standard distance

Applies when

non-negative weights, not all zero

Example

m.ask("Weighted mean center: x = [0, 10, 20], y = [0, 0, 10], w = [100, 50, 50]")

→ wmean_x = 7.5, wmean_y = 2.5, weighted_std_distance = 9.35414

Remote sensing 6

At-sensor brightness temperature from thermal radiancegeo_brightness_temperature

The at-sensor brightness temperature of a thermal-infrared band from its spectral radiance and the band's calibration constants K1 and K2 (as in Landsat metadata).

T_B = \frac{K_2}{\ln\left(\frac{K_1}{L_\lambda} + 1\right)}

Inputs

radiance spectral radiance (W/(m^2 sr um))
K1 band constant K1 (W/(m^2 sr um))
K2 band constant K2 (kelvin)

Returns

T brightness temperature [K]

Applies when

positive radiance and constants

Example

m.ask("Brightness temperature: radiance = 10.5, K1 = 774.8853, K2 = 1321.0789")

→ T = 306.174 K

Ground sample distance and swath of a nadir-looking camerageo_ground_sample_distance

The ground sample distance (pixel footprint) and swath width of a nadir-looking frame camera or scanner from its pixel pitch, focal length, flying height above ground and number of pixels across track.

GSD = \frac{p\,H}{f},\quad W = GSD \times N

Inputs

pixel_pitch detector pixel pitch [m]
focal_length focal length [m]
altitude height above the ground [m]
pixels number of pixels across track integer

Returns

gsd ground sample distance [m]
swath swath width [m]

Applies when

flat ground, nadir view

Example

m.ask("Ground sample distance: pixel_pitch = 0.0024 mm, focal_length = 8.8 mm, altitude = 120 m, pixels = 5472")

→ gsd = 0.0327273 m, swath = 179.084 m

Thematic map accuracy: overall accuracy, kappa, producer's and user's accuracygeo_classification_accuracy

The accuracy assessment of a classified map from its error (confusion) matrix: overall accuracy, Cohen's kappa, and the producer's and user's accuracy of each class.

OA = \frac{\sum_i n_{ii}}{N},\quad \kappa = \frac{p_o - p_e}{1 - p_e},\quad PA_j = \frac{n_{jj}}{n_{+j}},\quad UA_i = \frac{n_{ii}}{n_{i+}}

Inputs

matrix confusion matrix of sample counts: rows = map (classified) classes, columns = reference classes matrix

Returns

overall overall accuracy
kappa Cohen's kappa
producers producer's accuracy of each class (1 - omission error)
users user's accuracy of each class (1 - commission error)

Applies when

square matrix of whole counts
every class sampled in the map and the reference

Example

m.ask("Classification accuracy: matrix = [[65, 4, 22, 24], [6, 81, 5, 8], [0, 11, 85, 19], [4, 7, 3, 90]]")

→ overall = 0.739631, kappa = 0.653516, producers = [0.866666666667, 0.78640776699, 0.739130434783, 0.63829787234], users = [0.565217391304, 0.81, 0.739130434783, 0.865384615385]

Top-of-atmosphere reflectance from at-sensor radiancegeo_toa_reflectance

The top-of-atmosphere (planetary) reflectance of a pixel from its at-sensor spectral radiance, the band's mean exoatmospheric solar irradiance ESUN, the day of year and the sun elevation, with the Earth-Sun distance.

\rho_{TOA} = \frac{\pi L_\lambda d^2}{ESUN_\lambda\cos\theta_s},\quad \frac{1}{d^2} = 1.00011 + 0.034221\cos\Gamma + 0.00128\sin\Gamma + 0.000719\cos2\Gamma + 0.000077\sin2\Gamma

Inputs

radiance at-sensor spectral radiance (W/(m^2 sr um))
esun band mean exoatmospheric solar irradiance ESUN (W/(m^2 um))
doy day of year (1-366) integer
sun_elevation sun elevation angle (decimal degrees)

Returns

earth_sun_distance Earth-Sun distance (astronomical units)
reflectance top-of-atmosphere reflectance

Applies when

sun above the horizon

Example

m.ask("TOA reflectance: radiance = 85, esun = 1536, doy = 172, sun_elevation = 60")

→ earth_sun_distance = 1.01669, reflectance = 0.207502

Vegetation indices: NDVI, EVI, SAVI, OSAVI, MSAVI and EVI2 from band reflectancesgeo_vegetation_indices

The common spectral vegetation indices of a pixel from its red, near-infrared and blue surface reflectances, with a stated soil adjustment factor L for SAVI (EVI uses its standard coefficients G = 2.5, C1 = 6, C2 = 7.5, L = 1).

NDVI = \frac{\rho_{N} - \rho_R}{\rho_N + \rho_R},\ EVI = 2.5\frac{\rho_N - \rho_R}{\rho_N + 6\rho_R - 7.5\rho_B + 1},\ SAVI = (1+L)\frac{\rho_N - \rho_R}{\rho_N + \rho_R + L}

Inputs

red red reflectance (0-1)
nir near-infrared reflectance (0-1)
blue blue reflectance (0-1)
L SAVI soil adjustment factor (0.5 is usual)

Returns

NDVI (NIR - R)/(NIR + R)
EVI 2.5 (NIR - R)/(NIR + 6R - 7.5B + 1)
SAVI (1 + L)(NIR - R)/(NIR + R + L)
OSAVI (NIR - R)/(NIR + R + 0.16)
MSAVI (2NIR + 1 - sqrt((2NIR + 1)^2 - 8(NIR - R)))/2
EVI2 2.5 (NIR - R)/(NIR + 2.4R + 1)

Applies when

reflectances between 0 and 1
positive denominators

Example

m.ask("Vegetation indices: red = 0.08, nir = 0.45, blue = 0.04, L = 0.5")

→ NDVI = 0.698113, EVI = 0.567485, SAVI = 0.538835, OSAVI = 0.536232, MSAVI = 0.546887, EVI2 = 0.563337

Water, moisture, burn and built-up indices: NDWI, MNDWI, NDMI, NBR and NDBIgeo_water_burn_built_indices

Normalized-difference indices for open water (NDWI, MNDWI), vegetation moisture (NDMI), burn severity (NBR) and built-up land (NDBI) from green, near-infrared and two shortwave-infrared reflectances.

ND(a, b) = \frac{\rho_a - \rho_b}{\rho_a + \rho_b}

Inputs

green green reflectance (0-1)
nir near-infrared reflectance (0-1)
swir1 shortwave infrared 1 (about 1.6 um) reflectance
swir2 shortwave infrared 2 (about 2.2 um) reflectance

Returns

NDWI (G - NIR)/(G + NIR)
MNDWI (G - SWIR1)/(G + SWIR1)
NDMI (NIR - SWIR1)/(NIR + SWIR1)
NBR (NIR - SWIR2)/(NIR + SWIR2)
NDBI (SWIR1 - NIR)/(SWIR1 + NIR)

Applies when

reflectances between 0 and 1, band sums non-zero

Example

m.ask("Water indices: green = 0.08, nir = 0.03, swir1 = 0.02, swir2 = 0.01")

→ NDWI = 0.454545, MNDWI = 0.6, NDMI = 0.2, NBR = 0.5, NDBI = -0.2

Spatial statistics 10

Bivariate Moran's I (spatial cross-correlation)geo_bivariate_morans_i

The bivariate Moran's I between a variable x and the spatial lag of a second variable y with row-standardized weights (how x at a location relates to y around it).

I_{xy} = \frac{\sum_i z_{x,i}\sum_j w_{ij}z_{y,j}}{n-1}

Inputs

x first variable at the n locations vector
y second variable at the n locations vector
W spatial weights matrix (row-standardized internally) matrix

Returns

I_xy bivariate Moran's I

Applies when

n >= 4, neither variable constant
no observation without neighbours

Example

m.ask("Bivariate Moran's I: x = [8, 6, 6, 3, 2, 1], y = [7, 7, 5, 4, 1, 2], W = [[0, 1, 0, 0, 0, 0], [1, 0, 1, 0, 0, 0], [0, 1, 0, 1, 0, 0], [0, 0, 1, 0, 1, 0], [0, 0, 0, 1, 0, 1], [0, 0, 0, 0, 1, 0]]")

→ I_xy = 0.89176

Distance-band spatial weightsgeo_distance_band_weights

The binary distance-band weights matrix of a set of points (neighbours within a threshold distance), the number of neighbours of each point, the total weight and the number of islands.

w_{ij} = \begin{cases}1 & d_{ij} \le \delta,\ i \ne j\\ 0 & \text{otherwise}\end{cases}

Inputs

x x coordinates vector
y y coordinates vector
threshold distance threshold (coordinate units)

Returns

W weights matrix (n x n, 0/1)
neighbours number of neighbours of each point
s0 sum of all weights
islands number of points without neighbours

Applies when

threshold > 0

Example

m.ask("Distance band weights: x = [0, 1, 3, 6, 10], y = [0, 0, 0, 0, 0], threshold = 2.5")

→ W = [[0.0, 1.0, 0.0, 0.0, 0.0], [1.0, 0.0, 1.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0]], neighbours = [1.0, 2.0, 1.0, 0.0, 0.0], s0 = 4, islands = 2

Geary's C spatial autocorrelation with normal and randomization inferencegeo_gearys_c

Geary's contiguity ratio C of a variable with a given spatial weights matrix, its variance under normality and randomization, and the z-scores (C < 1 positive, C > 1 negative autocorrelation).

C = \frac{(n-1)\sum_i\sum_j w_{ij}(y_i - y_j)^2}{2S_0\sum_i (y_i - \bar y)^2},\quad V_N[C] = \frac{(2S_1 + S_2)(n-1) - 4S_0^2}{2(n+1)S_0^2}

Inputs

y values at the n locations vector
W spatial weights matrix (n x n, zero diagonal, row i lists the weights of i's neighbours) matrix
transform row-standardize the weights (row) or use them as given (none) choice

Returns

C Geary's C
VC_norm variance under normality
z_norm z-score under normality
VC_rand variance under randomization
z_rand z-score under randomization

Applies when

n >= 4, values not all equal
no observation without neighbours

Example

m.ask("Geary's C: y = [8, 6, 6, 3, 2, 1], W = [[0, 1, 0, 0, 0, 0], [1, 0, 1, 0, 0, 0], [0, 1, 0, 1, 0, 0], [0, 0, 1, 0, 1, 0], [0, 0, 0, 1, 0, 1], [0, 0, 0, 0, 1, 0]], transform = row")

→ C = 0.195312, VC_norm = 0.111111, z_norm = -2.41406, VC_rand = 0.132231, z_rand = -2.21289

Getis-Ord general G (global high/low clustering)geo_getis_ord_general_g

The Getis-Ord general G statistic of a positive variable with binary spatial weights, its expectation and variance under randomization and the z-score (G above its expectation: high values cluster).

G = \frac{\sum_i\sum_{j\neq i} w_{ij}y_iy_j}{\sum_i\sum_{j\neq i}y_iy_j},\quad E[G] = \frac{S_0}{n(n-1)}

Inputs

y positive values at the n locations vector
W binary spatial weights matrix (0/1, zero diagonal) matrix

Returns

G general G
EG expected G = S0/(n(n-1))
VG variance of G under randomization
z_norm z-score

Applies when

positive values, n >= 4
binary weights, no observation without neighbours

Example

m.ask("Getis-Ord general G: y = [10, 12, 11, 2, 3, 2], W = [[0, 1, 1, 0, 0, 0], [1, 0, 1, 0, 0, 0], [1, 1, 0, 1, 0, 0], [0, 0, 1, 0, 1, 1], [0, 0, 0, 1, 0, 1], [0, 0, 0, 1, 1, 0]]")

→ G = 0.656814, EG = 0.466667, VG = 0.00892117, z_norm = 2.01317

Global Moran's I spatial autocorrelation with normal and randomization inferencegeo_morans_i

Global Moran's I of a variable observed at n locations with a given spatial weights matrix, its expectation and variance under normality and under randomization, z-scores and two-sided p-values.

I = \frac{n}{S_0}\frac{\sum_i\sum_j w_{ij}z_iz_j}{\sum_i z_i^2},\quad E[I] = -\frac{1}{n-1},\quad V_N[I] = \frac{n^2S_1 - nS_2 + 3S_0^2}{(n^2-1)S_0^2} - E[I]^2

Inputs

y values at the n locations vector
W spatial weights matrix (n x n, zero diagonal, row i lists the weights of i's neighbours) matrix
transform row-standardize the weights (row) or use them as given (none) choice

Returns

I Moran's I
EI expected I = -1/(n-1)
VI_norm variance under normality
z_norm z-score under normality
p_norm two-sided p-value under normality
VI_rand variance under randomization
z_rand z-score under randomization
p_rand two-sided p-value under randomization

Applies when

n >= 4, values not all equal
no observation without neighbours

Example

m.ask("Moran's I: y = [8, 6, 6, 3, 2, 1], W = [[0, 1, 0, 0, 0, 0], [1, 0, 1, 0, 0, 0], [0, 1, 0, 1, 0, 0], [0, 0, 1, 0, 1, 0], [0, 0, 0, 1, 0, 1], [0, 0, 0, 0, 1, 0]], transform = row")

→ I = 0.65625, EI = -0.2, VI_norm = 0.140952, z_norm = 2.28068, p_norm = 0.0225674, VI_rand = 0.181793, z_rand = 2.00823, p_rand = 0.0446193

Join count statistics of a binary map (BB, BW, WW)geo_join_counts

The numbers of black-black, black-white and white-white joins between neighbouring areas of a binary (presence/absence) map, and the total number of joins.

BB = \tfrac12\sum_{i,j} w_{ij}y_iy_j,\ BW = \tfrac12\sum_{i,j} w_{ij}(y_i - y_j)^2,\ WW = \tfrac12\sum_{i,j}w_{ij}(1-y_i)(1-y_j)

Inputs

y binary variable at the n areas (1 black, 0 white) vector
W binary symmetric contiguity matrix (0/1, zero diagonal) matrix

Returns

BB black-black joins
BW black-white joins
WW white-white joins
J total joins

Applies when

binary symmetric weights

Example

m.ask("Join count statistics: y = [1, 1, 0, 0, 1, 0], W = [[0, 1, 1, 0, 0, 0], [1, 0, 1, 1, 0, 0], [1, 1, 0, 0, 1, 0], [0, 1, 0, 0, 1, 1], [0, 0, 1, 1, 0, 1], [0, 0, 0, 1, 1, 0]]")

→ BB = 1, BW = 6, WW = 1, J = 8

k-nearest-neighbour spatial weightsgeo_knn_weights

The binary k-nearest-neighbour weights matrix of a set of points (row i marks the k nearest points to point i), its total weight and whether it is symmetric.

w_{ij} = 1 \iff j \in N_k(i)

Inputs

x x coordinates vector
y y coordinates vector
k number of neighbours integer

Returns

W weights matrix (n x n, 0/1)
s0 sum of all weights (n k)
symmetric whether W is symmetric

Applies when

distinct points
1 <= k < n, no tie at the k-th neighbour

Example

m.ask("k nearest neighbour weights: x = [0, 1, 3.5, 6, 10], y = [0, 0, 0, 0, 0.5], k = 2")

→ W = [[0.0, 1.0, 1.0, 0.0, 0.0], [1.0, 0.0, 1.0, 0.0, 0.0], [0.0, 1.0, 0.0, 1.0, 0.0], [0.0, 0.0, 1.0, 0.0, 1.0], [0.0, 0.0, 1.0, 1.0, 0.0]], s0 = 10, symmetric = False

Local Getis-Ord Gi and Gi* hot-spot statisticsgeo_local_getis_ord_gi

The local Getis-Ord statistic of every location (Gi* includes the location itself, Gi excludes it) with binary weights, and its z-score for hot-spot (positive) and cold-spot (negative) mapping.

G_i^* = \frac{\sum_j w_{ij}y_j}{\sum_j y_j},\quad Z_i = \frac{G_i - W_i/N}{\sqrt{\frac{W_i(N - W_i)}{(N-1)N^2}\frac{s^2}{\bar y^2}}}

Inputs

y positive values at the n locations vector
W binary spatial weights matrix (0/1, zero diagonal) matrix
statistic Gi* (the location counts as its own neighbour) or Gi choice

Returns

Gs local G of each location
Zs z-score of each location

Applies when

positive values, n >= 4
binary weights, no observation without neighbours

Example

m.ask("Local Getis-Ord Gi: y = [10, 12, 11, 2, 3, 2], W = [[0, 1, 1, 0, 0, 0], [1, 0, 1, 0, 0, 0], [1, 1, 0, 1, 0, 0], [0, 0, 1, 0, 1, 1], [0, 0, 0, 1, 0, 1], [0, 0, 0, 1, 1, 0]], statistic = gi_star")

→ Gs = [0.825, 0.825, 0.875, 0.45, 0.175, 0.175], Zs = [2.21006630298, 2.21006630298, 1.50264699011, -1.56275286971, -2.21006630298, -2.21006630298]

Local Moran's I (LISA) and Moran scatterplot quadrantsgeo_local_morans_i

The local Moran statistic of every location with row-standardized weights, and its Moran scatterplot quadrant: 1 high-high, 2 low-high, 3 low-low, 4 high-low.

I_i = (n-1)\frac{z_i\sum_j w_{ij}z_j}{\sum_k z_k^2}

Inputs

y values at the n locations vector
W spatial weights matrix (row-standardized internally) matrix

Returns

Is local Moran's I of each location
quadrant quadrant of each location (1 HH, 2 LH, 3 LL, 4 HL)

Applies when

n >= 4, values not all equal
no observation without neighbours

Example

m.ask("Local Moran's I: y = [8, 6, 6, 3, 2, 1], W = [[0, 1, 0, 0, 0, 0], [1, 0, 1, 0, 0, 0], [0, 1, 0, 1, 0, 0], [0, 0, 1, 0, 1, 0], [0, 0, 0, 1, 0, 1], [0, 0, 0, 0, 1, 0]]")

→ Is = [0.818452380952, 0.595238095238, 0.0372023809524, 0.0595238095238, 0.729166666667, 1.04166666667], quadrant = [1.0, 1.0, 1.0, 3.0, 3.0, 3.0]

Spatial lag of a variablegeo_spatial_lag

The spatial lag W y of a variable: for row-standardized weights, the average of each location's neighbours; otherwise the weighted sum.

[Wy]_i = \sum_j w_{ij}y_j

Inputs

y values at the n locations vector
W spatial weights matrix (n x n, zero diagonal, row i lists the weights of i's neighbours) matrix
transform row-standardize the weights (row) or use them as given (none) choice

Returns

lag spatial lag at each location

Applies when

no observation without neighbours

Example

m.ask("Spatial lag: y = [8, 6, 6, 3], W = [[0, 1, 0, 1], [1, 0, 1, 0], [0, 1, 0, 1], [1, 0, 1, 0]], transform = row")

→ lag = [4.5, 7.0, 4.5, 7.0]

Terrain analysis 8

D8 flow accumulation of a small DEMgeo_d8_flow_accumulation

The number of cells draining through each cell (the cell itself included) under D8 flow directions, and the largest accumulation (the main outlet).

A(c) = 1 + \sum_{u\to c}A(u)

Inputs

dem elevations (row 0 to the north), at most 900 cells matrix
cell_size cell size [m]

Returns

accumulation contributing cells of each cell, itself included (matrix)
max_accumulation largest accumulation

Applies when

D8 directions as in geo_d8_flow_direction

Example

m.ask("Flow accumulation: dem = [[78, 72, 69, 71], [74, 67, 56, 49], [69, 53, 44, 37], [64, 58, 55, 22]], cell_size = 10 m")

→ accumulation = [[1.0, 1.0, 1.0, 1.0], [1.0, 2.0, 2.0, 3.0], [1.0, 4.0, 8.0, 6.0], [1.0, 1.0, 1.0, 16.0]], max_accumulation = 16

D8 flow direction of every cell of a small DEMgeo_d8_flow_direction

The D8 (steepest-descent) flow direction of every cell of a digital elevation model, as ESRI codes (1 E, 2 SE, 4 S, 8 SW, 16 W, 32 NW, 64 N, 128 NE; 0 where no neighbour inside the raster is lower).

d^* = \arg\max_{k}\frac{z_0 - z_k}{\ell_k},\quad \ell_k \in \{\Delta, \sqrt2\Delta\}

Inputs

dem elevations (row 0 to the north), at most 900 cells matrix
cell_size cell size [m]

Returns

direction D8 code of each cell (matrix)

Applies when

square cells; diagonal distance sqrt(2) cells
ties go to the first direction in the order E, SE, S, SW, W, NW, N, NE

Example

m.ask("D8 flow direction: dem = [[78, 72, 69, 71], [74, 67, 56, 49], [69, 53, 44, 37], [64, 58, 55, 22]], cell_size = 10 m")

→ direction = [[2.0, 2.0, 2.0, 4.0], [2.0, 2.0, 2.0, 4.0], [1.0, 1.0, 2.0, 4.0], [128.0, 128.0, 1.0, 0.0]]

Hillshade (analytical relief shading) of a 3x3 elevation windowgeo_hillshade

The illumination (cosine of the solar incidence angle on the terrain) and the 0-255 hillshade value at the centre of a 3 x 3 elevation window for a sun azimuth and altitude.

\cos i = \cos Z\cos S + \sin Z\sin S\cos(A - \alpha),\quad H = 255\max(0, \cos i)

Inputs

window 3 x 3 elevations in metres, row 0 to the north matrix
cell_size cell size [m]
sun_azimuth sun azimuth, clockwise from north (decimal degrees)
sun_altitude sun altitude above the horizon (decimal degrees)

Returns

illumination cos(incidence angle) (negative in self-shadow)
hillshade 255 x max(0, illumination)

Applies when

row 0 is the northern row
no cast shadows

Example

m.ask("Hillshade: window = [[2452.25, 2463, 2473.75], [2450, 2461, 2472], [2447.75, 2459, 2470.25]], cell_size = 5 m, sun_azimuth = 315, sun_altitude = 45")

→ illumination = 0.656099, hillshade = 167.305

Line-of-sight visibility along a terrain profile (with Earth curvature and refraction)geo_line_of_sight_profile

Which points along a terrain profile are visible from an observer at its start, for given observer and target heights above the ground, optionally with the Earth-curvature and refraction correction.

\text{visible}(i) \iff z_j < z_{eye} + (z_i + h_t - z_{eye})\frac{d_j}{d_i}\ \forall\, 0 < j < i,\quad \Delta z_{curv} = \frac{d^2(1-k)}{2R}

Inputs

d distances along the profile in metres, starting at 0 (the observer) vector
z ground elevations in metres vector
observer_height eye height above the ground [m]
target_height target height above the ground [m]
curvature none, or the drop d^2 (1 - k)/(2R) with R = 6371 km and k = 0.13 choice

Returns

visible 1 if the target over each point is visible, else 0 (the observer's own point is 1)
n_visible number of visible points beyond the observer

Applies when

distances increase from 0
positive observer and target heights

Example

m.ask("Line of sight: d = [0, 100, 200, 300, 400, 500], z = [100, 104, 112, 103, 101, 109], observer_height = 2 m, target_height = 2 m, curvature = none")

→ visible = [1.0, 1.0, 1.0, 0.0, 0.0, 0.0], n_visible = 2

Slope and aspect of a 3x3 elevation window (Horn's method)geo_slope_aspect_horn

The slope (degrees and per cent) and aspect (downslope direction, degrees clockwise from north) at the centre cell of a 3 x 3 window of a digital elevation model, by Horn's third-order finite differences.

\frac{\partial z}{\partial x} = \frac{(z_3 + 2z_6 + z_9) - (z_1 + 2z_4 + z_7)}{8\Delta},\ \frac{\partial z}{\partial y} = \frac{(z_1 + 2z_2 + z_3) - (z_7 + 2z_8 + z_9)}{8\Delta},\ S = \arctan\lVert\nabla z\rVert

Inputs

window 3 x 3 elevations in metres, row 0 to the north, column 0 to the west matrix
cell_size cell size (square cells) [m]

Returns

slope_deg slope (degrees)
slope_percent slope (per cent)
aspect_deg aspect: downslope direction clockwise from north (none for a flat window)

Applies when

square cells
row 0 is the northern row

Example

m.ask("Slope and aspect: window = [[37.2, 42, 47.2], [45.2, 50, 55.2], [53.2, 58, 63.2]], cell_size = 10 m")

→ slope_deg = 43.3317, slope_percent = 94.3398, aspect_deg = 327.995

Terrain curvature: total, profile and plan curvature (Zevenbergen-Thorne)geo_terrain_curvature

The total, profile (along the slope) and plan (across the slope) curvature at the centre of a 3 x 3 elevation window from the Zevenbergen-Thorne partial quartic surface.

K = -2(D+E),\quad K_{prof} = \frac{-2(DG^2 + EH^2 + FGH)}{G^2 + H^2},\quad K_{plan} = \frac{2(DH^2 + EG^2 - FGH)}{G^2 + H^2}

Inputs

window 3 x 3 elevations in metres, row 0 to the north matrix
cell_size cell size [m]

Returns

curvature total curvature -2(D + E) (1/m; positive upwardly convex)
profile profile curvature (1/m)
plan plan curvature (1/m)

Applies when

non-zero gradient at the centre for profile and plan curvature

Example

m.ask("Terrain curvature: window = [[10, 12, 14], [9, 10, 13], [7, 9, 11]], cell_size = 1 m")

→ curvature = -3, profile = -1.64, plan = 1.36

Terrain ruggedness index (Riley), topographic position index and local reliefgeo_terrain_ruggedness

Riley's terrain ruggedness index (root of the summed squared elevation differences to the 8 neighbours), the topographic position index (centre minus neighbour mean) and the local relief of a 3 x 3 window.

TRI = \sqrt{\sum_{i=1}^{8}(z_i - z_0)^2},\quad TPI = z_0 - \frac18\sum_{i=1}^{8}z_i

Inputs

window 3 x 3 elevations in metres matrix

Returns

TRI terrain ruggedness index (m)
TPI topographic position index (m)
relief maximum minus minimum elevation (m)

Applies when

Riley et al. (1999) form of the TRI

Example

m.ask("Terrain ruggedness index: window = [[10, 12, 14], [9, 10, 13], [7, 9, 11]]")

→ TRI = 6.40312, TPI = -0.625, relief = 7

Watershed (catchment) of an outlet cell on a small DEMgeo_watershed_delineation

The cells that drain to a chosen outlet cell under D8 flow directions: their number, the catchment area and the catchment mask.

W(o) = \{c : c \to \cdots \to o\}

Inputs

dem elevations (row 0 to the north), at most 900 cells matrix
cell_size cell size [m]
row outlet row (from 0, north) integer
col outlet column (from 0, west) integer

Returns

cells number of cells in the catchment (outlet included)
area catchment area [m^2]
mask 1 for cells in the catchment, 0 elsewhere (matrix)

Applies when

D8 directions as in geo_d8_flow_direction

Example

m.ask("Watershed delineation: dem = [[78, 72, 69, 71], [74, 67, 56, 49], [69, 53, 44, 37], [64, 58, 55, 22]], cell_size = 10 m, row = 3, col = 3")

→ cells = 16, area = 1600 m^2, mask = [[1.0, 1.0, 1.0, 1.0], [1.0, 1.0, 1.0, 1.0], [1.0, 1.0, 1.0, 1.0], [1.0, 1.0, 1.0, 1.0]]

Urban geography 6

Huff model: probabilities of shopping at each store from an origin (network distances)geo_huff_model

The Huff probabilities that a consumer at an origin node patronizes each store, from store attractiveness (floor space) and shortest-path distances on a road network.

P_{ij} = \frac{S_j^{\alpha}d_{ij}^{-\beta}}{\sum_k S_k^{\alpha}d_{ik}^{-\beta}}

Inputs

edges undirected network as weighted edges between nodes numbered 0 .. n-1 edges
origin origin node
stores store nodes vector
attractiveness attractiveness of each store (e.g. floor area) vector
alpha attractiveness exponent
beta distance-decay exponent

Returns

probabilities probability of each store
distances shortest-path distance to each store

Applies when

every store reachable and not at the origin

Example

m.ask("Huff model: edges = [0-1 2; 1-2 3; 0-2 6; 2-3 1; 1-3 5], origin = 0, stores = [1, 2, 3], attractiveness = [5000, 12000, 8000], alpha = 1, beta = 2")

→ probabilities = [0.640295959021, 0.245873648264, 0.113830392715], distances = [2.0, 5.0, 6.0]

Rank-size rule (Zipf) fit of a settlement systemgeo_rank_size_rule

The Zipf exponent q and fitted size of the largest city from a log-log regression of population on rank (P_r = P_1 r^-q), the R^2 of the fit and the two-city primacy ratio.

\ln P_r = \ln P_1 - q\ln r

Inputs

populations populations of the settlements (any order) vector

Returns

q rank-size exponent (1 is Zipf's law)
P1_fit fitted population of rank 1
r_squared R^2 of the log-log fit
primacy largest / second largest population

Applies when

positive populations

Example

m.ask("Rank-size rule: populations = [8400000, 3900000, 2700000, 2300000, 1600000, 1500000, 1400000, 1300000]")

→ q = 0.903102, P1_fit = 7.73565e+06, r_squared = 0.985307, primacy = 2.15385

Residential segregation (evenness): dissimilarity, Gini and Theil entropy indicesgeo_segregation_evenness

The evenness dimension of residential segregation of a group across areal units (census tracts): Duncan's index of dissimilarity, the segregation Gini index and Theil's entropy index H.

D = \frac{\sum_i t_i|p_i - P|}{2TP(1-P)},\quad G = \frac{\sum_i\sum_j t_it_j|p_i - p_j|}{2T^2P(1-P)},\quad H = \sum_i\frac{t_i(E - E_i)}{ET}

Inputs

x group count in each areal unit vector
t total population of each areal unit vector

Returns

dissimilarity index of dissimilarity D
gini Gini segregation index
theil_h Theil's information theory index H

Applies when

0 <= x <= t in each unit
the group present but not the whole population

Example

m.ask("Segregation evenness: x = [100, 200, 50, 800], t = [1000, 1000, 1000, 1000]")

→ dissimilarity = 0.625477, gini = 0.71701, theil_h = 0.364724

Residential segregation (exposure): isolation, interaction and correlation ratiogeo_segregation_exposure

The exposure dimension of residential segregation: the isolation index xPx (probability that a group member meets a member of the same group in their unit), the interaction (exposure) index xPy and the correlation ratio eta^2.

{}_xP^*_x = \sum_i\frac{x_i}{X}\frac{x_i}{t_i},\quad {}_xP^*_y = \sum_i\frac{x_i}{X}\frac{y_i}{t_i},\quad \eta^2 = \frac{{}_xP^*_x - P}{1-P}

Inputs

x group count in each areal unit vector
t total population of each areal unit vector

Returns

isolation isolation index xPx
exposure interaction index xPy
correlation_ratio correlation ratio (xPx - P)/(1 - P)

Applies when

0 <= x <= t in each unit
the group present but not the whole population

Example

m.ask("Segregation exposure: x = [100, 200, 50, 800], t = [1000, 1000, 1000, 1000]")

→ isolation = 0.602174, exposure = 0.397826, correlation_ratio = 0.441648

Spatial interaction (gravity) model fitted by Poisson regressiongeo_gravity_model_poisson_fit

The constant k, the origin and destination mass exponents and the distance-decay exponent of the unconstrained gravity model T = k O^alpha D^beta d^-gamma, fitted to observed flows by Poisson maximum likelihood.

T_{ij} \sim \mathrm{Poisson}(\mu_{ij}),\quad \ln\mu_{ij} = \ln k + \alpha\ln O_i + \beta\ln D_j - \gamma\ln d_{ij}

Inputs

T observed flows between the pairs vector
O origin masses (e.g. populations) vector
D destination masses vector
d distances between the pairs vector

Returns

k constant
alpha origin mass exponent
beta destination mass exponent
gamma distance-decay exponent

Applies when

non-negative flows, positive masses and distances that vary

Example

m.ask("Gravity model fit: T = [120, 45, 300, 80, 15, 60, 200, 33], O = [500000, 200000, 900000, 300000, 100000, 400000, 800000, 250000], D = [300000, 400000, 600000, 200000, 300000, 100000, 500000, 200000], d = [50, 120, 40, 90, 200, 70, 60, 110]")

→ k = 0.000370957, alpha = 0.802892, beta = 0.33475, gamma = 0.509079

Urban population density gradient (Clark's negative exponential model)geo_density_gradient

The central density D0, density gradient b and half-distance of Clark's model D(r) = D0 exp(-b r) fitted by log-linear least squares to densities at distances from the centre, with the R^2.

D(r) = D_0e^{-br},\quad \ln D = \ln D_0 - br

Inputs

r distances from the centre (any unit) vector
D population densities at those distances vector

Returns

D0 fitted central density
b density gradient (per distance unit)
r_squared R^2 of the log-linear fit
half_distance distance over which density halves (ln 2 / b), when b > 0

Applies when

positive densities
distances not all equal

Example

m.ask("Density gradient: r = [1, 2, 4, 6, 8, 10], D = [9000, 7600, 5300, 3900, 2800, 2000]")

→ D0 = 10546.4, b = 0.166383, r_squared = 0.999624, half_distance = 4.16598

Versions & files

Releases and repository files

Files

15 files, 44.1 MB in total. The largest:

File	Size
model.safetensors	18.9 MB
dist/geospatial_analytics_ai999x-1.0.0-py3-none-any.whl	14.9 MB
onnx/search/encoder.int8.onnx	4.9 MB
docs/geospatial_analytics_ai999x_capabilities.pdf	4.2 MB
onnx/search/index.json	403.3 KB
docs/geospatial_analytics_ai999x_capabilities.html	357 KB
docs/geospatial_analytics_ai999x_capabilities.md	278.7 KB
onnx/search/tokenizer.json	176.4 KB
README.md	28.9 KB
LICENSE	11.4 KB
evaluation/release_gate.md	2.3 KB
evaluation/KNOWLEDGE_CARD.md	1.4 KB

Browse every file ↗

Questions

Frequently asked

What does geospatial-analytics-ai999x do?

Geospatial Analytics answers questions in geospatial analytics: map projections and coordinates, geodesic distance and area, spatial statistics, remote-sensing indices, solar position and terrain, and urban and network analysis. Each answer gives the value with its unit, the steps and the formula in LaTeX.

How do I install geospatial-analytics-ai999x?

With pip: pip install geospatial-analytics-ai999x. The same files are on Hugging Face (hf download AI999X/geospatial-analytics-ai999x --local-dir geospatial-analytics-ai999x) and GitHub.

Does it run offline?

Yes. It is a Python package that runs on an ordinary CPU without a network connection (Python >=3.10); no GPU and no API key are needed.

How many capabilities does geospatial-analytics-ai999x have?

99 capabilities in 14 topics, including Map projections, Spatial statistics, Ellipsoidal geodesy, Terrain analysis, Coordinate systems and conversions.

Which field and tasks does it cover?

Field: Geospatial. Tasks: Question Answering, Token Classification, Text Classification, Sentence Similarity, Feature Extraction.

What does an answer look like?

Asked "Geodesic destination: lat1 = 52.2, lon1 = 0.12, azimuth = 45, distance = 100 km, ellipsoid = WGS84", it answers: lat2 = 52.8308, lon2 = 1.16915, azimuth2 = 45.8325.

What happens when a question cannot be answered?

It says so and names what is missing or out of range, instead of guessing a number.

Explore every model.

More in Geospatial All models Capability orb