Geo Plotting

{
  "type": "Feature",
  "geometry": {
    "type": "Point",
    "coordinates": [125.6, 10.1]
  },
  "properties": {
    "name": "Dinagat Islands"
  }
}

val usaStates =
    GeoDataFrame.readGeoJson("https://raw.githubusercontent.com/AndreiKingsley/datasets/refs/heads/main/USA.json")

usaStates.df

usaStates.df.geometry.type()

org.locationtech.jts.geom.Geometry

usaStates.df.geometry.map { it::class }.distinct().toList()

[class org.locationtech.jts.geom.Polygon, class org.locationtech.jts.geom.MultiPolygon]

usaStates.crs

null

val worldCities =
    GeoDataFrame.readShapefile("https://github.com/AndreiKingsley/datasets/raw/refs/heads/main/ne_10m_populated_places_simple/ne_10m_populated_places_simple.shp")

worldCities.df

worldCities.df.geometry.type()

org.locationtech.jts.geom.Point

worldCities.crs

GEOGCS["GCS_WGS_1984",
  DATUM["D_WGS_1984",
    SPHEROID["WGS_1984", 6378137.0, 298.257223563]],
  PRIMEM["Greenwich", 0.0],
  UNIT["degree", 0.017453292519943295],
  AXIS["Longitude", EAST],
  AXIS["Latitude", NORTH]]

usaStates.plot {
    // `geoPolygon` uses polygons and multipolygons
    // from the `geometry` column of `usaStates` inner DataFrame
    geoPolygon()
}

usaStates.plot {
    geoPolygon {
        fillColor(name) { legend.type = LegendType.None } // Hide legend
        borderLine {
            width = 0.1
            color = Color.BLACK
        }
    }
}

usaStates.plot {
    // `geoMap` takes polygons from the `geometry`
    // column of `usaStates` inner DataFrame
    geoMap()
    // Add a new dataset using the `worldCities` GeoDataFrame.
    // Layers created within this scope will use it as their base dataset
    // instead of the initial one
    withData(worldCities) {
        // `geoPoints` takes points from the `geometry`
        // column of `worldCities` inner DataFrame
        geoPoints {
            size = 1.5
        }
    }
}

// import mergePolygons utility
import org.jetbrains.kotlinx.kandy.letsplot.geo.util.mergePolygons

// Experimental function that merges a collection of polygons and
// multipolygons into a single multipolygon
val usaPolygon: MultiPolygon = usaStates.df.geometry.mergePolygons()

plot {
    // `geoPolygon` and `geoMap` can accept a single `Polygon` or `MultiPolygon`
    geoMap(usaPolygon)
}

val usaCities = worldCities.modify {
    // Filter the DataFrame to include only points inside the `usaPolygon`
    filter {
        // `usaPolygon.contains(geometry)` checks if the `geometry` (a Point)
        // from the current row of `worldCities` is within the `usaPolygon`
        usaPolygon.contains(geometry)
    }
        // Take 30 most populous cities.
        // Sort the remaining rows by population size in descending order
        .sortByDesc {
            pop_min
        }
        // Select the top 30 rows.
        .take(30)
}

usaStates.plot {
    geoMap()
    withData(usaCities) {
        geoPoints {
            tooltips(title = value(name)) {
                line("population", value(pop_min))
            }
        }
    }
}

val usa48 = usaStates.modify {
    filter {
        name !in listOf("Alaska", "Hawaii", "Puerto Rico")
    }
}

usa48.plot { geoMap() }

val usaAdjusted = usaStates.modify {
    // Custom extensions for `Geometry` based on JTS API.
    // Scale and move Alaska:
    update { geometry }.where { name == "Alaska" }.with {
        it.scaleAroundCenter(0.5).translate(40.0, -40.0)
    }
        // Move Hawaii and Puerto Rico:
        .update { geometry }.where { name == "Hawaii" }.with { it.translate(65.0, 0.0) }
        .update { geometry }.where { name == "Puerto Rico" }.with { it.translate(-10.0, 5.0) }
}

usaAdjusted.plot { geoMap() }

usa48.plot {
    geoMap()
    withData(usa48.modify {
        update { geometry }.with { it.centroid }
    }) {
        geoPoints()
    }
}

val conusAlbersCrs = CRS.decode("EPSG:5070", true)
val usaAlbers = usa48.applyCrs(conusAlbersCrs)
usaAlbers.crs
println("CRS.equalsIgnoreMetadata(usaAlbers.crs, CRS.decode("EPSG:5070", true)) is ${CRS.equalsIgnoreMetadata(usaAlbers.crs, CRS.decode("EPSG:5070", true))}") // true

PROJCS["NAD83 / Conus Albers",
  GEOGCS["NAD83",
    DATUM["North American Datum 1983",
      SPHEROID["GRS 1980", 6378137.0, 298.257222101, AUTHORITY["EPSG","7019"]],
...

usaAlbers.plot {
    // Polygons will work exactly the same -
    // no special coordinates transformation is applied
    // for GeoDF with unsupported CRS
    geoMap()
}

import org.locationtech.jts.geom.*
import kotlin.math.*

fun greatCircleLineString(start: Point, end: Point, n: Int = 100): LineString {
    val factory = GeometryFactory()

    val startLat = Math.toRadians(start.y)
    val startLon = Math.toRadians(start.x)
    val endLat = Math.toRadians(end.y)
    val endLon = Math.toRadians(end.x)

    val deltaLon = endLon - startLon
    val cosStartLat = cos(startLat)
    val cosEndLat = cos(endLat)
    val sinStartLat = sin(startLat)
    val sinEndLat = sin(endLat)
    val a = cosStartLat * cosEndLat * cos(deltaLon) + sinStartLat * sinEndLat
    val angularDistance = acos(a)

    if (angularDistance == 0.0) {
        return factory.createLineString(arrayOf(start.coordinate, end.coordinate))
    }

    val coordinates = mutableListOf<Coordinate>()
    for (i in 0..n) {
        val fraction = i.toDouble() / n
        val sinAngularDistance = sin(angularDistance)
        val A = sin((1 - fraction) * angularDistance) / sinAngularDistance
        val B = sin(fraction * angularDistance) / sinAngularDistance

        val x = A * cosStartLat * cos(startLon) + B * cosEndLat * cos(endLon)
        val y = A * cosStartLat * sin(startLon) + B * cosEndLat * sin(endLon)
        val z = A * sinStartLat + B * sinEndLat

        val lat = atan2(z, sqrt(x * x + y * y))
        val lon = atan2(y, x)

        coordinates.add(Coordinate(Math.toDegrees(lon), Math.toDegrees(lat)))
    }

    return factory.createLineString(coordinates.toTypedArray())
}

fun takeCity(name: String) = usaCities.df.filter { it.name == name }.single().geometry

val newYork = takeCity("New York")
val losAngeles = takeCity("Los Angeles")

val curveNY_LA = greatCircleLineString(newYork, losAngeles)

usa48.plot {
    geoMap { alpha = 0.5 }
    geoPath(curveNY_LA) { width = 1.5 }
    geoPoints(listOf(newYork, losAngeles)) {
        size = 8.0
        color = Color.RED
    }
}

// `.bounds()` function calculates the minimum bounding box
// of all geometries in the `geometry` column of a `GeoDataFrame`,
// returning it as an `Envelope`
val usa48Bounds: Envelope = usa48.bounds().also {
    // Use JTS API for in-place envelope expansion
    it.expandBy(1.0)
}

usa48.plot {
    geoMap()
    geoRectangles(usa48Bounds) {
        alpha = 0.0
        borderLine {
            width = 2.0
            color = Color.GREY
        }
    }
}

usa48.plot {
    geoMap()
    geoRectangles()
}

usaCities.writeGeoJson("usa_cities.geojson")

GeoDataFrame.readGeoJson("usa_cities.geojson").plot { geoPoints() }

// All geometries should be the same type (Shapefile restriction),
// but we have both `Polygon` and `MultiPolygon`.
// Cast them all into MultiPolygons
usa48.modify {
    convert { geometry }.with {
        when (it) {
            // Cast `Polygon` to `MultiPolygon` with a single entity
            is Polygon -> it.toMultiPolygon()
            is MultiPolygon -> it
            else -> error("not a polygonal geometry")
        }
    }
}
    // All files comprising the Shapefile will be saved to
    // a directory named "usa_48" and will have the same base name
    .writeShapefile("usa_48")

GeoDataFrame.readShapefile("usa_48/usa_48.shp").plot { geoMap() }