🧭 Goal
Compute NDVI values for each pixel by aligning RED and NIR GeoTIFFs spatially.

✅ Assumptions
You already have:

Two GeoTIFFs: red.tif, nir.tif

Each image:

width, height

tiepoint: top-left world coordinate [x, y]

pixelScale: [xSize, ySize] (usually [30, 30])

1D pixel value arrays:

redPixels: Uint16Array, length = width × height

nirPixels: Uint16Array, same

CRS is the same for both (e.g., UTM), and there's no rotation/skew

🧠 Conceptual Plan
Choose one image as the reference grid (e.g. RED).

For each pixel (i, j) in RED:

Convert (i, j) to world coordinates (x, y)

Convert (x, y) to NIR pixel (i', j')

If in bounds, get both RED and NIR values

Compute NDVI = (NIR - RED) / (NIR + RED)

Store NDVI

🧱 Code Implementation Outline
1. Helper: Get world coordinate from raster index
js
Copier
Modifier
function rasterToWorld(i, j, tiepoint, scale) {
  const x = tiepoint.x + i * scale.x;
  const y = tiepoint.y - j * scale.y;  // minus because y increases downward in raster
  return { x, y };
}
2. Helper: Get raster index from world coordinate
js
Copier
Modifier
function worldToRaster(x, y, tiepoint, scale) {
  const i = Math.floor((x - tiepoint.x) / scale.x);
  const j = Math.floor((tiepoint.y - y) / scale.y);
  return { i, j };
}
3. Main loop: Compute NDVI
js
Copier
Modifier
const ndviArray = new Float32Array(redWidth * redHeight);

for (let j = 0; j < redHeight; j++) {
  for (let i = 0; i < redWidth; i++) {
    const redIdx = j * redWidth + i;
    const redDN = redPixels[redIdx];

    // 1. RED pixel → world coordinates
    const { x, y } = rasterToWorld(i, j, redTiepoint, redScale);

    // 2. World coordinates → NIR pixel
    const { i: iNir, j: jNir } = worldToRaster(x, y, nirTiepoint, nirScale);

    // 3. Skip if outside NIR image bounds
    if (
      iNir < 0 || iNir >= nirWidth ||
      jNir < 0 || jNir >= nirHeight
    ) {
      ndviArray[redIdx] = NaN;
      continue;
    }

    const nirIdx = jNir * nirWidth + iNir;
    const nirDN = nirPixels[nirIdx];

    // 4. Compute NDVI
    const denom = nirDN + redDN;
    const ndvi = denom === 0 ? 0 : (nirDN - redDN) / denom;

    ndviArray[redIdx] = ndvi;
  }
}
4. Post-process NDVI (optional)
You can:

Clamp or normalize NDVI for grayscale or color visualization.

Write out to a canvas or export to a new GeoTIFF.

Example: grayscale mapping for canvas (from NDVI range -1..1 to 0..255)

js
Copier
Modifier
const grayscale = new Uint8ClampedArray(redWidth * redHeight * 4);
for (let i = 0; i < ndviArray.length; i++) {
  const ndvi = ndviArray[i];
  const gray = isNaN(ndvi) ? 0 : Math.round(((ndvi + 1) / 2) * 255);
  grayscale[i * 4 + 0] = gray;
  grayscale[i * 4 + 1] = gray;
  grayscale[i * 4 + 2] = gray;
  grayscale[i * 4 + 3] = 255;
}
🔄 What If Images Are in Opposite Orientations?
You may need to flip Y if one image’s pixel_height is negative and the other isn’t. Always double-check this by printing out a few worldToRaster and rasterToWorld mappings.

✅ You Now Have:
A fully aligned NDVI computation pipeline

No assumptions about tile count or pixel alignment

Easy integration with canvas or export to raster

Would you like:

Help writing out a new NDVI GeoTIFF?

A geotiff.js wrapper that automates alignment?

NDVI color mapping for vegetation classification (e.g. healthy vs stressed)?



Vous avez dit :
For each pixel (i, j) in RED

It ends up as a 1D array though, there is only i

ChatGPT a dit :
Great observation — and yes, you’re absolutely right.

While the NDVI processing logically operates on 2D pixels (row j, column i), when you're working with the raw data from GeoTIFFs or tiles, it’s typically a flattened 1D array. So instead of looping over (i, j), you loop over a single index idx and calculate i and j from it.

✅ How to Work with the 1D Array
Let’s say you have:

width = redWidth

height = redHeight

redPixels.length === width * height

Then for any 1D index, you compute (i, j) like this:

js
Copier
Modifier
const j = Math.floor(index / width); // row (Y)
const i = index % width;             // column (X)
So this lets you still treat the flattened array like a 2D grid — just via math.

🔁 Revised Loop Using 1D Indexing
Here’s the NDVI loop fully flattened for your current setup:

js
Copier
Modifier
const ndviArray = new Float32Array(redPixels.length);

for (let index = 0; index < redPixels.length; index++) {
  const redDN = redPixels[index];

  // Convert 1D index to 2D pixel location in RED
  const j = Math.floor(index / redWidth);
  const i = index % redWidth;

  // Get world coordinates for RED pixel
  const { x, y } = rasterToWorld(i, j, redTiepoint, redScale);

  // Convert world coords to NIR (i', j') index
  const { i: iNir, j: jNir } = worldToRaster(x, y, nirTiepoint, nirScale);

  // Bounds check for NIR
  if (
    iNir < 0 || iNir >= nirWidth ||
    jNir < 0 || jNir >= nirHeight
  ) {
    ndviArray[index] = NaN;
    continue;
  }

  const nirIndex = jNir * nirWidth + iNir;
  const nirDN = nirPixels[nirIndex];

  const denom = nirDN + redDN;
  const ndvi = denom === 0 ? 0 : (nirDN - redDN) / denom;

  ndviArray[index] = ndvi;
}
🔧 Recap
Task	Method
Loop over RED pixels	Use a flat for (let idx = 0; idx < width * height; idx++)
Convert 1D to (i, j)	j = Math.floor(idx / width), i = idx % width
Use world alignment	Same as before, using raster ↔ world conversions

Let me know if you want this to support bilinear interpolation (for subpixel alignment), or output to canvas or a new raster!



Option D (GPU-based NDVI with WebGL shaders) is actually a very smart move if:

You're working in the browser

You want fast NDVI rendering (especially on large rasters)

You’re okay handling GPU texture logistics and GLSL shaders

Let’s walk through exactly how to build it. You don’t need a full GIS engine — just WebGL, RED and NIR textures, and a small shader.

🎯 Goal
Use WebGL to:

Load two grayscale GeoTIFF bands (RED and NIR) as textures

Send them to a fragment shader

Compute NDVI per pixel in parallel on GPU

Render result to canvas

🧱 Step-by-Step Blueprint
✅ Step 1: Preprocess GeoTIFFs into same size in JS
You already:

Decompress RED and NIR GeoTIFFs

Get the pixels as 1D Uint16Arrays

Know their width/height

For WebGL:

Convert those arrays into Uint8Arrays (or Float32Arrays if using OES_texture_float)

Normalize DNs to [0, 1] if needed

js
Copier
Modifier
function normalizeToUint8(srcArray, maxValue = 65535) {
  const out = new Uint8Array(srcArray.length);
  for (let i = 0; i < srcArray.length; i++) {
    out[i] = Math.min(255, Math.round((srcArray[i] / maxValue) * 255));
  }
  return out;
}
✅ Step 2: Upload RED and NIR to WebGL as grayscale textures
Create two LUMINANCE or R8 textures:

js
Copier
Modifier
const redTex = gl.createTexture();
gl.bindTexture(gl.TEXTURE_2D, redTex);
gl.texImage2D(
  gl.TEXTURE_2D, 0, gl.LUMINANCE, width, height, 0,
  gl.LUMINANCE, gl.UNSIGNED_BYTE, redData
);

const nirTex = gl.createTexture();
gl.bindTexture(gl.TEXTURE_2D, nirTex);
gl.texImage2D(
  gl.TEXTURE_2D, 0, gl.LUMINANCE, width, height, 0,
  gl.LUMINANCE, gl.UNSIGNED_BYTE, nirData
);
✅ Step 3: Fragment shader that computes NDVI
glsl
Copier
Modifier
precision mediump float;

uniform sampler2D redTex;
uniform sampler2D nirTex;
varying vec2 vTexCoord;

void main() {
  float red = texture2D(redTex, vTexCoord).r;
  float nir = texture2D(nirTex, vTexCoord).r;

  float denom = nir + red;
  float ndvi = denom == 0.0 ? 0.0 : (nir - red) / denom;

  // Map NDVI [-1, 1] → grayscale [0, 1]
  float gray = (ndvi + 1.0) / 2.0;

  gl_FragColor = vec4(vec3(gray), 1.0);
}
✅ Step 4: Setup WebGL quad and pipeline
Standard full-screen quad + vTexCoord passthrough from vertex shader:

glsl
Copier
Modifier
attribute vec2 aPosition;
varying vec2 vTexCoord;

void main() {
  vTexCoord = (aPosition + 1.0) / 2.0; // from clip space to [0,1]
  gl_Position = vec4(aPosition, 0.0, 1.0);
}
Then draw once and let the GPU compute the NDVI for every pixel instantly.

✅ Step 5: Optional – Download result or apply colormap
You can:

Read the rendered canvas using gl.readPixels() if needed

Apply NDVI-based colormap in shader

Overlay on Leaflet, Mapbox, etc. via WebGL canvas