---

name: geomaster description: Comprehensive geospatial science skill covering remote sensing, GIS, spatial analysis, machine learning for earth observation, and 30+ scientific domains. Supports satellite imagery processing (Sentinel, Landsat, MODIS, SAR, hyperspectral), vector and raster data operations, spatial statistics, point cloud processing, network analysis, and 7 programming languages (Python, R, Julia, JavaScript, C++, Java, Go) with 500+ code examples. Use for remote sensing workflows, GIS analysis, spatial ML, Earth observation data processing, terrain analysis, hydrological modeling, marine spatial analysis, atmospheric science, and any geospatial computation task. license: MIT License metadata: skill-author: K-Dense Inc.

GeoMaster

GeoMaster is a comprehensive geospatial science skill covering the full spectrum of geographic information systems, remote sensing, spatial analysis, and machine learning for Earth observation. This skill provides expert knowledge across 70+ topics with 500+ code examples in 7 programming languages.

Installation

Core Python Geospatial Stack

# Install via conda (recommended for geospatial dependencies)
conda install -c conda-forge gdal rasterio fiona shapely pyproj geopandas

# Or via uv
uv pip install geopandas rasterio fiona shapely pyproj

Remote Sensing & Image Processing

# Core remote sensing libraries
uv pip install rsgislib torchgeo eo-learn

# For Google Earth Engine
uv pip install earthengine-api

# For SNAP integration
# Download from: https://step.esa.int/main/download/

GIS Software Integration

# QGIS Python bindings (usually installed with QGIS)
# ArcPy requires ArcGIS Pro installation

# GRASS GIS
conda install -c conda-forge grassgrass

# SAGA GIS
conda install -c conda-forge saga-gis

Machine Learning for Geospatial

# Deep learning for remote sensing
uv pip install torch-geometric tensorflow-caney

# Spatial machine learning
uv pip install libpysal esda mgwr
uv pip install scikit-learn xgboost lightgbm

Point Cloud & 3D

# LiDAR processing
uv pip install laspy pylas

# Point cloud manipulation
uv pip install open3d pdal

# Photogrammetry
uv pip install opendm

Network & Routing

# Street network analysis
uv pip install osmnx networkx

# Routing engines
uv pip install osrm pyrouting

Visualization

# Static mapping
uv pip install cartopy contextily mapclassify

# Interactive web maps
uv pip install folium ipyleaflet keplergl

# 3D visualization
uv pip install pydeck pythreejs

Big Data & Cloud

# Distributed geospatial processing
uv pip install dask-geopandas

# Xarray for multidimensional arrays
uv pip install xarray rioxarray

# Planetary Computer
uv pip install pystac-client planetary-computer

Database Support

# PostGIS
conda install -c conda-forge postgis

# SpatiaLite
conda install -c conda-forge spatialite

# GeoAlchemy2 for SQLAlchemy
uv pip install geoalchemy2

Additional Programming Languages

# R geospatial packages
# install.packages(c("sf", "terra", "raster", "terra", "stars"))

# Julia geospatial packages
# import Pkg; Pkg.add(["ArchGDAL", "GeoInterface", "GeoStats.jl"])

# JavaScript (Node.js)
# npm install @turf/turf terraformer-arcgis-parser

# Java
# Maven: org.geotools:gt-main

Quick Start

Reading Satellite Imagery and Calculating NDVI

import rasterio
import numpy as np

# Open Sentinel-2 imagery
with rasterio.open('sentinel2.tif') as src:
    # Read red (B04) and NIR (B08) bands
    red = src.read(4)
    nir = src.read(8)

    # Calculate NDVI
    ndvi = (nir.astype(float) - red.astype(float)) / (nir + red)
    ndvi = np.nan_to_num(ndvi, nan=0)

    # Save result
    profile = src.profile
    profile.update(count=1, dtype=rasterio.float32)

    with rasterio.open('ndvi.tif', 'w', **profile) as dst:
        dst.write(ndvi.astype(rasterio.float32), 1)

print(f"NDVI range: {ndvi.min():.3f} to {ndvi.max():.3f}")

Spatial Analysis with GeoPandas

import geopandas as gpd

# Load spatial data
zones = gpd.read_file('zones.geojson')
points = gpd.read_file('points.geojson')

# Ensure same CRS
if zones.crs != points.crs:
    points = points.to_crs(zones.crs)

# Spatial join (points within zones)
joined = gpd.sjoin(points, zones, how='inner', predicate='within')

# Calculate statistics per zone
stats = joined.groupby('zone_id').agg({
    'value': ['count', 'mean', 'std', 'min', 'max']
}).round(2)

print(stats)

Google Earth Engine Time Series

import ee
import pandas as pd

# Initialize Earth Engine
ee.Initialize(project='your-project-id')

# Define region of interest
roi = ee.Geometry.Point([-122.4, 37.7]).buffer(10000)

# Get Sentinel-2 collection
s2 = (ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
      .filterBounds(roi)
      .filterDate('2020-01-01', '2023-12-31')
      .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20)))

# Add NDVI band
def add_ndvi(image):
    ndvi = image.normalizedDifference(['B8', 'B4']).rename('NDVI')
    return image.addBands(ndvi)

s2_ndvi = s2.map(add_ndvi)

# Extract time series
def extract_series(image):
    stats = image.reduceRegion(
        reducer=ee.Reducer.mean(),
        geometry=roi.centroid(),
        scale=10,
        maxPixels=1e9
    )
    return ee.Feature(None, {
        'date': image.date().format('YYYY-MM-dd'),
        'ndvi': stats.get('NDVI')
    })

series = s2_ndvi.map(extract_series).getInfo()
df = pd.DataFrame([f['properties'] for f in series['features']])
df['date'] = pd.to_datetime(df['date'])
print(df.head())

Core Concepts

Coordinate Reference Systems (CRS)

Understanding CRS is fundamental to geospatial work:

• Geographic CRS: EPSG:4326 (WGS 84) - uses lat/lon degrees
• Projected CRS: EPSG:3857 (Web Mercator) - uses meters
• UTM Zones: EPSG:326xx (North), EPSG:327xx (South) - minimizes distortion

See coordinate-systems.md for comprehensive CRS reference.

Vector vs Raster Data

Vector Data: Points, lines, polygons with discrete boundaries

• Shapefiles, GeoJSON, GeoPackage, PostGIS
• Best for: administrative boundaries, roads, infrastructure

Raster Data: Grid of cells with continuous values

• GeoTIFF, NetCDF, HDF5, COG
• Best for: satellite imagery, elevation, climate data

Spatial Data Types

| Type | Examples | Libraries | |------|----------|-----------| | Vector | Shapefiles, GeoJSON, GeoPackage | GeoPandas, Fiona, GDAL | | Raster | GeoTIFF, NetCDF, IMG | Rasterio, GDAL, Xarray | | Point Cloud | LAZ, LAS, PCD | Laspy, PDAL, Open3D | | Topology | TopoJSON, TopoArchive | TopoJSON, NetworkX | | Spatiotemporal | Trajectories, Time-series | MovingPandas, PyTorch Geometric |

OGC Standards

Key Open Geospatial Consortium standards:

• WMS: Web Map Service - raster maps
• WFS: Web Feature Service - vector data
• WCS: Web Coverage Service - raster coverage
• WPS: Web Processing Service - geoprocessing
• WMTS: Web Map Tile Service - tiled maps

Common Operations

Remote Sensing Operations

Spectral Indices Calculation

import rasterio
import numpy as np

def calculate_indices(image_path, output_path):
    """Calculate NDVI, EVI, SAVI, and NDWI from Sentinel-2."""
    with rasterio.open(image_path) as src:
        # Read bands: B2=Blue, B3=Green, B4=Red, B8=NIR, B11=SWIR1
        blue = src.read(2).astype(float)
        green = src.read(3).astype(float)
        red = src.read(4).astype(float)
        nir = src.read(8).astype(float)
        swir1 = src.read(11).astype(float)

        # Calculate indices
        ndvi = (nir - red) / (nir + red + 1e-8)
        evi = 2.5 * (nir - red) / (nir + 6*red - 7.5*blue + 1)
        savi = ((nir - red) / (nir + red + 0.5)) * 1.5
        ndwi = (green - nir) / (green + nir + 1e-8)

        # Stack and save
        indices = np.stack([ndvi, evi, savi, ndwi])
        profile = src.profile
        profile.update(count=4, dtype=rasterio.float32)

        with rasterio.open(output_path, 'w', **profile) as dst:
            dst.write(indices)

# Usage
calculate_indices('sentinel2.tif', 'indices.tif')

Image Classification

from sklearn.ensemble import RandomForestClassifier
import geopandas as gpd
import rasterio
from rasterio.features import rasterize
import numpy as np

def classify_imagery(raster_path, training_gdf, output_path):
    """Train Random Forest classifier and classify imagery."""
    # Load imagery
    with rasterio.open(raster_path) as src:
        image = src.read()
        profile = src.profile
        transform = src.transform

    # Extract training data
    X_train, y_train = [], []

    for _, row in training_gdf.iterrows():
        mask = rasterize(
            [(row.geometry, 1)],
            out_shape=(profile['height'], profile['width']),
            transform=transform,
            fill=0,
            dtype=np.uint8
        )
        pixels = image[:, mask > 0].T
        X_train.extend(pixels)
        y_train.extend([row['class_id']] * len(pixels))

    X_train = np.array(X_train)
    y_train = np.array(y_train)

    # Train classifier
    rf = RandomForestClassifier(n_estimators=100, max_depth=20, n_jobs=-1)
    rf.fit(X_train, y_train)

    # Predict full image
    image_reshaped = image.reshape(image.shape[0], -1).T
    prediction = rf.predict(image_reshaped)
    prediction = prediction.reshape(profile['height'], profile['width'])

    # Save result
    profile.update(dtype=rasterio.uint8, count=1)
    with rasterio.open(output_path, 'w', **profile) as dst:
        dst.write(prediction.astype(rasterio.uint8), 1)

    return rf

Vector Operations

import geopandas as gpd
from shapely.ops import unary_union

# Buffer analysis
gdf['buffer_1km'] = gdf.geometry.to_crs(epsg=32633).buffer(1000)

# Spatial relationships
intersects = gdf[gdf.geometry.intersects(other_geometry)]
contains = gdf[gdf.geometry.contains(point_geometry)]

# Geometric operations
gdf['centroid'] = gdf.geometry.centroid
gdf['convex_hull'] = gdf.geometry.convex_hull
gdf['simplified'] = gdf.geometry.simplify(tolerance=0.001)

# Overlay operations
intersection = gpd.overlay(gdf1, gdf2, how='intersection')
union = gpd.overlay(gdf1, gdf2, how='union')
difference = gpd.overlay(gdf1, gdf2, how='difference')

Terrain Analysis

import rasterio
from rasterio.features import shapes
import numpy as np

def calculate_terrain_metrics(dem_path):
    """Calculate slope, aspect, hillshade from DEM."""
    with rasterio.open(dem_path) as src:
        dem = src.read(1)
        transform = src.transform

    # Calculate gradients
    dy, dx = np.gradient(dem)

    # Slope (in degrees)
    slope = np.arctan(np.sqrt(dx**2 + dy**2)) * 180 / np.pi

    # Aspect (in degrees, clockwise from north)
    aspect = np.arctan2(-dy, dx) * 180 / np.pi
    aspect = (90 - aspect) % 360

    # Hillshade
    azimuth = 315
    altitude = 45
    azimuth_rad = np.radians(azimuth)
    altitude_rad = np.radians(altitude)

    hillshade = (np.sin(altitude_rad) * np.sin(np.radians(slope)) +
                 np.cos(altitude_rad) * np.cos(np.radians(slope)) *
                 np.cos(np.radians(aspect) - azimuth_rad))

    return slope, aspect, hillshade

Network Analysis

import osmnx as ox
import networkx as nx

# Download street network
G = ox.graph_from_place('San Francisco, CA', network_type='drive')

# Add speeds and travel times
G = ox.add_edge_speeds(G)
G = ox.add_edge_travel_times(G)

# Find shortest path
orig_node = ox.distance.nearest_nodes(G, -122.4, 37.7)
dest_node = ox.distance.nearest_nodes(G, -122.3, 37.8)
route = nx.shortest_path(G, orig_node, dest_node, weight='travel_time')

# Calculate accessibility
accessibility = {}
for node in G.nodes():
    subgraph = nx.ego_graph(G, node, radius=5, distance='time')
    accessibility[node]