Merge pull request #4 from gregory-halverson/main

aligning `Modland.jl` package with former `MODLAND` sub-module of `STARSDataFusion.jl`

Author: gregory-halverson

Project.toml

@@ -1,7 +1,7 @@
 name = "Modland"
 uuid = "2a2d4310-a4da-4416-b9e4-81e65019a994"
 authors = ["Gregory Halverson <[email protected]>"]
-version = "0.1.1"
+version = "1.0.0"
 
 [deps]
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"

src/Modland.jl

@@ -1,116 +1,167 @@
 module Modland
 
+# MODLAND.jl: Utilities for working with MODIS land grid and projections
+# Provides functions to convert between latitude/longitude and MODIS sinusoidal projection,
+# and to determine MODIS tile indices for given coordinates or polygons.
+
+
+# Import required packages for geospatial operations and data handling
 import ArchGDAL as AG
 import ArchGDAL
 import GeoDataFrames as GDF
 import GeoFormatTypes as GFT
 using DataFrames
 using Rasters
-using DimensionalData.Dimensions.LookupArrays
 import JSON
 
-# boundaries of sinusodial projection
-UPPER_LEFT_X_METERS = -20015109.355798
-UPPER_LEFT_Y_METERS = 10007554.677899
-LOWER_RIGHT_X_METERS = 20015109.355798
-LOWER_RIGHT_Y_METERS = -10007554.677899
 
-# size across (width or height) of any equal-area sinusoidal target
+# Constants defining the boundaries of the MODIS sinusoidal projection (in meters)
+UPPER_LEFT_X_METERS = -20015109.355798   # Westernmost x coordinate
+UPPER_LEFT_Y_METERS = 10007554.677899    # Northernmost y coordinate
+LOWER_RIGHT_X_METERS = 20015109.355798   # Easternmost x coordinate
+LOWER_RIGHT_Y_METERS = -10007554.677899  # Southernmost y coordinate
+
+# Size (width or height) of a MODIS tile in the sinusoidal projection (in meters)
 TILE_SIZE_METERS = 1111950.5197665554
 
-# boundaries of MODIS land grid
+# Number of rows and columns in the MODIS land grid
 TOTAL_ROWS = 18
 TOTAL_COLUMNS = 36
 
+# Coordinate reference systems (CRS) for WGS84 and MODIS sinusoidal projection
 WGS84 = ProjString("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
 SINUSOIDAL = ProjString("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs")
 
-export UPPER_LEFT_X_METERS, UPPER_LEFT_Y_METERS, LOWER_RIGHT_X_METERS, LOWER_RIGHT_Y_METERS, TILE_SIZE_METERS, TOTAL_ROWS, TOTAL_COLUMNS, WGS84, SINUSOIDAL
-
 
 """
-    latlon_to_sinusoidal(lat::Float64, lon::Float64)::Tuple{Float64,Float64}
+    latlon_to_sinusoidal(lat::Float64, lon::Float64) -> (x::Float64, y::Float64)
 
-Convert latitude and longitude coordinates to sinusoidal projection coordinates.
+Convert latitude and longitude (in degrees, WGS84) to MODIS sinusoidal projection coordinates (in meters).
 
-# Arguments
-- `lat::Float64`: Latitude in degrees.
-- `lon::Float64`: Longitude in degrees.
+Arguments:
+    lat::Float64: Latitude in degrees (-90 to 90)
+    lon::Float64: Longitude in degrees (-180 to 180)
 
-# Returns
-- `Tuple{Float64, Float64}`: Sinusoidal projection coordinates (x, y).
+Returns:
+    Tuple (x, y): Sinusoidal projection coordinates in meters
+
+Throws:
+    Error if latitude or longitude is out of bounds.
 """
 function latlon_to_sinusoidal(lat::Float64, lon::Float64)::Tuple{Float64,Float64}
-    # Check if latitude is within valid bounds
+    # Validate latitude
     if lat < -90 || lat > 90
         error("latitude ($(lat)) out of bounds")
     end
-
-    # Check if longitude is within valid bounds
+    # Validate longitude
     if lon < -180 || lon > 180
         error("longitude ($(lon))) out of bounds")
     end
 
-    # Create a point in latitude and longitude coordinates
+    # Create a point in WGS84 and reproject to sinusoidal
     point_latlon = AG.createpoint(lon, lat)
-    # Reproject the point from WGS84 to sinusoidal projection
     point_sinusoidal = AG.reproject(point_latlon, WGS84, SINUSOIDAL)
-    # Extract the x coordinate from the reprojected point
     x = AG.getx(point_sinusoidal, 0)
-    # Extract the y coordinate from the reprojected point
     y = AG.gety(point_sinusoidal, 0)
 
-    # Return the sinusoidal coordinates as a tuple
     return x, y
 end
 
 export latlon_to_sinusoidal
 
+
 """
-    sinusoidal_to_modland(x::Float64, y::Float64)::String
+    sinusoidal_to_MODLAND(x::Float64, y::Float64) -> String
+
+Given sinusoidal projection coordinates (x, y), return the MODIS land tile index (e.g., "h12v04")
+that contains the point.
 
-Convert sinusoidal projection coordinates to MODIS land tile indices.
+Arguments:
+    x::Float64: Sinusoidal x coordinate (meters)
+    y::Float64: Sinusoidal y coordinate (meters)
 
-# Arguments
-- `x::Float64`: Sinusoidal x coordinate.
-- `y::Float64`: Sinusoidal y coordinate.
+Returns:
+    String: MODIS tile index in the format "hXXvYY"
 
-# Returns
-- `String`: MODIS land tile index in the format "hXXvYY".
+Throws:
+    Error if x or y is out of bounds for the MODIS grid.
 """
-function sinusoidal_to_modland(x::Float64, y::Float64)::String
-    # Check if the x coordinate is within the valid bounds of the sinusoidal projection
+function sinusoidal_to_MODLAND(x::Float64, y::Float64)::String
+    # Check if x is within the valid range
     if x < UPPER_LEFT_X_METERS || x > LOWER_RIGHT_X_METERS
         error("sinusoidal x coordinate ($(x)) out of bounds")
     end
-
-    # Check if the y coordinate is within the valid bounds of the sinusoidal projection
+    # Check if y is within the valid range
     if y < LOWER_RIGHT_Y_METERS || y > UPPER_LEFT_Y_METERS
         error("sinusoidal y ($(y))) coordinate out of bounds")
     end
 
-    # Calculate the horizontal tile index
+    # Calculate horizontal (column) and vertical (row) indices
     horizontal_index = Int(floor((x - UPPER_LEFT_X_METERS) / TILE_SIZE_METERS))
-    # Calculate the vertical tile index
     vertical_index = Int(floor((-1 * (y + LOWER_RIGHT_Y_METERS)) / TILE_SIZE_METERS))
 
-    # Adjust the horizontal index if it is at the boundary
+    # Clamp indices to valid range if on the edge
     if horizontal_index == TOTAL_COLUMNS
         horizontal_index -= 1
     end
-
-    # Adjust the vertical index if it is at the boundary
     if vertical_index == TOTAL_ROWS
         vertical_index -= 1
     end
 
-    # Format the tile index as a string in the format "hXXvYY"
+    # Format as MODIS tile string (e.g., h12v04)
     tile = "h$(lpad(horizontal_index, 2, '0'))v$(lpad(vertical_index, 2, '0'))"
+    return tile
+end
+
+export sinusoidal_to_MODLAND
+
+
+"""
+    latlon_to_MODLAND(lat, lon) -> String
 
-    # Return the formatted tile index
+Convert latitude and longitude (in degrees, WGS84) directly to the MODIS land tile index (e.g., "h12v04").
+
+Arguments:
+    lat: Latitude in degrees
+    lon: Longitude in degrees
+
+Returns:
+    String: MODIS tile index in the format "hXXvYY"
+"""
+function latlon_to_MODLAND(lat, lon)
+    # Convert lat/lon to sinusoidal coordinates
+    x, y = latlon_to_sinusoidal(lat, lon)
+    # Get MODIS tile index for the coordinates
+    tile = sinusoidal_to_MODLAND(x, y)
     return tile
 end
 
-export sinusoidal_to_modland
+export latlon_to_MODLAND
+
+
+"""
+    MODLAND_tiles_in_polygon(poly::ArchGDAL.IGeometry{ArchGDAL.wkbPolygon25D}) -> Set{String}
+    MODLAND_tiles_in_polygon(poly::AG.IGeometry{AG.wkbPolygon}) -> Set{String}
+
+Given a polygon geometry, return the set of MODIS land tile indices that contain the polygon's vertices.
+Note: Only the tiles containing the polygon's vertices are returned, not all tiles intersecting the polygon.
+
+Arguments:
+    poly: Polygon geometry (ArchGDAL.IGeometry)
 
-end
+Returns:
+    Set{String}: Set of MODIS tile indices (e.g., Set(["h12v04", ...]))
+"""
+function MODLAND_tiles_in_polygon(poly::ArchGDAL.IGeometry{ArchGDAL.wkbPolygon25D})::Set{String}
+    # Convert polygon to JSON, extract coordinates, and map to MODIS tiles
+    Set([latlon_to_MODLAND(lat, lon) for (lon, lat) in JSON.parse(AG.toJSON(poly))["coordinates"][1]])
+end
+
+function MODLAND_tiles_in_polygon(poly::AG.IGeometry{AG.wkbPolygon})::Set{String}
+    # Convert polygon to JSON, extract coordinates, and map to MODIS tiles
+    Set([latlon_to_MODLAND(lat, lon) for (lon, lat) in JSON.parse(AG.toJSON(poly))["coordinates"][1]])
+end
+
+export MODLAND_tiles_in_polygon
+
+end

test/runtests.jl

@@ -17,14 +17,14 @@ using Modland
 end
 
 # Test sinusoidal_to_modland
-@testset "sinusoidal_to_modland tests" begin
+@testset "sinusoidal_to_MODLAND tests" begin
     # Test with valid sinusoidal coordinates
-    @test sinusoidal_to_modland(0.0, 0.0) == "h18v09"
-    @test sinusoidal_to_modland(3.538204887918666e6, 5.003777338949354e6) == "h21v04"
+    @test sinusoidal_to_MODLAND(0.0, 0.0) == "h18v09"
+    @test sinusoidal_to_MODLAND(3.538204887918666e6, 5.003777338949354e6) == "h21v04"
 
     # Test with out-of-bounds sinusoidal coordinates
-    @test_throws ErrorException sinusoidal_to_modland(UPPER_LEFT_X_METERS - 1, 0.0)
-    @test_throws ErrorException sinusoidal_to_modland(LOWER_RIGHT_X_METERS + 1, 0.0)
-    @test_throws ErrorException sinusoidal_to_modland(0.0, LOWER_RIGHT_Y_METERS - 1)
-    @test_throws ErrorException sinusoidal_to_modland(0.0, UPPER_LEFT_Y_METERS + 1)
+    @test_throws ErrorException sinusoidal_to_MODLAND(Modland.UPPER_LEFT_X_METERS - 1, 0.0)
+    @test_throws ErrorException sinusoidal_to_MODLAND(Modland.LOWER_RIGHT_X_METERS + 1, 0.0)
+    @test_throws ErrorException sinusoidal_to_MODLAND(0.0, Modland.LOWER_RIGHT_Y_METERS - 1)
+    @test_throws ErrorException sinusoidal_to_MODLAND(0.0, Modland.UPPER_LEFT_Y_METERS + 1)
 end