moisture,latitude,longitude,sampled_at\n0.1000,50.451201,30.501118,2026-06-11T08:00:00Z\n0.8100,50.451802,30.502045,2026-06-11T08:01:00Z\n0.1842,50.452118,30.500752,2026-06-11T08:02:00Z\n<\/code><\/pre>\n\n\n\nTo display farm soil moisture level on map we need to turn those point measurements into a single-band GeoTIFF. At the final step of this post we will have GeoTIFF visualization on the map using leaflet like this :<\/p>\n\n\n\n![\"\"](\"https:\/\/mshaeri.com\/blog\/wp-content\/uploads\/2026\/06\/image-1024x817.png\")
<\/a>GeoTIFF visualiation on the map using leaflet<\/figcaption><\/figure>\n\n\n\nThe only issue with this example is that the latitude and longitude values were randomly generated, so the sample farm happens to appear inside an urban area\ud83d\ude05. You can think of it as a farm area.<\/p>\n\n\n\n
What you will read in this post is:<\/p>\n\n\n\n
\n- Generating GeoTIFF from CSV file data<\/a><\/li>\n\n\n\n
- Displaying the geotif on map using leaflet<\/a><\/li>\n<\/ol>\n\n\n\n
Also, here’s the GitHub repo for the full GeoTIFF-Generator project: https:\/\/github.com\/birddevelper\/GeoTIFF-Generator<\/a><\/p>\n\n\n\nGenerating GeoTIFF from CSV file data<\/h2>\n\n\n\n
Assume our CSV file has these columns:<\/p>\n\n\n\n
\nmoisture<\/code>: decimal value, usually 0.10<\/code> to 0.81<\/code> in the sample data<\/li>\n\n\n\nlatitude<\/code>: WGS84 latitude<\/li>\n\n\n\nlongitude<\/code>: WGS84 longitude<\/li>\n\n\n\nsampled_at<\/code>: ISO timestamp<\/li>\n<\/ul>\n\n\n\nThe output GeoTIFF <\/strong>would be a single-band byte raster. I choosed single band beacause it’s compact, easy to render in a browser, and still useful in GIS tools. Let’s start implementation:<\/p>\n\n\n\nThe Entry Point<\/h2>\n\n\n\n
In the main function of make_geotiff.py<\/code> (our main module) we get user arguments and create a RasterConfig<\/code>, and generates the GeoTIFF using make_geotiff<\/code> function. I’ll explain the make_geotiff <\/code>in later sections:<\/p>\n\n\n\nfrom __future__ import annotations\n\nimport argparse\nfrom pathlib import Path\n\nfrom core import RasterConfig\n\n\ndef default_output_dir() -> Path:\n return Path(__file__).resolve().parent \/ \"output\"\n\n\ndef build_parser() -> argparse.ArgumentParser:\n parser = argparse.ArgumentParser(\n description=\"Generate a single-band GeoTIFF from a moisture CSV file.\"\n )\n parser.add_argument(\"csv_path\", type=Path, help=\"Path to the input CSV file.\")\n parser.add_argument(\n \"--boundary-geojson\",\n dest=\"boundary_geojson_path\",\n type=Path,\n required=True,\n help=\"Path to a GeoJSON file containing the farm boundary polygon.\",\n )\n parser.add_argument(\n \"-o\",\n \"--output\",\n type=Path,\n default=default_output_dir(),\n help=\"Output directory or GeoTIFF file path. Defaults to the package output folder.\",\n )\n parser.add_argument(\n \"--resolution-m\",\n type=float,\n default=10.0,\n help=\"Approximate raster resolution in meters.\",\n )\n parser.add_argument(\n \"--max-distance-m\",\n type=float,\n default=None,\n help=\"Optional maximum distance from a sensor point before a pixel is marked NoData.\",\n )\n parser.add_argument(\n \"--nodata\",\n type=int,\n default=255,\n help=\"NoData byte value to write into the GeoTIFF. Must be 101-255.\",\n )\n return parser\n\n\ndef main(argv: list[str] | None = None) -> int:\n parser = build_parser()\n args = parser.parse_args(argv)\n\n from pipeline import make_geotiff\n\n config = RasterConfig(\n resolution_m=args.resolution_m,\n max_distance_m=args.max_distance_m,\n nodata_value=args.nodata,\n )\n output_path = make_geotiff(\n args.csv_path, args.boundary_geojson_path, args.output, config\n )\n print(f\"GeoTIFF written to: {output_path}\")\n return 0\n\n\nif __name__ == \"__main__\":\n raise SystemExit(main())\n<\/code><\/pre>\n\n\n\nIn the code above, the RasterConfig <\/code><\/strong>holds the GeoTIFF generation settings: pixel resolution, optional max interpolation distance, and NoData byte value.<\/p>\n\n\n\nCore Data And Input Loading<\/h2>\n\n\n\n
In core.py<\/code><\/strong> we declare data classes and input parsers.<\/p>\n\n\n\nfrom __future__ import annotations\n\nimport csv\nimport json\nfrom dataclasses import dataclass\nfrom datetime import datetime\nfrom pathlib import Path\n\n\n@dataclass(frozen=True)\nclass MoistureReading:\n moisture: float\n latitude: float\n longitude: float\n sampled_at: datetime\n\n\n@dataclass(frozen=True)\nclass RasterConfig:\n resolution_m: float = 10.0\n max_distance_m: float | None = None\n nodata_value: int = 255\n<\/code><\/pre>\n\n\n\nMoistureReading<\/code> <\/strong>data class stores one measured point, and RasterConfig<\/code> <\/strong>stores output settings, we created one instance of it in main function in make_geotiff.py<\/code> .<\/p>\n\n\n\ndef _parse_sampled_at(value: str) -> datetime:\n return datetime.fromisoformat(value.strip().replace(\"Z\", \"+00:00\"))\n\n\ndef read_moisture_csv(csv_path: Path) -> list[MoistureReading]:\n if not csv_path.exists():\n raise FileNotFoundError(f\"CSV file does not exist: {csv_path}\")\n\n readings: list[MoistureReading] = []\n with csv_path.open(\"r\", encoding=\"utf-8\", newline=\"\") as file:\n reader = csv.DictReader(file)\n required_columns = {\"moisture\", \"latitude\", \"longitude\", \"sampled_at\"}\n missing_columns = required_columns.difference(reader.fieldnames or [])\n if missing_columns:\n raise ValueError(f\"CSV file is missing required columns: {sorted(missing_columns)}\")\n\n for row in reader:\n moisture_raw = (row.get(\"moisture\") or \"\").strip()\n latitude_raw = (row.get(\"latitude\") or \"\").strip()\n longitude_raw = (row.get(\"longitude\") or \"\").strip()\n sampled_at_raw = (row.get(\"sampled_at\") or \"\").strip()\n\n if not moisture_raw or not latitude_raw or not longitude_raw or not sampled_at_raw:\n continue\n\n readings.append(\n MoistureReading(\n moisture=float(moisture_raw),\n latitude=float(latitude_raw),\n longitude=float(longitude_raw),\n sampled_at=_parse_sampled_at(sampled_at_raw),\n )\n )\n\n return readings\n<\/code><\/pre>\n\n\n\n_parse_sampled_at()<\/code> parses ISO timestamps and supports Z<\/code>. read_moisture_csv()<\/code> validates required columns and returns typed readings.<\/p>\n\n\n\ndef _extract_geojson_geometry(data: object) -> dict:\n if not isinstance(data, dict):\n raise ValueError(\"Invalid GeoJSON content.\")\n\n geometry_type = data.get(\"type\")\n if geometry_type == \"FeatureCollection\":\n features = data.get(\"features\") or []\n if not features:\n raise ValueError(\"GeoJSON FeatureCollection does not contain any features.\")\n return _extract_geojson_geometry(features[0])\n\n if geometry_type == \"Feature\":\n geometry = data.get(\"geometry\")\n if geometry is None:\n raise ValueError(\"GeoJSON Feature does not contain a geometry.\")\n return geometry\n\n return data\n\n\ndef load_boundary_geojson(geojson_path: Path):\n if not geojson_path.exists():\n raise FileNotFoundError(f\"Boundary GeoJSON file does not exist: {geojson_path}\")\n\n from shapely.geometry import MultiPolygon, Polygon, shape\n\n with geojson_path.open(\"r\", encoding=\"utf-8\") as file:\n data = json.load(file)\n\n geometry = _extract_geojson_geometry(data)\n polygon = shape(geometry)\n\n if polygon.is_empty:\n raise ValueError(f\"Boundary GeoJSON is empty: {geojson_path}\")\n\n if isinstance(polygon, MultiPolygon):\n polygon = max(polygon.geoms, key=lambda geom: geom.area)\n\n if not isinstance(polygon, Polygon):\n raise ValueError(\"Boundary GeoJSON must describe a Polygon or MultiPolygon.\")\n\n return polygon\n<\/code><\/pre>\n\n\n\n_extract_geojson_geometry()<\/strong><\/code> normalizes different valid GeoJSON shapes into one plain geometry object. If the input is a FeatureCollection, it takes the first feature; if it is a Feature, it extracts its geometry; if it is already a raw geometry like Polygon or MultiPolygon, it returns it unchanged. Its output is the geometry dictionary that Shapely can convert into a polygon for clipping the raster. load_boundary_geojson()<\/strong> reads the boundary GeoJSON file using _extract_geojson_geometry()<\/strong><\/code>, extracts its geometry, converts it into a Shapely shape, validates that it is not empty, and ensures the result is a usable farm boundary polygon. For simplicity, if the GeoJSON contains a MultiPolygon, it keeps the largest polygon. Its output is a Shapely Polygon, which the raster step uses to mask pixels outside the farm boundary as NoData.<\/p>\n\n\n\n
\n\n\n\nRaster Generation<\/h2>\n\n\n\n
In raster.py<\/strong><\/code> we convert point measurements into a georeferenced raster. A raster is a grid of pixels or cells where each cell stores a value. For a moisture GeoTIFF, the raster is like a spreadsheet laid over the farm. The raster remains in EPSG:4326, but interpolation and distance filtering operations need meter-like distances. The code approximates meters per degree at the farm latitude.<\/p>\n\n\n\nfrom __future__ import annotations\n\nimport math\nfrom dataclasses import dataclass\nfrom pathlib import Path\n\nimport numpy as np\nfrom osgeo import gdal, osr\nfrom scipy.spatial import cKDTree\nfrom core import MoistureReading\n\ndef meters_per_degree(latitude: float) -> tuple[float, float]:\n latitude_radians = math.radians(latitude)\n meters_per_degree_latitude = (\n 111132.92\n - 559.82 * math.cos(2 * latitude_radians)\n + 1.175 * math.cos(4 * latitude_radians)\n - 0.0023 * math.cos(6 * latitude_radians)\n )\n meters_per_degree_longitude = (\n 111412.84 * math.cos(latitude_radians)\n - 93.5 * math.cos(3 * latitude_radians)\n + 0.118 * math.cos(5 * latitude_radians)\n )\n return meters_per_degree_latitude, meters_per_degree_longitude\n\n\ndef to_metric_xy(\n longitudes: np.ndarray,\n latitudes: np.ndarray,\n reference_latitude: float,\n) -> tuple[np.ndarray, np.ndarray]:\n meters_lat, meters_lon = meters_per_degree(reference_latitude)\n x_coords = np.asarray(longitudes, dtype=np.float64) * meters_lon\n y_coords = np.asarray(latitudes, dtype=np.float64) * meters_lat\n return x_coords, y_coords\n\n\ndef degree_steps_for_resolution(resolution_m: float, reference_latitude: float) -> tuple[float, float]:\n meters_lat, meters_lon = meters_per_degree(reference_latitude)\n lon_step = resolution_m \/ meters_lon\n lat_step = resolution_m \/ meters_lat\n return lon_step, lat_step\n<\/code><\/pre>\n\n\n\nmeters_per_degree()<\/strong><\/code> estimates latitude and longitude scale at a reference latitude. to_metric_xy()<\/strong><\/code> converts lon\/lat arrays into approximate meter coordinates, and degree_steps_for_resolution()<\/strong><\/code> converts requested meter resolution into lon\/lat pixel steps.<\/p>\n\n\n\nGrid Definition<\/h3>\n\n\n\n@dataclass(frozen=True)\nclass GridSpec:\n width: int\n height: int\n lon_step: float\n lat_step: float\n geotransform: tuple[float, float, float, float, float, float]\n lon_grid: np.ndarray\n lat_grid: np.ndarray\n reference_latitude: float\n\n\ndef build_grid(\n bounds: tuple[float, float, float, float],\n resolution_m: float,\n reference_latitude: float,\n) -> GridSpec:\n min_lon, min_lat, max_lon, max_lat = bounds\n lon_step, lat_step = degree_steps_for_resolution(resolution_m, reference_latitude)\n\n width = max(1, int(math.ceil((max_lon - min_lon) \/ lon_step)))\n height = max(1, int(math.ceil((max_lat - min_lat) \/ lat_step)))\n\n lon_centers = min_lon + (np.arange(width) + 0.5) * lon_step\n lat_centers = max_lat - (np.arange(height) + 0.5) * lat_step\n lon_grid, lat_grid = np.meshgrid(lon_centers, lat_centers)\n geotransform = (min_lon, lon_step, 0.0, max_lat, 0.0, -lat_step)\n\n return GridSpec(\n width=width,\n height=height,\n lon_step=lon_step,\n lat_step=lat_step,\n geotransform=geotransform,\n lon_grid=lon_grid,\n lat_grid=lat_grid,\n reference_latitude=reference_latitude,\n )\n<\/code><\/pre>\n\n\n\nGridSpec<\/code> <\/strong>keeps raster size, pixel centers, and GDAL geotransform together. build_grid()<\/code><\/strong> creates a lon\/lat grid over the farm boundary bounds. It creates the empty raster grid that the moisture values will be written into. It takes the farm boundary bounds, the desired pixel size, like 10<\/code> meters, a reference latitude for meter-to-degree conversion, then it converts the requested meter resolution into longitude\/latitude steps, and calculates how many columns and rows the raster needs. Then, it creates longitude and latitude coordinates for each pixel center, and builds lon_grid<\/code> and lat_grid<\/code>, which represent every raster cell location. Finally it creates the GDAL geotransform<\/code>, which tells GIS tools where the raster sits on Earth.<\/p>\n\n\n\nContinuous Moisture Surface<\/h3>\n\n\n\n
For interpolation we use inverse-distance weighting. Each raster cell queries nearby sample points with cKDTree<\/strong><\/code>; closer samples receive higher weights. This way we will have smooth transation between tails in moisture visulization instead of disecret tiles.<\/p>\n\n\n\n![\"\"](\"https:\/\/mshaeri.com\/blog\/wp-content\/uploads\/2026\/06\/image-2.png\")
<\/a>Example of disecret tiles geotiff generation<\/figcaption><\/figure><\/div>\n\n\ndef interpolate_moisture(\n readings: list[MoistureReading],\n grid: GridSpec,\n *,\n neighbors: int = 12,\n power: float = 2.0,\n smoothing: float = 1e-9,\n) -> np.ndarray:\n longitudes = np.array([reading.longitude for reading in readings], dtype=np.float64)\n latitudes = np.array([reading.latitude for reading in readings], dtype=np.float64)\n moisture = np.array([reading.moisture for reading in readings], dtype=np.float64)\n\n point_x, point_y = to_metric_xy(longitudes, latitudes, grid.reference_latitude)\n grid_x, grid_y = to_metric_xy(grid.lon_grid.ravel(), grid.lat_grid.ravel(), grid.reference_latitude)\n\n query_points = np.column_stack([grid_x, grid_y])\n sample_points = np.column_stack([point_x, point_y])\n tree = cKDTree(sample_points)\n\n k = min(max(1, neighbors), len(readings))\n distances, indices = tree.query(query_points, k=k)\n distances = np.asarray(distances, dtype=np.float64)\n indices = np.asarray(indices)\n\n if k == 1:\n surface = moisture[indices]\n return surface.reshape(grid.lon_grid.shape).astype(np.float32)\n\n exact_matches = distances <= smoothing\n weights = 1.0 \/ np.power(distances + smoothing, power)\n weighted_moisture = moisture[indices] * weights\n weighted_sum = weighted_moisture.sum(axis=1)\n weight_total = weights.sum(axis=1)\n surface = np.divide(\n weighted_sum,\n weight_total,\n out=np.zeros_like(weighted_sum),\n where=weight_total > 0,\n )\n\n if np.any(exact_matches):\n exact_rows = np.where(np.any(exact_matches, axis=1))[0]\n exact_cols = np.argmax(exact_matches[exact_rows], axis=1)\n surface[exact_rows] = moisture[indices[exact_rows, exact_cols]]\n\n return surface.reshape(grid.lon_grid.shape).astype(np.float32)\n<\/code><\/pre>\n\n\n\nFunction summary:<\/p>\n\n\n\n
\ninterpolate_moisture()<\/code> creates a continuous decimal moisture grid.<\/li>\n\n\n\nneighbors<\/code> limits how many sample points affect each cell.<\/li>\n\n\n\npower<\/code> controls how quickly influence drops with distance.<\/li>\n\n\n\nsmoothing<\/code> prevents division by zero and supports exact point matches.<\/li>\n<\/ul>\n\n\n\nBoundary And Distance Masks<\/h3>\n\n\n\n
Displaying moisture levels outside the farm boundary is undesirable, so we need to clip the map to the farm’s perimeter. The farm boundary is available from a GeoJSON file provided via command-line arguments.<\/p>\n\n\n\n
def clip_to_boundary(grid_values: np.ndarray, grid: GridSpec, boundary_polygon) -> np.ndarray:\n from shapely import contains_xy\n\n mask = contains_xy(boundary_polygon, grid.lon_grid, grid.lat_grid)\n clipped = np.array(grid_values, copy=True)\n clipped[~mask] = np.nan\n return clipped\n\n\ndef limit_by_distance(\n readings: list[MoistureReading],\n grid_values: np.ndarray,\n grid: GridSpec,\n max_distance_m: float,\n) -> np.ndarray:\n point_lons = np.array([reading.longitude for reading in readings], dtype=np.float64)\n point_lats = np.array([reading.latitude for reading in readings], dtype=np.float64)\n point_x, point_y = to_metric_xy(point_lons, point_lats, grid.reference_latitude)\n\n grid_x, grid_y = to_metric_xy(grid.lon_grid.ravel(), grid.lat_grid.ravel(), grid.reference_latitude)\n tree = cKDTree(np.column_stack([point_x, point_y]))\n distances, _ = tree.query(np.column_stack([grid_x, grid_y]))\n\n limited = np.array(grid_values, copy=True).ravel()\n limited[distances > max_distance_m] = np.nan\n return limited.reshape(grid_values.shape)\n<\/code><\/pre>\n\n\n\nIn the code above, the clip_to_boundary()<\/code><\/strong> sets cells outside the farm polygon to np.nan<\/code>, and limit_by_distance()<\/strong><\/code> optionally sets cells too far from samples to np.nan<\/code>. For example, if the nearest sample point is 80 meters away and you set –max-distance-m 30, that cell becomes np.nan (NoData) in the GeoTIFF.<\/p>\n\n\n\nGeoTIFF Writer<\/h3>\n\n\n\n
The writer converts decimal moisture to byte percent values. GDAL writes the raster as GDT_Byte<\/code>.<\/p>\n\n\n\ndef write_single_band_geotiff(\n output_path: Path,\n grid_values: np.ndarray,\n grid: GridSpec,\n nodata_value: int,\n) -> Path:\n if not 101 <= nodata_value <= 255:\n raise ValueError(\"UInt8 GeoTIFF NoData value must be between 101 and 255.\")\n\n output_path.parent.mkdir(parents=True, exist_ok=True)\n\n driver = gdal.GetDriverByName(\"GTiff\")\n dataset = driver.Create(\n str(output_path),\n grid.width,\n grid.height,\n 1,\n gdal.GDT_Byte,\n options=[\"COMPRESS=LZW\"],\n )\n if dataset is None:\n raise RuntimeError(f\"Could not create GeoTIFF at {output_path}\")\n\n dataset.SetGeoTransform(grid.geotransform)\n spatial_reference = osr.SpatialReference()\n spatial_reference.ImportFromEPSG(4326)\n dataset.SetProjection(spatial_reference.ExportToWkt())\n\n band = dataset.GetRasterBand(1)\n band.SetNoDataValue(nodata_value)\n\n percent_values = np.where(\n np.isfinite(grid_values),\n np.rint(np.clip(grid_values * 100.0, 0.0, 100.0)),\n nodata_value,\n ).astype(np.uint8)\n band.WriteRaster(\n 0,\n 0,\n grid.width,\n grid.height,\n percent_values.tobytes(order=\"C\"),\n buf_xsize=grid.width,\n buf_ysize=grid.height,\n buf_type=gdal.GDT_Byte,\n )\n band.FlushCache()\n\n dataset.FlushCache()\n dataset = None\n return output_path\n\n\n<\/code><\/pre>\n\n\n\n\nwrite_single_band_geotiff<\/strong>()<\/code> writes one UInt8<\/code> band with LZW compression.<\/li>\n\n\n\n- Valid decimal moisture is scaled by
* 100<\/code>, rounded, and clipped to 0..100<\/code>.<\/li>\n\n\n\n- Missing cells become
nodata_value<\/code>, normally 255<\/code>.<\/li>\n\n\n\nWriteRaster<\/strong>()<\/code> avoids requiring GDAL’s gdal_array<\/code> <\/strong>extension.<\/li>\n<\/ul>\n\n\n\nPipeline<\/h2>\n\n\n\n
Finally, pipeline.py<\/code> connects all pieces.<\/p>\n\n\n\nfrom __future__ import annotations\n\nfrom pathlib import Path\n\nfrom core import RasterConfig, load_boundary_geojson, read_moisture_csv\nfrom raster import build_grid, clip_to_boundary, interpolate_moisture, limit_by_distance, write_single_band_geotiff\n\n\ndef default_output_dir() -> Path:\n return Path(__file__).resolve().parent \/ \"output\"\n\n\ndef resolve_output_path(csv_path: Path, output_target: Path | None) -> Path:\n if output_target is None:\n return default_output_dir() \/ f\"{csv_path.stem}_moisture.tif\"\n\n if output_target.suffix.lower() in {\".tif\", \".tiff\"}:\n return output_target\n\n return output_target \/ f\"{csv_path.stem}_moisture.tif\"\n\n\ndef make_geotiff(\n csv_path: Path,\n boundary_geojson_path: Path,\n output_target: Path | None = None,\n config: RasterConfig | None = None,\n) -> Path:\n job_config = config or RasterConfig()\n readings = read_moisture_csv(csv_path)\n if not readings:\n raise ValueError(\"No valid moisture readings were found in the CSV file.\")\n\n boundary = load_boundary_geojson(boundary_geojson_path)\n reference_latitude = float(sum(reading.latitude for reading in readings) \/ len(readings))\n\n grid = build_grid(boundary.bounds, job_config.resolution_m, reference_latitude)\n interpolated = interpolate_moisture(readings, grid)\n clipped = clip_to_boundary(interpolated, grid, boundary)\n\n if job_config.max_distance_m is not None:\n clipped = limit_by_distance(readings, clipped, grid, job_config.max_distance_m)\n\n output_path = resolve_output_path(csv_path, output_target)\n return write_single_band_geotiff(output_path, clipped, grid, job_config.nodata_value)\n\n\n\n\n<\/code><\/pre>\n\n\n\nmake_geotiff()<\/strong><\/code> is the main pipeline function that turns the input CSV and farm boundary GeoJSON into the final GeoTIFF file. To sum up the flow, the generation order is like following:<\/p>\n\n\n\n\n- read CSV<\/li>\n\n\n\n
- load boundary<\/li>\n\n\n\n
- compute average latitude for meter conversion<\/li>\n\n\n\n
- build grid<\/li>\n\n\n\n
- interpolate moisture<\/li>\n\n\n\n
- clip to boundary<\/li>\n\n\n\n
- optionally apply max-distance mask<\/li>\n\n\n\n
- write GeoTIFF<\/li>\n<\/ol>\n\n\n\n
To generate GeoTIFF, run the make_geotiff.py<\/code> like this:<\/p>\n\n\n\npython make_geotiff.py input\/farm_moisture_200.csv --boundary-geojson input\/farm_boundary.geojson\n<\/code><\/pre>\n\n\n\n<\/p>\n\n\n\n
Displaying the GeoTIFF on map using leaflet<\/h2>\n\n\n\n
The GeoTIFF stays single-band. So in the browser we read its pixel value and maps it to a color.<\/p>\n\n\n\n
HTML<\/h3>\n\n\n\n<!DOCTYPE html>\n<html lang=\"en\">\n<head>\n <meta charset=\"UTF-8\" \/>\n <meta name=\"viewport\" content=\"width=device-width, initial-scale=1.0\" \/>\n <title>Farm Moisture Map<\/title>\n <link rel=\"stylesheet\" href=\"https:\/\/unpkg.com\/leaflet@1.9.4\/dist\/leaflet.css\" \/>\n <style>\n html, body, #map {\n margin: 0;\n width: 100%;\n height: 100%;\n font-family: Arial, sans-serif;\n }\n <\/style>\n<\/head>\n<body>\n <div id=\"map\"><\/div>\n <script src=\"https:\/\/unpkg.com\/leaflet@1.9.4\/dist\/leaflet.js\"><\/script>\n <script src=\"https:\/\/unpkg.com\/geotiff\/dist-browser\/geotiff.js\"><\/script>\n <script src=\"https:\/\/cdn.jsdelivr.net\/npm\/georaster@1.6.0\/dist\/georaster.browser.bundle.min.js\"><\/script>\n <script src=\"https:\/\/unpkg.com\/georaster-layer-for-leaflet\/dist\/georaster-layer-for-leaflet.min.js\"><\/script>\n <script src=\"https:\/\/unpkg.com\/chroma-js@2.4.2\/chroma.min.js\"><\/script>\n <script src=\"app.js\"><\/script>\n<\/body>\n<\/html>\n<\/code><\/pre>\n\n\n\nJavaScript<\/h3>\n\n\n\nconst GEOTIFF_URL = \"output\/farm_moisture_200_moisture.tif\";\nconst MIN_VALUE = 0;\nconst MAX_VALUE = 100;\n\nconst map = L.map(\"map\", { zoomControl: true }).setView([0, 0], 2);\n\nL.tileLayer(\"https:\/\/tile.openstreetmap.org\/{z}\/{x}\/{y}.png\", {\n attribution: \"© OpenStreetMap contributors\",\n maxZoom: 19,\n}).addTo(map);\n\nconst moistureScale = chroma\n .scale([\"#7b4f2a\", \"#d9c27d\", \"#4caf50\"])\n .domain([MIN_VALUE, MAX_VALUE]);\n\nasync function loadMoistureRaster() {\n try {\n const response = await fetch(GEOTIFF_URL);\n if (!response.ok) {\n throw new Error(`Failed to load GeoTIFF: ${response.status}`);\n }\n\n const arrayBuffer = await response.arrayBuffer();\n const georaster = await parseGeoraster(arrayBuffer);\n const noData = georaster.noDataValue;\n\n const rasterLayer = new GeoRasterLayer({\n georaster,\n opacity: 0.85,\n resolution: 256,\n pixelValuesToColorFn: (values) => {\n const moisture = values?.[0];\n\n if (\n moisture === undefined ||\n moisture === null ||\n Number.isNaN(moisture) ||\n moisture === noData\n ) {\n return undefined;\n }\n\n const clamped = Math.max(MIN_VALUE, Math.min(MAX_VALUE, moisture));\n return moistureScale(clamped).hex();\n },\n });\n\n rasterLayer.addTo(map);\n map.fitBounds(rasterLayer.getBounds());\n } catch (error) {\n console.error(error);\n alert(\"Failed to load GeoTIFF. Check the console for details.\");\n }\n}\n\nloadMoistureRaster();\n\n\n<\/code><\/pre>\n\n\n\nIn the javascript after we load the GeoTIFF bytes using fetch()<\/code>, we parse it by parseGeoraster()<\/code>. Then GeoRasterLayer<\/code> renders the raster on Leaflet. The pixelValuesToColorFn()<\/code> function skips NoData and maps 0..100<\/code> moisture to a corresponding color. The fitBounds()<\/code> zooms the map to the raster extent. So when page loads, it automatically zoom in the farm area.<\/p>\n\n\n\nServe the HTML over a HTTP server like apache, or ngnix, not file:\/\/<\/code>, so the browser can fetch the GeoTIFF and openstreet, otherwise you face CORS error. BTW, you can find source code of the project in this github repository: https:\/\/github.com\/birddevelper\/GeoTIFF-Generator<\/a><\/p>\n\n\n\nI hope you found this post helpful. If you enjoyed it or have any questions, leave a comment, I\u2019d love to hear your thoughts.<\/p>\n","protected":false},"excerpt":{"rendered":"
Assume we have a collection of soil-moisture measurements from a farm. Each measurement has a decimal moisture value, latitude, longitude, and sample time. We want … <\/p>\n","protected":false},"author":1,"featured_media":2536,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[1,41],"tags":[364,353,363,351,361,360,358,352,359,39,354,362,357],"_links":{"self":[{"href":"https:\/\/mshaeri.com\/blog\/wp-json\/wp\/v2\/posts\/2510"}],"collection":[{"href":"https:\/\/mshaeri.com\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/mshaeri.com\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/mshaeri.com\/blog\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/mshaeri.com\/blog\/wp-json\/wp\/v2\/comments?post=2510"}],"version-history":[{"count":4,"href":"https:\/\/mshaeri.com\/blog\/wp-json\/wp\/v2\/posts\/2510\/revisions"}],"predecessor-version":[{"id":2566,"href":"https:\/\/mshaeri.com\/blog\/wp-json\/wp\/v2\/posts\/2510\/revisions\/2566"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/mshaeri.com\/blog\/wp-json\/wp\/v2\/media\/2536"}],"wp:attachment":[{"href":"https:\/\/mshaeri.com\/blog\/wp-json\/wp\/v2\/media?parent=2510"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/mshaeri.com\/blog\/wp-json\/wp\/v2\/categories?post=2510"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/mshaeri.com\/blog\/wp-json\/wp\/v2\/tags?post=2510"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}
![\"\"](\"https:\/\/mshaeri.com\/blog\/wp-content\/uploads\/2026\/06\/image.png\")
<\/a>