!pip install pystac-client planetary-computer odc-stac geopandas pandas

"""
Search Sentinel-2 imagery from Planetary Computer
and return JSON.
"""

import json
from pystac_client import Client

STAC_URL = "https://planetarycomputer.microsoft.com/api/stac/v1"


def search_sentinel2(
    geojson_fc: dict,
    start_date: str,
    end_date: str,
    max_cloud_cover: int = 20,
):
    # --------------------------------
    # convert FeatureCollection → geometry
    # --------------------------------
    if geojson_fc["type"] == "FeatureCollection":
        geometry = geojson_fc["features"][0]["geometry"]

    elif geojson_fc["type"] == "Feature":
        geometry = geojson_fc["geometry"]

    else:
        geometry = geojson_fc

    # --------------------------------
    # connect stac
    # --------------------------------
    catalog = Client.open(STAC_URL)

    search = catalog.search(
        collections=["sentinel-2-l2a"],
        intersects=geometry,
        datetime=f"{start_date}/{end_date}",
        query={
            "eo:cloud_cover": {
                "lte": max_cloud_cover
            }
        },
    )

    items = list(search.items())

    # --------------------------------
    # no results
    # --------------------------------
    if not items:
        return {
            "count": 0,
            "results": []
        }

    # --------------------------------
    # results
    # --------------------------------
    results = []

    for item in items:
        results.append(
            {
                "id": item.id,
                "date": item.datetime.date().isoformat(),
                "datetime": item.datetime.isoformat(),
                "cloud_cover": item.properties.get(
                    "eo:cloud_cover"
                ),
                "product": item.properties.get(
                    "s2:product_uri"
                ),
                "platform": item.properties.get(
                    "platform"
                ),
            }
        )

    # sort by date
    results = sorted(
        results,
        key=lambda x: x["date"]
    )

    return {
        "count": len(results),
        "results": results,
    }


# --------------------------------
# Example
# --------------------------------

geojson_fc = {
    "type": "FeatureCollection",
    "features": [
        {
            "type": "Feature",
            "properties": {},
            "geometry": {
                "type": "Polygon",
                "coordinates": [[
                    [76.8149, 17.6984],
                    [76.8131, 17.6972],
                    [76.8155, 17.6961],
                    [76.8172, 17.6976],
                    [76.8149, 17.6984],
                ]],
            },
        }
    ],
}

data = search_sentinel2(
    geojson_fc=geojson_fc,
    start_date="2025-12-01",
    end_date="2025-12-31",
    max_cloud_cover=15,
)

print(json.dumps(data, indent=2))

!pip install pystac-client planetary-computer odc-stac geopandas pandas rasterio rioxarray

import os
import rasterio
from shapely.geometry import shape
from pystac_client import Client
import planetary_computer as pc
from odc.stac import load
import rioxarray


STAC_URL = "https://planetarycomputer.microsoft.com/api/stac/v1"


def download_clipped_raw_bands(
    geojson_data: dict,
    item_id: str,
    bands: list,
    output_dir: str,
):
    """
    Download clipped raw Sentinel-2 bands.

    Parameters
    ----------
    geojson_data : dict
        GeoJSON Polygon / Feature / FeatureCollection

    item_id : str
        Sentinel-2 item id

    bands : list
        Example:
        ["B02", "B03", "B04", "B08"]

    output_dir : str
        Folder to save TIFFs
    """

    # --------------------------
    # geometry
    # --------------------------
    if geojson_data["type"] == "FeatureCollection":
        geom_geojson = geojson_data["features"][0]["geometry"]

    elif geojson_data["type"] == "Feature":
        geom_geojson = geojson_data["geometry"]

    else:
        geom_geojson = geojson_data

    geom = shape(geom_geojson)

    # --------------------------
    # stac
    # --------------------------
    client = Client.open(
        STAC_URL,
        modifier=pc.sign_inplace,
    )

    search = client.search(
        collections=["sentinel-2-l2a"],
        ids=[item_id],
    )

    items = list(search.items())

    if not items:
        raise ValueError(
            f"No item found for {item_id}"
        )

    signed_items = [
        pc.sign(item)
        for item in items
    ]

    # --------------------------
    # load only required bands
    # --------------------------
    ds = load(
        items=signed_items,
        geopolygon=geom,
        groupby="solar_day",
        bands=bands,
    )

    # exact polygon clip
    ds = ds.rio.clip(
        [geom],
        crs="EPSG:4326",
    )

    # first image only
    ds = ds.isel(time=0)

    # create folder
    os.makedirs(
        output_dir,
        exist_ok=True,
    )

    saved_files = []

    # --------------------------
    # save each band
    # --------------------------
    for band in bands:

        output_path = os.path.join(
            output_dir,
            f"{band}.tif",
        )

        ds[band].rio.to_raster(
            output_path
        )

        saved_files.append(
            output_path
        )

    return {
        "status": "success",
        "item_id": item_id,
        "bands": bands,
        "files": saved_files,
    }


# -----------------------------------
# Example
# -----------------------------------

res = download_clipped_raw_bands(
    geojson_data=geojson_fc,
    item_id="S2C_MSIL2A_20251214T053241_R105_T43QDF_20251214T090710",
    bands=[
        "B02",
        "B03",
        "B04",
        "B08",
    ],
    output_dir="/content/raw_bands",
)

print(res)

{
    "status": "success",
    "item_id": "S2C_MSIL2A_20251214T053241_R105_T43QDF_20251214T090710",
    "bands": [
        "B02",
        "B03",
        "B04",
        "B08"
    ],
    "files": [
        "/content/raw_bands/B02.tif",
        "/content/raw_bands/B03.tif",
        "/content/raw_bands/B04.tif",
        "/content/raw_bands/B08.tif"
    ]
}

!pip install pystac-client planetary-computer odc-stac geopandas pandas rasterio rioxarray

import numpy as np
import rasterio
from shapely.geometry import shape
from pystac_client import Client
import planetary_computer as pc
from odc.stac import load
import rioxarray


STAC_URL = "https://planetarycomputer.microsoft.com/api/stac/v1"


def download_clipped_index_tiff(
    geojson_data: dict,
    item_id: str,
    index: str,
    output_path: str,
):
    """
    Download clipped Sentinel-2 NDVI/EVI GeoTIFF.

    Parameters
    ----------
    geojson_data : dict
    item_id : str
    index : str ("ndvi" or "evi")
    output_path : str
    """

    # ---------------------------------
    # geometry
    # ---------------------------------
    if geojson_data["type"] == "FeatureCollection":
        geom_geojson = geojson_data["features"][0]["geometry"]
    elif geojson_data["type"] == "Feature":
        geom_geojson = geojson_data["geometry"]
    else:
        geom_geojson = geojson_data

    geom = shape(geom_geojson)

    # ---------------------------------
    # stac search by id
    # ---------------------------------
    client = Client.open(
        STAC_URL,
        modifier=pc.sign_inplace,
    )

    search = client.search(
        collections=["sentinel-2-l2a"],
        ids=[item_id],
    )

    items = list(search.items())

    if not items:
        raise ValueError(
            f"No item found for {item_id}"
        )

    signed_items = [
        pc.sign(item)
        for item in items
    ]

    # ---------------------------------
    # required bands
    # ---------------------------------
    bands = ["B04", "B08"]

    if index.lower() == "evi":
        bands.append("B02")

    # ---------------------------------
    # load clipped to polygon
    # ---------------------------------
    ds = load(
        items=signed_items,
        geopolygon=geom,
        groupby="solar_day",
        bands=bands,
    )

    # exact polygon clip
    ds = ds.rio.clip(
        [geom],
        crs="EPSG:4326",
    )

    # ---------------------------------
    # calculate index
    # ---------------------------------
    red = ds["B04"]
    nir = ds["B08"]

    if index.lower() == "ndvi":
        result = (
            (nir - red)
            / (nir + red + 1e-6)
        )

    elif index.lower() == "evi":
        blue = ds["B02"]

        result = (
            2.5
            * (nir - red)
            / (
                nir
                + 6 * red
                - 7.5 * blue
                + 1
            )
        )

    else:
        raise ValueError(
            "Supported: ndvi or evi"
        )

    # first image only
    result = result.isel(time=0)

    # ---------------------------------
    # save
    # ---------------------------------
    result.rio.to_raster(output_path)

    return {
        "status": "success",
        "item_id": item_id,
        "index": index,
        "output": output_path,
    }


# -------------------------------------
# Example
# -------------------------------------

res = download_clipped_index_tiff(
    geojson_data=geojson_fc,
    item_id="S2C_MSIL2A_20251204T053221_R105_T43QDF_20251204T091915",
    index="ndvi",
    output_path="/content/ndvi_clip.tif",
)

print(res)

{
    "status": "success",
    "item_id": "S2C_MSIL2A_20251204T053221_R105_T43QDF_20251204T091915",
    "index": "ndvi",
    "output": "/content/ndvi_clip.tif"
}

!pip install pystac-client planetary-computer odc-stac geopandas pandas rasterio rioxarray scipy

import numpy as np
import rasterio
import geopandas as gpd
from shapely.geometry import shape
from rasterio.mask import mask as rio_mask
from pystac_client import Client
import planetary_computer as pc
from odc.stac import load
from scipy.ndimage import median_filter


STAC_URL = "https://planetarycomputer.microsoft.com/api/stac/v1"


def download_lulc_tiff(
    geojson_data: dict,
    date_range: str,
    cloud_cover: int,
    output_path: str,
    max_items: int = 10,
    resolution: int = 20,
):
    """
    Download clipped Sentinel-2 LULC GeoTIFF.

    Classes
    -------
    0 -> Water
    1 -> Built-up
    2 -> Barren/Open soil
    3 -> Sparse vegetation
    4 -> Dense vegetation
    255 -> NoData
    """

    # --------------------------------
    # geometry
    # --------------------------------
    if geojson_data["type"] == "FeatureCollection":
        geom_geojson = geojson_data["features"][0]["geometry"]

    elif geojson_data["type"] == "Feature":
        geom_geojson = geojson_data["geometry"]

    else:
        geom_geojson = geojson_data

    geom = shape(geom_geojson)

    gdf = gpd.GeoDataFrame(
        geometry=[geom],
        crs="EPSG:4326",
    )

    bbox = tuple(gdf.total_bounds)

    utm_crs = gdf.estimate_utm_crs()

    # --------------------------------
    # search stac
    # --------------------------------
    client = Client.open(STAC_URL)

    search = client.search(
        collections=["sentinel-2-l2a"],
        intersects=geom.__geo_interface__,
        datetime=date_range,
        query={
            "eo:cloud_cover": {
                "lt": cloud_cover
            }
        },
    )

    items = sorted(
        search.items(),
        key=lambda x: x.properties.get(
            "eo:cloud_cover",
            100,
        ),
    )[:max_items]

    if not items:
        raise ValueError(
            "No scenes found"
        )

    signed_items = [
        pc.sign(item)
        for item in items
    ]

    # --------------------------------
    # load bands
    # --------------------------------
    ds = load(
        items=signed_items,
        bands=[
            "B02",
            "B03",
            "B04",
            "B08",
            "B11",
            "SCL",
        ],
        crs=utm_crs,
        resolution=resolution,
        bbox=bbox,
        groupby="solar_day",
    )

    # --------------------------------
    # cloud mask
    # --------------------------------
    scl = ds["SCL"]

    cloud_mask = scl.isin(
        [1, 3, 8, 9, 10, 11]
    )

    blue = ds["B02"].where(
        ~cloud_mask
    ).astype("float32")

    green = ds["B03"].where(
        ~cloud_mask
    ).astype("float32")

    red = ds["B04"].where(
        ~cloud_mask
    ).astype("float32")

    nir = ds["B08"].where(
        ~cloud_mask
    ).astype("float32")

    swir = ds["B11"].where(
        ~cloud_mask
    ).astype("float32")

    # median composite
    blue = blue.median(dim="time")
    green = green.median(dim="time")
    red = red.median(dim="time")
    nir = nir.median(dim="time")
    swir = swir.median(dim="time")

    # --------------------------------
    # indices
    # --------------------------------
    def safe_index(a, b):
        denom = a + b

        return np.where(
            denom == 0,
            np.nan,
            (a - b) / denom,
        )

    b = blue.values
    g = green.values
    r = red.values
    n = nir.values
    s = swir.values

    ndvi = safe_index(n, r)

    ndbi = safe_index(s, n)

    mndwi = safe_index(g, s)

    # --------------------------------
    # classification
    # --------------------------------
    lulc = np.full(
        ndvi.shape,
        255,
        dtype=np.uint8,
    )

    valid = (
        np.isfinite(ndvi)
        & np.isfinite(ndbi)
        & np.isfinite(mndwi)
    )

    water = (
        valid
        & (mndwi > 0.1)
        & (mndwi > ndvi)
    )
    lulc[water] = 0

    dense_veg = (
        valid
        & ~water
        & (ndvi > 0.5)
    )
    lulc[dense_veg] = 4

    sparse_veg = (
        valid
        & ~water
        & ~dense_veg
        & (ndvi >= 0.2)
        & (ndvi <= 0.5)
    )
    lulc[sparse_veg] = 3

    builtup = (
        valid
        & ~water
        & ~dense_veg
        & ~sparse_veg
        & (ndbi > 0)
        & (ndvi < 0.2)
    )

    lulc[builtup] = 1

    barren = (
        valid
        & (lulc == 255)
    )

    lulc[barren] = 2

    # --------------------------------
    # smoothing
    # --------------------------------
    filtered = median_filter(
        lulc,
        size=3,
    )

    filtered[
        lulc == 255
    ] = 255

    lulc = filtered

    # --------------------------------
    # save temp
    # --------------------------------
    geobox = ds.odc.geobox

    temp_path = (
        output_path
        + ".tmp.tif"
    )

    with rasterio.open(
        temp_path,
        "w",
        driver="GTiff",
        height=lulc.shape[0],
        width=lulc.shape[1],
        count=1,
        dtype=rasterio.uint8,
        crs=geobox.crs,
        transform=geobox.transform,
        nodata=255,
    ) as dst:
        dst.write(
            lulc,
            1,
        )

    # --------------------------------
    # exact clip
    # --------------------------------
    gdf_utm = gdf.to_crs(
        utm_crs
    )

    geom_utm = (
        gdf_utm.geometry.values[0]
    )

    with rasterio.open(
        temp_path
    ) as src:

        clipped, transform = rio_mask(
            src,
            [geom_utm.__geo_interface__],
            crop=True,
            nodata=255,
        )

        meta = src.meta.copy()

    meta.update(
        {
            "height": clipped.shape[1],
            "width": clipped.shape[2],
            "transform": transform,
            "nodata": 255,
        }
    )

    with rasterio.open(
        output_path,
        "w",
        **meta,
    ) as dst:
        dst.write(
            clipped[0],
            1,
        )

    return {
        "status": "success",
        "date_range": date_range,
        "cloud_cover": cloud_cover,
        "output": output_path,
        "classes": {
            0: "Water",
            1: "Built-up",
            2: "Barren",
            3: "Sparse vegetation",
            4: "Dense vegetation",
            255: "NoData",
        },
    }

geojson_fc = {
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "properties": {},
      "geometry": {
        "type": "Polygon",
        "coordinates": [[
          [73.658151, 20.034861],
          [73.658151, 19.977499],
          [73.771634, 19.977499],
          [73.771634, 20.034861],
          [73.658151, 20.034861]
        ]]
      }
    }
  ]
}

res = download_lulc_tiff(
    geojson_data=geojson_fc,
    date_range="2025-01-01/2025-03-31",
    cloud_cover=20,
    output_path="/content/lulc_2025.tif",
)

print(res)

import rasterio

path = "assets/tiff/Tiff_1.tif"  # bundled practice GeoTIFF
with rasterio.open(path) as src:
    band1 = src.read(1)          # 2D NumPy array, first band (height × width)
    profile = src.profile        # driver, dtype, nodata, transform, crs, size
    crs_epsg = src.crs.to_epsg() if src.crs else None
    crs_label = f"EPSG:{crs_epsg}" if crs_epsg else str(src.crs)
    print(src.shape, crs_label, src.dtypes[0])

import numpy as np
import rasterio

with rasterio.open("assets/tiff/Tiff_1.tif") as src:
    arr = src.read(1).astype("float64")
    nodata = src.nodata

    if nodata is not None:
        arr = np.where(arr == nodata, np.nan, arr)

# Basic stats
print("Min:", np.nanmin(arr))
print("Max:", np.nanmax(arr))
print("Mean:", np.nanmean(arr))
print("Median:", np.nanmedian(arr))
print("Std Dev:", np.nanstd(arr))
print("Variance:", np.nanvar(arr))

# Percentiles
print("25th Percentile:", np.nanpercentile(arr, 25))
print("50th Percentile (Median):", np.nanpercentile(arr, 50))
print("75th Percentile:", np.nanpercentile(arr, 75))

# Count info
print("Total Pixels:", arr.size)
print("Valid Pixels:", np.count_nonzero(~np.isnan(arr)))
print("NoData Pixels:", np.count_nonzero(np.isnan(arr)))

import numpy as np
import rasterio
from rasterio.enums import Resampling

with rasterio.open("/content/Tiff_1.tif") as src:
    arr = src.read(1).astype("float32")
    profile = src.profile.copy()
    mask = (arr >= 0.1) & (arr <= 0.5)   # keep elevations 100–500 (example)
    out = np.where(mask, arr, src.nodata).astype(profile["dtype"])

profile.update(dtype=out.dtype, count=1, compress="deflate")
with rasterio.open("/content/Tiff_1_filter.tif", "w", **profile) as dst:
    dst.write(out, 1)

import rasterio
from rasterio.mask import mask

# GeoJSON Polygon: outer ring must close (first point = last point)
clip_poly = {
    "type": "Polygon",
    "coordinates": [
        [
            [
              73.71042841436719,
              20.09023103896574
            ],
            [
              73.69665791730623,
              20.074091366211377
            ],
            [
              73.71214563805651,
              20.063933096227032
            ],
            [
              73.7259168550645,
              20.066013154915694
            ],
            [
              73.72492563666324,
              20.077312229882054
            ],
            [
              73.72123577059406,
              20.087440203990184
            ],
            [
              73.71042841436719,
              20.09023103896574
            ]
        ]
    ],
}
geom = [clip_poly]

with rasterio.open("/content/Tiff_1.tif") as src:
    clipped, transform = mask(src, geom, crop=True, nodata=src.nodata)
    meta = src.meta.copy()
    meta.update({"height": clipped.shape[1], "width": clipped.shape[2], "transform": transform})

with rasterio.open("/content/Tiff_1_clipped.tif", "w", **meta) as dst:
    dst.write(clipped)

import rasterio
from rasterio.merge import merge

src_files = ["/content/Tiff_1.tif", "/content/Tiff_2.tif"]
srcs = [rasterio.open(p) for p in src_files]
mosaic, out_transform = merge(srcs)

out_meta = srcs[0].meta.copy()
out_meta.update({"height": mosaic.shape[1], "width": mosaic.shape[2], "transform": out_transform})
for s in srcs:
    s.close()

with rasterio.open("/content/Tiff_merge.tif", "w", **out_meta) as dst:
    dst.write(mosaic)

import rasterio
from rasterio.enums import Resampling
from rasterio.transform import Affine

src_path = "/content/Tiff_1.tif"
dst_path = "/content/Tiff_1_resampled.tif"
factor = 2  # >1 = downsample, <1 = upsample

with rasterio.open(src_path) as src:
    # New dimensions
    new_height = int(src.height / factor)
    new_width = int(src.width / factor)

    # Read & resample
    data = src.read(
        out_shape=(src.count, new_height, new_width),
        resampling=Resampling.average
    )

    # Update transform
    transform = src.transform * Affine.scale(
        src.width / new_width,
        src.height / new_height
    )

    # Update metadata
    profile = src.profile
    profile.update({
        "height": new_height,
        "width": new_width,
        "transform": transform
    })

# Save output
with rasterio.open(dst_path, "w", **profile) as dst:
    dst.write(data)

print(f"Resampled: {src.height}x{src.width} → {new_height}x{new_width}")

import rasterio
import rasterio.plot
import rasterio.mask
import numpy as np
import matplotlib.pyplot as plt
import geopandas as gpd
from rasterio.warp import calculate_default_transform, reproject, Resampling
from rasterio.windows import from_bounds
import warnings
warnings.filterwarnings('ignore')

# Set up plotting
plt.style.use('default')

# Data directory path - change this if your TIFF files are in a different folder
DATA_DIR = 'data/'

# Load available raster data
dem_path = DATA_DIR + 'dem.tif'
b04_path = DATA_DIR + 'B04.tif'  # Red band
b08_path = DATA_DIR + 'B08.tif'  # NIR band
true_color_path = DATA_DIR + 'true_color.tiff'

print("=== LOADING RASTER DATA ===")
print(f"DEM: {dem_path}")
print(f"Red Band (B04): {b04_path}")
print(f"NIR Band (B08): {b08_path}")
print(f"True Color: {true_color_path}")

# Load the DEM (Digital Elevation Model)
with rasterio.open(dem_path) as dem_src:
    elevation_data = dem_src.read(1)  # Read first band
    dem_transform = dem_src.transform
    dem_crs = dem_src.crs
    dem_bounds = dem_src.bounds
    dem_nodata = dem_src.nodata

print("=== DEM DATA PROPERTIES ===")
print(f"Shape: {elevation_data.shape}")
print(f"Data type: {elevation_data.dtype}")
print(f"CRS: {dem_crs}")
print(f"Bounds: {dem_bounds}")
print(f"NoData value: {dem_nodata}")

# Handle NoData values
if dem_nodata is not None:
    elevation_masked = np.ma.masked_equal(elevation_data, dem_nodata)
    print(f"Valid pixels: {elevation_masked.count():,}")
    print(f"NoData pixels: {elevation_masked.mask.sum():,}")
    # Use masked array for statistics
    print(f"Min elevation: {elevation_masked.min():.1f}m")
    print(f"Max elevation: {elevation_masked.max():.1f}m")
    print(f"Mean elevation: {elevation_masked.mean():.1f}m")
else:
    print(f"Min elevation: {elevation_data.min():.1f}m")
    print(f"Max elevation: {elevation_data.max():.1f}m")
    print(f"Mean elevation: {elevation_data.mean():.1f}m")

# Use the actual transform from the DEM
transform = dem_transform
bounds = dem_bounds
height, width = elevation_data.shape

print("=== RASTER SPATIAL PROPERTIES ===")
print(f"Bounds: {bounds}")
print(f"Transform: {transform}")
print(f"Pixel size X: {transform[0]:.4f}")
print(f"Pixel size Y: {abs(transform[4]):.4f}")
print(f"Width: {width} pixels")
print(f"Height: {height} pixels")

# Calculate pixel coordinates
def pixel_to_coord(row, col, transform):
    """Convert pixel coordinates to geographic coordinates"""
    x = transform[0] * col + transform[1] * row + transform[2]
    y = transform[3] * col + transform[4] * row + transform[5]
    return x, y

# Example: center pixel coordinates
center_row, center_col = height // 2, width // 2
center_x, center_y = pixel_to_coord(center_row, center_col, transform)
print(f"Center pixel ({center_row}, {center_col}) = ({center_x:.2f}, {center_y:.2f})")

# Calculate extent for plotting
extent = [bounds.left, bounds.right, bounds.bottom, bounds.top]
print(f"Plot extent: {extent}")