<!DOCTYPE html>
<html lang="en">
<head>
  <meta charset="UTF-8">
  <title>USA LiDAR Extractor - 3DEP + PLSS</title>
  <meta name="viewport" content="width=device-width, initial-scale=1.0">
  <!-- OpenLayers CSS -->
  <link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/ol@7.3.0/ol.css" />
  <style>
    html, body {
      margin: 0; padding: 0; height: 100%;
      font-family: Arial, sans-serif;
    }
    #map { width:100%; height:80%; }
    #controls {
      width:100%; height:20%; padding:10px;
      background:#f3f3f3; box-sizing:border-box;
    }
    .control-button {
      background:#007BFF; color:#fff; padding:8px 12px;
      border:none; border-radius:4px; cursor:pointer;
      margin-right:10px;
    }
    .control-button:hover { background:#0056b3; }

    /* Spinner overlay styling */
    #spinnerOverlay {
      position: fixed; top: 0; left: 0;
      width:100%; height:100%;
      background-color: rgba(0,0,0,0.4);
      display:none; /* hidden by default */
      align-items:center; justify-content:center;
      z-index:9999; /* top */
    }
    #spinnerOverlay .spinnerMessage {
      color:#fff; font-size:1.2em; text-align:center; margin-top:10px;
    }
    .loader {
      border:8px solid #f3f3f3; border-top:8px solid #3498db;
      border-radius:50%; width:60px; height:60px;
      animation:spin 1s linear infinite; margin:auto;
    }
    @keyframes spin {
      0%{transform:rotate(0deg);}
      100%{transform:rotate(360deg);}
    }
  </style>
</head>
<body>

<div id="map"></div>

<div id="controls">
  <button class="control-button" id="drawAOIButton">Draw AOI Polygon</button>
  <button class="control-button" id="clearAOIButton">Clear AOI</button>

  <button class="control-button" id="download4326Button">Download LiDAR - EPSG:4326</button>
  <button class="control-button" id="download3857Button">Download LiDAR - EPSG:3857</button>

  <p style="margin-top:8px; color:#666;">
    1) Click "Draw AOI Polygon" to define your area.<br/>
    2) Double-click (or click first vertex) to finish.<br/>
    3) For smaller AOIs (&lt;2.5×2.5 km) we request ~1 m resolution.<br/>
    4) Larger AOIs are clamped at 2500 px in the bigger dimension.<br/>
    5) If area &gt; 200 km² we tile automatically (~200 km² each).
  </p>
</div>

<!-- Spinner Overlay -->
<div id="spinnerOverlay">
  <div>
    <div class="loader"></div>
    <div class="spinnerMessage">Downloading… please wait.</div>
  </div>
</div>

<!-- OpenLayers -->
<script src="https://cdn.jsdelivr.net/npm/ol@7.3.0/dist/ol.js"></script>

<script>
//////////////////
// Spinner
//////////////////
function showSpinner(){
  document.getElementById('spinnerOverlay').style.display='flex';
}
function hideSpinner(){
  document.getElementById('spinnerOverlay').style.display='none';
}

//////////////////
// Layers & Map
//////////////////

// US Imagery+Topo
const usImageryTopoLayer = new ol.layer.Tile({
  source: new ol.source.TileWMS({
    url:'https://basemap.nationalmap.gov/arcgis/services/USGSImageryTopo/MapServer/WMSServer',
    params:{
      'LAYERS':'0','TILED':true,
      'FORMAT':'image/png','VERSION':'1.3.0'
    },
    crossOrigin:'anonymous'
  }),
  visible:true
});

// US Hillshade
const usHillshadeLayer = new ol.layer.Tile({
  source: new ol.source.TileWMS({
    url:'https://elevation.nationalmap.gov/arcgis/services/3DEPElevation/ImageServer/WMSServer',
    params:{
      'LAYERS':'3DEPElevation:Hillshade Gray','TILED':true,
      'FORMAT':'image/png','VERSION':'1.3.0'
    },
    crossOrigin:'anonymous'
  }),
  opacity:0.5,
  visible:true
});

// PLSS
const plssLayer = new ol.layer.Tile({
  source: new ol.source.TileWMS({
    url:'https://gis.blm.gov/arcgis/services/Cadastral/BLM_Natl_PLSS_CadNSDI/MapServer/WMSServer',
    params:{
      'LAYERS':'1,2','TILED':false,
      'FORMAT':'image/png','VERSION':'1.3.0'
    },
    crossOrigin:'anonymous'
  }),
  visible:true
});

// View
const usView = new ol.View({
  center: ol.proj.fromLonLat([-96,38]),
  zoom:5
});

// Map
const map = new ol.Map({
  target:'map',
  layers:[
    usImageryTopoLayer,
    usHillshadeLayer,
    plssLayer
  ],
  view:usView
});

// AOI Vector
const drawSource= new ol.source.Vector();
const drawLayer= new ol.layer.Vector({
  source:drawSource,
  style: new ol.style.Style({
    stroke:new ol.style.Stroke({ color:'rgba(255,0,0,0.8)', width:2 }),
    fill:new ol.style.Fill({ color:'rgba(255,0,0,0.2)' })
  })
});
map.addLayer(drawLayer);

//////////////////
// Draw
//////////////////
let drawInteraction=null;
document.getElementById('drawAOIButton').addEventListener('click',()=>{
  if(drawInteraction){
    map.removeInteraction(drawInteraction);
    drawInteraction=null;
  }
  drawSource.clear();

  drawInteraction=new ol.interaction.Draw({
    source:drawSource,
    type:'Polygon'
  });
  map.addInteraction(drawInteraction);

  // Turn off double-click zoom while drawing
  const dblClickZoom= map.getInteractions().getArray().find(i=>i instanceof ol.interaction.DoubleClickZoom);
  if(dblClickZoom) dblClickZoom.setActive(false);

  drawInteraction.on('drawend',()=>{
    setTimeout(()=>{
      map.removeInteraction(drawInteraction);
      drawInteraction=null;
      if(dblClickZoom) dblClickZoom.setActive(true);
      alert('AOI polygon drawn. Now pick a download button.');
    },100);
  });
});

// Clear
document.getElementById('clearAOIButton').addEventListener('click',()=>{
  drawSource.clear();
});

//////////////////
// Tiling config
//////////////////
const MAX_AREA_SQKM=200; // tile if >200
const TILE_SIZE_METERS= Math.sqrt(MAX_AREA_SQKM*1e6); // ~14142m side

////////////////////////////
// EPSG:4326 Download
////////////////////////////
document.getElementById('download4326Button').addEventListener('click', async()=>{
  if(drawSource.getFeatures().length===0){
    alert('No AOI polygon drawn!');
    return;
  }
  const feature= drawSource.getFeatures()[0];
  const geom= feature.getGeometry();
  if(geom.getType()!=='Polygon'){
    alert('AOI must be polygon.');
    return;
  }

  // measure area in EPSG:3857
  const extent3857= geom.getExtent();
  const widthM = Math.abs(extent3857[2]-extent3857[0]);
  const heightM= Math.abs(extent3857[3]-extent3857[1]);
  const areaM2= widthM*heightM;

  if(areaM2 <= MAX_AREA_SQKM*1e6){
    await downloadSingleCoverage4326(geom);
  } else {
    await downloadTiledCoverage4326(geom, areaM2);
  }
});

async function downloadSingleCoverage4326(geom){
  // transform to EPSG:4326
  const ext3857= geom.getExtent();
  const ext4326= ol.proj.transformExtent(ext3857,'EPSG:3857','EPSG:4326');
  const [minLon,minLat,maxLon,maxLat] = ext4326;

  // figure out bounding box in meters (EPSG:3857) to see if <2.5km
  const xMinMax= ol.proj.transform([minLon,minLat], 'EPSG:4326','EPSG:3857');
  const xMaxMax= ol.proj.transform([maxLon,maxLat], 'EPSG:4326','EPSG:3857');
  const xSpan= Math.abs(xMaxMax[0] - xMinMax[0]);
  const ySpan= Math.abs(xMaxMax[1] - xMinMax[1]);

  // If both xSpan,ySpan <2500 => 1m resolution => width=round(xSpan), height=round(ySpan)
  // else clamp to 2500 in bigger dimension
  let [width,height]= pickPixelDimensions(xSpan,ySpan);

  showSpinner();
  const wcsUrl= `https://elevation.nationalmap.gov/arcgis/services/3DEPElevation/ImageServer/WCSServer`
    + `?service=WCS&version=1.0.0&request=GetCoverage`
    + `&coverage=DEP3Elevation`
    + `&CRS=EPSG:4326`
    + `&BBOX=${minLon},${minLat},${maxLon},${maxLat}`
    + `&width=${width}&height=${height}`
    + `&format=GeoTIFF`;

  try{
    const resp= await fetch(wcsUrl);
    if(!resp.ok){
      hideSpinner();
      alert(`WCS request failed: ${resp.statusText}`);
      return;
    }
    const arrayBuf= await resp.arrayBuffer();

    hideSpinner();

    const blob= new Blob([arrayBuf], { type:'image/tiff' });
    const urlObj= URL.createObjectURL(blob);

    const link= document.createElement('a');
    link.href= urlObj;
    link.download= 'lidar_usa_4326.tif';
    link.click();

    URL.revokeObjectURL(urlObj);
    alert('LiDAR (EPSG:4326) downloaded!');
  }catch(err){
    hideSpinner();
    console.error(err);
    alert('Error: '+ err);
  }
}

async function downloadTiledCoverage4326(geom, areaM2){
  alert(`AOI ~${(areaM2/1e6).toFixed(2)} km². Splitting into ~200 km² tiles.`);

  const ext3857= geom.getExtent();
  const [xMin,yMin,xMax,yMax]= ext3857;

  const totalW= Math.abs(xMax - xMin);
  const totalH= Math.abs(yMax - yMin);

  const tileCountX= Math.ceil(totalW / TILE_SIZE_METERS);
  const tileCountY= Math.ceil(totalH / TILE_SIZE_METERS);
  const totalTiles= tileCountX * tileCountY;

  alert(`Tiling into ${tileCountX} × ${tileCountY} = ${totalTiles} tiles...`);

  let tileIndex=0;
  for(let row=0; row<tileCountY; row++){
    for(let col=0; col<tileCountX; col++){
      tileIndex++;
      const subXMin= xMin + col*TILE_SIZE_METERS;
      const subXMax= Math.min(subXMin+TILE_SIZE_METERS, xMax);
      const subYMin= yMin + row*TILE_SIZE_METERS;
      const subYMax= Math.min(subYMin+TILE_SIZE_METERS, yMax);

      // transform to EPSG:4326
      const sub4326= ol.proj.transformExtent([subXMin,subYMin,subXMax,subYMax],'EPSG:3857','EPSG:4326');
      await downloadOneTile4326(sub4326, tileIndex, totalTiles);
      await new Promise(r=>setTimeout(r,500));
    }
  }
  alert('All tiles downloaded (EPSG:4326).');
}

async function downloadOneTile4326(bbox4326, tileIndex, totalTiles){
  const [minLon,minLat,maxLon,maxLat]= bbox4326;

  // figure out bounding box in meters for sizing
  const minXY= ol.proj.transform([minLon,minLat],'EPSG:4326','EPSG:3857');
  const maxXY= ol.proj.transform([maxLon,maxLat],'EPSG:4326','EPSG:3857');
  const xSpan= Math.abs(maxXY[0] - minXY[0]);
  const ySpan= Math.abs(maxXY[1] - minXY[1]);

  let [width,height]= pickPixelDimensions(xSpan,ySpan);

  showSpinner();
  alert(`Downloading tile ${tileIndex}/${totalTiles} in EPSG:4326...`);

  const wcsUrl= `https://elevation.nationalmap.gov/arcgis/services/3DEPElevation/ImageServer/WCSServer`
    + `?service=WCS&version=1.0.0&request=GetCoverage`
    + `&coverage=DEP3Elevation`
    + `&CRS=EPSG:4326`
    + `&BBOX=${minLon},${minLat},${maxLon},${maxLat}`
    + `&width=${width}&height=${height}`
    + `&format=GeoTIFF`;

  try{
    const resp= await fetch(wcsUrl);
    if(!resp.ok){
      hideSpinner();
      alert(`Tile request failed: ${resp.statusText}`);
      return;
    }
    const arrayBuf= await resp.arrayBuffer();

    hideSpinner();
    const blob= new Blob([arrayBuf],{type:'image/tiff'});
    const urlObj= URL.createObjectURL(blob);

    const link= document.createElement('a');
    link.href= urlObj;
    link.download= `lidar_tile_${tileIndex}_4326.tif`;
    link.click();

    URL.revokeObjectURL(urlObj);
  }catch(err){
    hideSpinner();
    console.error(err);
    alert('Error fetching tile: '+ err);
  }
}

////////////////////////////
// EPSG:3857 DOWNLOAD
////////////////////////////
document.getElementById('download3857Button').addEventListener('click', async()=>{
  if(drawSource.getFeatures().length===0){
    alert('No AOI polygon drawn!');
    return;
  }
  const feature= drawSource.getFeatures()[0];
  const geom= feature.getGeometry();
  if(geom.getType()!=='Polygon'){
    alert('AOI must be a polygon.');
    return;
  }

  const extent3857= geom.getExtent();
  const widthM= Math.abs(extent3857[2]-extent3857[0]);
  const heightM= Math.abs(extent3857[3]-extent3857[1]);
  const areaM2= widthM*heightM;

  if(areaM2 <= MAX_AREA_SQKM*1e6){
    await downloadSingleCoverage3857(geom);
  } else {
    await downloadTiledCoverage3857(geom, areaM2);
  }
});

async function downloadSingleCoverage3857(geom){
  const ext3857= geom.getExtent();
  const [xMin,yMin,xMax,yMax]= ext3857;

  const xSpan= Math.abs(xMax-xMin);
  const ySpan= Math.abs(yMax-yMin);

  let [width,height]= pickPixelDimensions(xSpan,ySpan);

  showSpinner();

  const wcsUrl= `https://elevation.nationalmap.gov/arcgis/services/3DEPElevation/ImageServer/WCSServer`
    + `?service=WCS&version=1.0.0&request=GetCoverage`
    + `&coverage=DEP3Elevation`
    + `&CRS=EPSG:3857`
    + `&BBOX=${xMin},${yMin},${xMax},${yMax}`
    + `&width=${width}&height=${height}`
    + `&format=GeoTIFF`;

  try{
    const resp= await fetch(wcsUrl);
    if(!resp.ok){
      hideSpinner();
      alert(`WCS request failed: ${resp.statusText}`);
      return;
    }
    const arrayBuf= await resp.arrayBuffer();

    hideSpinner();

    const blob= new Blob([arrayBuf],{type:'image/tiff'});
    const urlObj= URL.createObjectURL(blob);

    const link= document.createElement('a');
    link.href= urlObj;
    link.download= 'lidar_usa_3857.tif';
    link.click();

    URL.revokeObjectURL(urlObj);
    alert('LiDAR (EPSG:3857) downloaded!');
  }catch(err){
    hideSpinner();
    console.error(err);
    alert('Error: '+ err);
  }
}

async function downloadTiledCoverage3857(geom, areaM2){
  alert(`AOI ~${(areaM2/1e6).toFixed(2)} km² => tiling into ~200 km² blocks.`);

  const ext3857= geom.getExtent();
  const [xMin,yMin,xMax,yMax]= ext3857;

  const totalW= Math.abs(xMax-xMin);
  const totalH= Math.abs(yMax-yMin);

  const tileCountX= Math.ceil(totalW / TILE_SIZE_METERS);
  const tileCountY= Math.ceil(totalH / TILE_SIZE_METERS);
  const totalTiles= tileCountX*tileCountY;

  alert(`Tiling into ${tileCountX} × ${tileCountY} = ${totalTiles} tiles...`);

  let tileIndex=0;
  for(let row=0; row<tileCountY; row++){
    for(let col=0; col<tileCountX; col++){
      tileIndex++;
      const subXMin= xMin + col*TILE_SIZE_METERS;
      const subXMax= Math.min(subXMin+TILE_SIZE_METERS, xMax);
      const subYMin= yMin + row*TILE_SIZE_METERS;
      const subYMax= Math.min(subYMin+TILE_SIZE_METERS, yMax);

      await downloadOneTile3857([subXMin,subYMin,subXMax,subYMax], tileIndex, totalTiles);
      await new Promise(r=>setTimeout(r,500));
    }
  }
  alert('All tiles downloaded (EPSG:3857).');
}

async function downloadOneTile3857(bbox3857, tileIndex, totalTiles){
  const [xMin,yMin,xMax,yMax]= bbox3857;

  const xSpan= Math.abs(xMax-xMin);
  const ySpan= Math.abs(yMax-yMin);

  let [width,height]= pickPixelDimensions(xSpan,ySpan);

  showSpinner();
  alert(`Downloading tile ${tileIndex}/${totalTiles} (EPSG:3857)...`);

  const wcsUrl= `https://elevation.nationalmap.gov/arcgis/services/3DEPElevation/ImageServer/WCSServer`
    + `?service=WCS&version=1.0.0&request=GetCoverage`
    + `&coverage=DEP3Elevation`
    + `&CRS=EPSG:3857`
    + `&BBOX=${xMin},${yMin},${xMax},${yMax}`
    + `&width=${width}&height=${height}`
    + `&format=GeoTIFF`;

  try{
    const resp= await fetch(wcsUrl);
    if(!resp.ok){
      hideSpinner();
      alert(`Tile request failed: ${resp.statusText}`);
      return;
    }
    const arrayBuf= await resp.arrayBuffer();

    hideSpinner();
    const blob= new Blob([arrayBuf],{type:'image/tiff'});
    const urlObj= URL.createObjectURL(blob);

    const link= document.createElement('a');
    link.href= urlObj;
    link.download= `lidar_tile_${tileIndex}_3857.tif`;
    link.click();

    URL.revokeObjectURL(urlObj);
  }catch(err){
    hideSpinner();
    console.error(err);
    alert('Error fetching tile: '+ err);
  }
}

//////////////////////////////////////
// Helper function for pixel dimension
//////////////////////////////////////
//  1) If both xSpan,ySpan <2500 => use xSpan,ySpan directly (one pixel per meter).
//  2) Otherwise, clamp the bigger dimension to 2500, scale the other dimension proportionally.
function pickPixelDimensions(xSpan, ySpan){
  let width=0, height=0;
  // if both <2500 => 1m resolution
  if(xSpan<=2500 && ySpan<=2500){
    width= Math.max(1, Math.round(xSpan));
    height= Math.max(1, Math.round(ySpan));
  } else {
    // clamp bigger dimension to 2500
    if(xSpan>ySpan){
      width=2500;
      // scale y
      height=Math.max(1,Math.round((ySpan/xSpan)*2500));
    } else {
      height=2500;
      // scale x
      width=Math.max(1,Math.round((xSpan/ySpan)*2500));
    }
  }
  return [width, height];
}
</script>
</body>
</html>
