brycewang-stanford/gis-remote-sensing-guide

GIS and Remote Sensing Guide

A comprehensive skill for conducting geospatial analysis and remote sensing research. Covers data acquisition from satellite platforms, spatial analysis with open-source tools, and publication-quality map production.

Satellite Data Sources

Key Earth Observation Platforms

| Platform | Provider | Spatial Res. | Revisit | Free? | Use Case | |----------|----------|-------------|---------|-------|----------| | Landsat 8/9 | USGS/NASA | 30m (MS), 15m (pan) | 16 days | Yes | Land cover, NDVI time series | | Sentinel-2 | ESA/Copernicus | 10m | 5 days | Yes | Agriculture, urban mapping | | MODIS | NASA | 250m-1km | 1-2 days | Yes | Large-scale vegetation, fire | | Sentinel-1 | ESA | 5-20m | 6 days | Yes | SAR, flood mapping, deformation | | SRTM/ASTER | NASA | 30m | N/A | Yes | Digital elevation models |

Data Download with Python

import ee

# Initialize Google Earth Engine
ee.Initialize()

def get_sentinel2_composite(aoi: ee.Geometry, start: str, end: str,
                             cloud_max: int = 20) -> ee.Image:
    """
    Create a cloud-free Sentinel-2 composite.

    Args:
        aoi: Area of interest as ee.Geometry
        start: Start date (YYYY-MM-DD)
        end: End date (YYYY-MM-DD)
        cloud_max: Maximum cloud cover percentage
    """
    collection = (ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
                  .filterBounds(aoi)
                  .filterDate(start, end)
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', cloud_max)))

    # Cloud masking using SCL band
    def mask_clouds(image):
        scl = image.select('SCL')
        mask = scl.neq(3).And(scl.neq(8)).And(scl.neq(9)).And(scl.neq(10))
        return image.updateMask(mask)

    return collection.map(mask_clouds).median().clip(aoi)

# Define study area
study_area = ee.Geometry.Rectangle([116.0, 39.5, 117.0, 40.5])  # Beijing region
composite = get_sentinel2_composite(study_area, '2024-06-01', '2024-09-30')

Spatial Analysis with GeoPandas

Vector Data Processing

import geopandas as gpd
from shapely.geometry import Point

def spatial_join_analysis(points_gdf: gpd.GeoDataFrame,
                          polygons_gdf: gpd.GeoDataFrame,
                          agg_col: str) -> gpd.GeoDataFrame:
    """
    Perform spatial join and aggregate point data within polygons.
    """
    joined = gpd.sjoin(points_gdf, polygons_gdf, how='inner', predicate='within')
    summary = joined.groupby('index_right').agg(
        count=(agg_col, 'count'),
        mean_value=(agg_col, 'mean'),
        std_value=(agg_col, 'std')
    ).reset_index()
    result = polygons_gdf.merge(summary, left_index=True, right_on='index_right')
    return result

# Example: aggregate soil samples within administrative boundaries
soil_samples = gpd.read_file('soil_data.geojson')
admin_bounds = gpd.read_file('admin_boundaries.shp')
result = spatial_join_analysis(soil_samples, admin_bounds, 'pH_value')

Remote Sensing Indices

Vegetation and Water Indices

import rasterio
import numpy as np

def compute_indices(image_path: str) -> dict:
    """Compute common remote sensing spectral indices."""
    with rasterio.open(image_path) as src:
        red = src.read(3).astype(float)    # Band 4 in Sentinel-2
        nir = src.read(4).astype(float)    # Band 8
        green = src.read(2).astype(float)  # Band 3
        swir = src.read(5).astype(float)   # Band 11

    # Normalized Difference Vegetation Index
    ndvi = (nir - red) / (nir + red + 1e-10)

    # Normalized Difference Water Index
    ndwi = (green - nir) / (green + nir + 1e-10)

    # Normalized Burn Ratio
    nbr = (nir - swir) / (nir + swir + 1e-10)

    return {'NDVI': ndvi, 'NDWI': ndwi, 'NBR': nbr}

Map Production

For publication-quality maps, always include: scale bar, north arrow, coordinate reference system label, legend, and data source attribution. Use matplotlib with cartopy for projected maps, or folium for interactive web maps. Export at 300 DPI minimum for journal submissions.

Coordinate Reference Systems

Always verify and document the CRS. Use EPSG codes (e.g., EPSG:4326 for WGS84, EPSG:32650 for UTM Zone 50N). Reproject all layers to a common CRS before spatial operations to avoid misalignment errors.