diff --git a/api_server.py b/api_server.py
index 7fb1b91..756c6e1 100644
--- a/api_server.py
+++ b/api_server.py
@@ -224,6 +224,9 @@ class PredictionConfig(BaseModel):
     # Cloud removal strategy
     cloud_removal_method: str = "classic"
     cloud_removal_model: Optional[str] = None  # Optional: .pth model filename for deep learning cloud removal
+    
+    # Shapefile overlay for visualization
+    shapefile_overlay: Optional[str] = None  # Path to shapefile for overlaying boundaries
 
 
 class TrainingStatus(BaseModel):
@@ -286,6 +289,7 @@ class PredictionWithNDVIConfig(BaseModel):
     export_classification: bool = True  # Export classification raster
     cloud_removal_method: str = "classic"
     cloud_removal_model: Optional[str] = None  # Optional: .pth model filename for deep learning cloud removal
+    shapefile_overlay: Optional[str] = None  # Path to shapefile for overlaying boundaries
 
 
 class CloudRemovalTrainingConfig(BaseModel):
@@ -1228,6 +1232,87 @@ async def list_training_files():
     }
 
 
+@app.get("/api/overlay/shapefiles")
+async def list_overlay_shapefiles():
+    """Liệt kê các shapefile có sẵn cho overlay trên prediction"""
+    overlay_dirs = ["region", "ChauThanh", "ThuanHoa"]
+    shapefiles = []
+    
+    for overlay_dir in overlay_dirs:
+        dir_path = Path(overlay_dir)
+        if not dir_path.exists():
+            continue
+            
+        # Find all .shp files in this directory and subdirectories
+        for shp_file in dir_path.rglob("*.shp"):
+            try:
+                file_size = shp_file.stat().st_size
+                file_modified = datetime.fromtimestamp(shp_file.stat().st_mtime).isoformat()
+                
+                # Try to read shapefile to get feature count and bbox
+                import geopandas as gpd
+                gdf = gpd.read_file(str(shp_file))
+                feature_count = len(gdf)
+                
+                # Calculate bbox (always in EPSG:4326 for consistency)
+                bbox = None
+                if not gdf.empty and gdf.crs:
+                    try:
+                        # Reproject to EPSG:4326 if needed
+                        if gdf.crs != "EPSG:4326":
+                            gdf_4326 = gdf.to_crs("EPSG:4326")
+                        else:
+                            gdf_4326 = gdf
+                        
+                        # Get total bounds [minx, miny, maxx, maxy]
+                        bounds = gdf_4326.total_bounds
+                        if len(bounds) == 4:
+                            bbox = [
+                                float(bounds[0]),  # min_lon
+                                float(bounds[1]),  # min_lat
+                                float(bounds[2]),  # max_lon
+                                float(bounds[3])   # max_lat
+                            ]
+                    except Exception as bbox_error:
+                        print(f"[WARNING] Cannot calculate bbox for {shp_file}: {bbox_error}")
+                
+                # Get relative path from workspace root
+                relative_path = str(shp_file)
+                
+                shapefiles.append({
+                    "filename": shp_file.name,
+                    "path": relative_path,
+                    "directory": overlay_dir,
+                    "size_bytes": file_size,
+                    "size_mb": round(file_size / 1024 / 1024, 2),
+                    "modified": file_modified,
+                    "feature_count": feature_count,
+                    "crs": str(gdf.crs) if gdf.crs else "Unknown",
+                    "bbox": bbox  # [min_lon, min_lat, max_lon, max_lat] in EPSG:4326
+                })
+            except Exception as e:
+                # If cannot read shapefile, just add basic info
+                file_size = shp_file.stat().st_size
+                file_modified = datetime.fromtimestamp(shp_file.stat().st_mtime).isoformat()
+                relative_path = str(shp_file)
+                
+                shapefiles.append({
+                    "filename": shp_file.name,
+                    "path": relative_path,
+                    "directory": overlay_dir,
+                    "size_bytes": file_size,
+                    "size_mb": round(file_size / 1024 / 1024, 2),
+                    "modified": file_modified,
+                    "error": f"Cannot read shapefile: {str(e)}"
+                })
+    
+    return {
+        "shapefiles": shapefiles,
+        "count": len(shapefiles),
+        "directories": overlay_dirs
+    }
+
+
 @app.get("/api/training/shapefile/{filename}/labels")
 async def get_shapefile_labels(filename: str):
     """Lấy các label từ một shapefile cụ thể"""
@@ -1729,6 +1814,85 @@ def update_progress(message: str):
     print(f"[PROGRESS] {message}")
 
 
+def rasterize_shapefile_overlay(shapefile_path, reference_raster, boundary_value=255):
+    """
+    Rasterize shapefile boundaries to overlay on prediction result.
+    
+    Args:
+        shapefile_path: Path to shapefile
+        reference_raster: xarray DataArray to match dimensions and CRS
+        boundary_value: Value to use for boundaries (default 255 for white)
+    
+    Returns:
+        numpy array with boundaries, same shape as reference_raster
+    """
+    try:
+        import geopandas as gpd
+        from rasterio.features import rasterize
+        import numpy as np
+        
+        # Read shapefile
+        gdf = gpd.read_file(shapefile_path)
+        print(f"[OVERLAY] Loaded shapefile with {len(gdf)} features, CRS: {gdf.crs}")
+        
+        # Ensure CRS matches
+        target_crs = reference_raster.rio.crs
+        if gdf.crs != target_crs:
+            print(f"[OVERLAY] Reprojecting from {gdf.crs} to {target_crs}")
+            gdf = gdf.to_crs(target_crs)
+        
+        # Get raster dimensions and transform
+        height, width = reference_raster.shape
+        transform = reference_raster.rio.transform()
+        print(f"[OVERLAY] Raster dimensions: {height}x{width}")
+        print(f"[OVERLAY] Transform: {transform}")
+        
+        # Calculate appropriate buffer size based on pixel resolution
+        # Get pixel size from transform (transform[0] is x resolution)
+        pixel_size = abs(transform[0])  # in CRS units
+        # Very thin boundary - only 0.2 pixels wide for 1px line
+        buffer_distance = pixel_size * 0.2
+        print(f"[OVERLAY] Pixel size: {pixel_size}, Buffer distance: {buffer_distance} (thin 1px line)")
+        
+        # Create boundary geometries with minimal buffering
+        boundary_geoms = []
+        for idx, geom in enumerate(gdf.geometry):
+            if geom is not None and geom.is_valid:
+                # Get boundary of each polygon
+                boundary = geom.boundary
+                if boundary is not None:
+                    # Minimal buffer for 1-pixel thin line
+                    buffered = boundary.buffer(buffer_distance)
+                    boundary_geoms.append((buffered, boundary_value))
+        
+        if not boundary_geoms:
+            print(f"[WARNING] No valid boundary geometries found in {shapefile_path}")
+            return np.zeros((height, width), dtype=np.uint8)
+        
+        print(f"[OVERLAY] Rasterizing {len(boundary_geoms)} boundaries...")
+        
+        # Rasterize boundaries
+        boundary_mask = rasterize(
+            shapes=boundary_geoms,
+            out_shape=(height, width),
+            transform=transform,
+            fill=0,  # Background
+            dtype=np.uint8
+        )
+        
+        boundary_count = np.count_nonzero(boundary_mask)
+        print(f"[OVERLAY] Boundary pixels: {boundary_count} / {height*width} ({boundary_count/(height*width)*100:.2f}%)")
+        
+        if boundary_count == 0:
+            print(f"[OVERLAY WARNING] No boundary pixels were rasterized! Check CRS and geometry overlap.")
+        
+        return boundary_mask
+        
+    except Exception as e:
+        print(f"[ERROR] Failed to rasterize shapefile {shapefile_path}: {e}")
+        return None
+
+
 def update_prediction_progress(message: str):
     """Cập nhật prediction progress message"""
     global prediction_status
@@ -2083,6 +2247,33 @@ async def run_prediction(config: PredictionConfig):
         
         prediction_da.rio.to_raster(str(output_file), driver="GTiff")
         
+        # ============ SHAPEFILE OVERLAY ============
+        overlay_mask = None
+        if config.shapefile_overlay:
+            prediction_status["progress"] = f"Đang overlay shapefile: {config.shapefile_overlay}..."
+            print(f"[OVERLAY] Shapefile overlay requested: {config.shapefile_overlay}")
+            
+            # Validate shapefile path exists
+            shapefile_path = Path(config.shapefile_overlay)
+            if not shapefile_path.exists():
+                print(f"[OVERLAY ERROR] Shapefile not found: {shapefile_path}")
+                print(f"[OVERLAY ERROR] Absolute path: {shapefile_path.absolute()}")
+                print(f"[OVERLAY ERROR] Current working directory: {Path.cwd()}")
+            else:
+                print(f"[OVERLAY] Shapefile exists: {shapefile_path.absolute()}")
+                try:
+                    overlay_mask = rasterize_shapefile_overlay(str(shapefile_path), prediction_da)
+                    if overlay_mask is not None and np.count_nonzero(overlay_mask) > 0:
+                        print(f"[OVERLAY] Successfully rasterized shapefile boundaries ({np.count_nonzero(overlay_mask)} pixels)")
+                    else:
+                        print(f"[OVERLAY WARNING] Shapefile rasterized but no boundary pixels found")
+                except Exception as overlay_error:
+                    print(f"[OVERLAY ERROR] Exception: {overlay_error}")
+                    import traceback
+                    traceback.print_exc()
+        else:
+            print(f"[OVERLAY] No shapefile overlay requested")
+        
         # Generate PNG preview for web display
         prediction_status["progress"] = "Đang tạo PNG preview..."
         png_file = output_dir / f"prediction_{timestamp}.png"
@@ -2097,6 +2288,42 @@ async def run_prediction(config: PredictionConfig):
             
             # Plot prediction with colormap
             im = ax.imshow(predictions_2d, cmap='tab20', interpolation='nearest')
+            
+            # Overlay shapefile boundaries if available
+            if overlay_mask is not None and np.count_nonzero(overlay_mask) > 0:
+                print(f"[PNG OVERLAY] Overlaying {np.count_nonzero(overlay_mask)} boundary pixels")
+                
+                # Create a mask for boundaries (where overlay_mask > 0)
+                boundary_mask = overlay_mask > 0
+                
+                # Method: Direct overlay with high-contrast colors
+                # Create RGBA overlay image
+                overlay_rgba = np.zeros((*predictions_2d.shape, 4))
+                overlay_rgba[boundary_mask, 0] = 1.0  # Red = 1.0 (white)
+                overlay_rgba[boundary_mask, 1] = 1.0  # Green = 1.0 (white)
+                overlay_rgba[boundary_mask, 2] = 1.0  # Blue = 1.0 (white)
+                overlay_rgba[boundary_mask, 3] = 1.0  # Alpha = 1.0 (fully opaque)
+                
+                # Overlay on top of prediction
+                ax.imshow(overlay_rgba, interpolation='nearest')
+                
+                # Also add a black outline for better contrast
+                from scipy import ndimage
+                boundary_dilated = ndimage.binary_dilation(boundary_mask, iterations=1)
+                boundary_outline = boundary_dilated & ~boundary_mask
+                
+                outline_rgba = np.zeros((*predictions_2d.shape, 4))
+                outline_rgba[boundary_outline, 0] = 0.0  # Black outline
+                outline_rgba[boundary_outline, 1] = 0.0
+                outline_rgba[boundary_outline, 2] = 0.0
+                outline_rgba[boundary_outline, 3] = 0.8
+                
+                ax.imshow(outline_rgba, interpolation='nearest')
+                
+                print(f"[PNG OVERLAY] Added shapefile boundaries to visualization (direct overlay method)")
+            else:
+                print(f"[PNG OVERLAY] No overlay mask or empty mask (pixels: {np.count_nonzero(overlay_mask) if overlay_mask is not None else 0})")
+            
             ax.set_title(f'Land Classification - {timestamp}', fontsize=16, fontweight='bold', pad=20)
             ax.set_xlabel('X (pixels)', fontsize=11)
             ax.set_ylabel('Y (pixels)', fontsize=11)
@@ -2176,7 +2403,9 @@ async def run_prediction(config: PredictionConfig):
             "n_features": features.shape[1],
             "feature_mode": feature_mode,
             "used_radar": use_radar,
-            "model_used": config.model_filename
+            "model_used": config.model_filename,
+            "shapefile_overlay": config.shapefile_overlay,
+            "overlay_applied": overlay_mask is not None
         }
         
         # Auto generate prediction report
@@ -4355,8 +4584,69 @@ async def predict_with_ndvi(config: PredictionWithNDVIConfig, background_tasks:
                 import matplotlib.pyplot as plt
                 from matplotlib.patches import Patch
                 
+                # ============ SHAPEFILE OVERLAY ============
+                overlay_mask = None
+                if config.shapefile_overlay:
+                    print(f"[OVERLAY] Shapefile overlay requested: {config.shapefile_overlay}")
+                    
+                    # Validate shapefile path exists
+                    shapefile_path = Path(config.shapefile_overlay)
+                    if not shapefile_path.exists():
+                        print(f"[OVERLAY ERROR] Shapefile not found: {shapefile_path}")
+                        print(f"[OVERLAY ERROR] Absolute path: {shapefile_path.absolute()}")
+                    else:
+                        print(f"[OVERLAY] Shapefile exists: {shapefile_path.absolute()}")
+                        try:
+                            # Import rioxarray for rio accessor
+                            import rioxarray
+                            
+                            # Create temporary DataArray for rasterization
+                            temp_da = xr.DataArray(
+                                prediction_raster,
+                                coords={
+                                    "y": np.linspace(bbox[3], bbox[1], height),
+                                    "x": np.linspace(bbox[0], bbox[2], width)
+                                },
+                                dims=["y", "x"]
+                            )
+                            temp_da.rio.write_crs("EPSG:4326", inplace=True)
+                            temp_da.rio.write_transform(transform, inplace=True)
+                            
+                            overlay_mask = rasterize_shapefile_overlay(str(shapefile_path), temp_da)
+                            if overlay_mask is not None and np.count_nonzero(overlay_mask) > 0:
+                                print(f"[OVERLAY] Successfully rasterized shapefile boundaries ({np.count_nonzero(overlay_mask)} pixels)")
+                            else:
+                                print(f"[OVERLAY WARNING] Shapefile rasterized but no boundary pixels found")
+                        except Exception as overlay_error:
+                            print(f"[OVERLAY ERROR] Exception: {overlay_error}")
+                            import traceback
+                            traceback.print_exc()
+                else:
+                    print(f"[OVERLAY] No shapefile overlay requested")
+                
                 fig, ax = plt.subplots(figsize=(14, 10), dpi=150)
                 im = ax.imshow(prediction_raster, cmap='tab20', interpolation='nearest')
+                
+                # Overlay shapefile boundaries if available
+                if overlay_mask is not None and np.count_nonzero(overlay_mask) > 0:
+                    print(f"[PNG OVERLAY] Overlaying {np.count_nonzero(overlay_mask)} boundary pixels (1px thin line)")
+                    
+                    # Create a mask for boundaries
+                    boundary_mask = overlay_mask > 0
+                    
+                    # Create RGBA overlay image - thin 1px white line only
+                    overlay_rgba = np.zeros((*prediction_raster.shape, 4))
+                    overlay_rgba[boundary_mask, 0] = 1.0  # White (R=1)
+                    overlay_rgba[boundary_mask, 1] = 1.0  # White (G=1)
+                    overlay_rgba[boundary_mask, 2] = 1.0  # White (B=1)
+                    overlay_rgba[boundary_mask, 3] = 1.0  # Fully opaque
+                    
+                    ax.imshow(overlay_rgba, interpolation='nearest')
+                    
+                    print(f"[PNG OVERLAY] Added thin 1px shapefile boundaries to visualization")
+                else:
+                    print(f"[PNG OVERLAY] No overlay mask or empty mask")
+                
                 ax.set_title(f'Land Classification - {timestamp}', fontsize=16, fontweight='bold', pad=20)
                 ax.set_xlabel('X (pixels)', fontsize=11)
                 ax.set_ylabel('Y (pixels)', fontsize=11)
diff --git a/prediction_interface.html b/prediction_interface.html
index 3274b8b..8ec5467 100644
--- a/prediction_interface.html
+++ b/prediction_interface.html
@@ -774,6 +774,23 @@
                             </div>
                         </label>
                     </div>
+                    
+                    <!-- Shapefile Overlay Option -->
+                    <div class="form-group" style="margin-bottom: 15px;">
+                        <label style="font-weight: 600; color: #c2410c; margin-bottom: 8px; display: block;">
+                            🗺️ Overlay Shapefile (Hiển thị ranh giới lô đất)
+                        </label>
+                        <select id="shapefileOverlay" onchange="onShapefileSelected(event)" style="width: 100%; padding: 10px; border: 2px solid #fdba74; border-radius: 8px; font-size: 14px; background: white;">
+                            <option value="">-- Không overlay --</option>
+                            <!-- Shapefiles will be loaded here -->
+                        </select>
+                        <div style="font-size: 0.85em; color: #9a3412; margin-top: 4px; line-height: 1.4;">
+                            <b>🎯 Tự động cập nhật vùng prediction:</b><br>
+                            ✅ Khi chọn shapefile → <b>Bbox trên bản đồ tự động thay đổi</b> theo vùng shapefile<br>
+                            ✅ <b>CRS sẽ tự động chuyển đổi</b> - không cần lo về EPSG:4326/9209/32648<br>
+                            💡 Không chọn shapefile → Dùng bbox tùy chỉnh do bạn vẽ trên bản đồ
+                        </div>
+                    </div>
                 </div>
                 
                 <button class="btn btn-primary" onclick="startPrediction()" id="predictBtn" style="margin-top: 5px; width: 100%; font-size: 1.1em; padding: 16px;">
@@ -1651,6 +1668,29 @@
                 }
             }
             
+            // Get shapefile overlay option
+            const shapefileOverlay = document.getElementById('shapefileOverlay').value;
+            
+            // Warn if no shapefile selected (optional but recommended)
+            if (!shapefileOverlay) {
+                const confirmWithoutShapefile = confirm(
+                    '⚠️ CẢNH BÁO: Bạn chưa chọn shapefile!\n\n' +
+                    '❌ Kết quả sẽ KHÔNG có đường ranh giới lô đất.\n\n' +
+                    '💡 Để có đường phân lô trên ảnh kết quả:\n' +
+                    '   - Hủy bỏ\n' +
+                    '   - Chọn shapefile trong dropdown "Overlay Shapefile"\n' +
+                    '   - Chạy lại prediction\n\n' +
+                    'Bạn có muốn tiếp tục KHÔNG CÓ ranh giới lô đất không?'
+                );
+                
+                if (!confirmWithoutShapefile) {
+                    console.log('[PREDICTION] User cancelled to select shapefile');
+                    return;
+                }
+            } else {
+                console.log(`[PREDICTION] Shapefile selected: ${shapefileOverlay}`);
+            }
+            
             const config = {
                 model_filename: modelFilename,
                 min_lon: selectedBbox.min_lon,
@@ -1666,7 +1706,8 @@
                 export_ndvi: exportNDVI,
                 export_classification: true,
                 cloud_removal_method: cloudRemovalConfig.method,
-                cloud_removal_model: cloudRemovalConfig.model_filename || null
+                cloud_removal_model: cloudRemovalConfig.model_filename || null,
+                shapefile_overlay: shapefileOverlay || null
             };
             
             try {
@@ -2598,6 +2639,7 @@
             loadPredProvinces();  // Load provinces list
             loadNDVIProvinces();  // Load NDVI provinces list
             loadNDVIModels();  // Load models for NDVI
+            loadOverlayShapefiles();  // Load shapefiles for overlay
             
             // Add event listener for model selection
             document.getElementById('modelSelect').addEventListener('change', updateModelInfo);
@@ -2993,6 +3035,153 @@
             
             alert(`✅ Đã áp dụng preset: ${config.name}\n\nBbox: [${config.bbox.join(', ')}]\nThời gian: ${config.start_date} → ${config.end_date}\nSample points: ${config.sample_points}`);
         }
+        
+        // Load overlay shapefiles
+        async function loadOverlayShapefiles() {
+            try {
+                const response = await fetch('/api/overlay/shapefiles');
+                const data = await response.json();
+                
+                const select = document.getElementById('shapefileOverlay');
+                select.innerHTML = '<option value="">-- Không overlay --</option>';
+                
+                if (data.shapefiles && data.shapefiles.length > 0) {
+                    data.shapefiles.forEach(shp => {
+                        const option = document.createElement('option');
+                        option.value = shp.path;
+                        
+                        // Build detailed label with CRS and bbox info
+                        let label = `${shp.filename} - ${shp.feature_count} features`;
+                        
+                        // Add CRS info (important for matching!)
+                        if (shp.crs) {
+                            const crsCode = shp.crs.split(':').pop(); // Extract code from "EPSG:4326"
+                            label += ` | CRS: ${crsCode}`;
+                        }
+                        
+                        // Add bbox info for easy matching with prediction area
+                        if (shp.bbox && shp.bbox.length === 4) {
+                            const [minLon, minLat, maxLon, maxLat] = shp.bbox;
+                            label += ` | Vùng: [${minLon.toFixed(2)}, ${minLat.toFixed(2)}, ${maxLon.toFixed(2)}, ${maxLat.toFixed(2)}]`;
+                        }
+                        
+                        option.textContent = label;
+                        
+                        // Store full shapefile info as data attributes for later use
+                        option.dataset.crs = shp.crs || '';
+                        option.dataset.bbox = JSON.stringify(shp.bbox || []);
+                        option.dataset.featureCount = shp.feature_count;
+                        
+                        select.appendChild(option);
+                    });
+                    
+                    console.log(`[Overlay Shapefiles] Loaded ${data.shapefiles.length} shapefiles`);
+                } else {
+                    console.log('[Overlay Shapefiles] No shapefiles found');
+                }
+                
+                // Add event listener for shapefile selection change (OUTSIDE the if block)
+                // Remove old listener first to prevent duplicates
+                select.removeEventListener('change', onShapefileSelected);
+                select.addEventListener('change', onShapefileSelected);
+                console.log('[Overlay Shapefiles] Event listener attached');
+                
+            } catch (error) {
+                console.error('[Overlay Shapefiles] Error loading shapefiles:', error);
+                const select = document.getElementById('shapefileOverlay');
+                select.innerHTML = '<option value="">Error loading shapefiles</option>';
+            }
+        }
+        
+        // Handle shapefile selection - auto update bbox on map
+        function onShapefileSelected(event) {
+            console.log('[Shapefile Select] Event triggered');
+            console.log('[Shapefile Select] map exists:', typeof map !== 'undefined');
+            console.log('[Shapefile Select] drawnItems exists:', typeof drawnItems !== 'undefined');
+            
+            const selectedOption = event.target.selectedOptions[0];
+            
+            // If no shapefile selected (empty value), keep current bbox
+            if (!selectedOption || !selectedOption.value) {
+                console.log('[Shapefile Select] No shapefile selected, keeping current bbox');
+                return;
+            }
+            
+            console.log('[Shapefile Select] Selected shapefile:', selectedOption.value);
+            
+            // Get bbox from data attribute
+            const bboxData = selectedOption.dataset.bbox;
+            console.log('[Shapefile Select] Bbox data:', bboxData);
+            
+            if (!bboxData || bboxData === '[]') {
+                console.warn('[Shapefile Select] Selected shapefile has no bbox data');
+                return;
+            }
+            
+            try {
+                const bbox = JSON.parse(bboxData);
+                console.log('[Shapefile Select] Parsed bbox:', bbox);
+                
+                if (bbox.length !== 4) {
+                    console.warn('[Shapefile Select] Invalid bbox format:', bbox);
+                    return;
+                }
+                
+                const [minLon, minLat, maxLon, maxLat] = bbox;
+                
+                // Validate bbox
+                if (minLon < -180 || maxLon > 180 || minLat < -90 || maxLat > 90) {
+                    alert('❌ Bbox của shapefile không hợp lệ!');
+                    return;
+                }
+                
+                console.log('[Shapefile Select] Creating rectangle with bounds:', [[minLat, minLon], [maxLat, maxLon]]);
+                
+                // Update prediction map bbox
+                const bounds = [
+                    [minLat, minLon],
+                    [maxLat, maxLon]
+                ];
+                
+                const rectangle = L.rectangle(bounds, {
+                    color: '#667eea',
+                    weight: 3,
+                    fillOpacity: 0.2
+                });
+                
+                // Clear old bbox and add new one
+                console.log('[Shapefile Select] Clearing old layers...');
+                drawnItems.clearLayers();
+                console.log('[Shapefile Select] Adding new rectangle...');
+                drawnItems.addLayer(rectangle);
+                console.log('[Shapefile Select] Fitting map to bounds...');
+                map.fitBounds(bounds, { padding: [50, 50] });
+                
+                // Update selected bbox variable
+                selectedBbox = {
+                    min_lon: minLon,
+                    min_lat: minLat,
+                    max_lon: maxLon,
+                    max_lat: maxLat
+                };
+                
+                // Save to localStorage
+                localStorage.setItem('prediction_bbox', JSON.stringify(selectedBbox));
+                
+                console.log(`[Shapefile Select] Auto-updated bbox from shapefile:`, selectedBbox);
+                
+                // Show notification
+                const crs = selectedOption.dataset.crs || 'Unknown';
+                alert(`✅ Đã tự động cập nhật vùng prediction theo shapefile!\n\n` +
+                      `📍 Bbox: [${minLon.toFixed(4)}, ${minLat.toFixed(4)}, ${maxLon.toFixed(4)}, ${maxLat.toFixed(4)}]\n` +
+                      `🗺️ CRS: ${crs}\n\n` +
+                      `💡 Bạn có thể điều chỉnh lại bằng cách vẽ lại trên bản đồ nếu muốn.`);
+                
+            } catch (e) {
+                console.error('[Shapefile Select] Error parsing bbox:', e);
+                alert(`❌ Lỗi khi xử lý bbox: ${e.message}`);
+            }
+        }
     </script>
 </body>
 </html>
diff --git a/test_shapefile_api.py b/test_shapefile_api.py
new file mode 100644
index 0000000..30089ff
--- /dev/null
+++ b/test_shapefile_api.py
@@ -0,0 +1,48 @@
+#!/usr/bin/env python3
+"""
+Test script to verify shapefile overlay API returns correct bbox data
+"""
+
+import requests
+import json
+
+def test_shapefile_api():
+    """Test /api/overlay/shapefiles endpoint"""
+    print("Testing /api/overlay/shapefiles endpoint...")
+    
+    try:
+        response = requests.get('http://localhost:8000/api/overlay/shapefiles')
+        
+        if response.status_code == 200:
+            data = response.json()
+            print(f"\n✅ API Response successful")
+            print(f"Total shapefiles: {data.get('count', 0)}")
+            
+            if data.get('shapefiles'):
+                print("\n📋 Shapefile details:")
+                for idx, shp in enumerate(data['shapefiles'], 1):
+                    print(f"\n{idx}. {shp.get('filename')}")
+                    print(f"   Path: {shp.get('path')}")
+                    print(f"   CRS: {shp.get('crs')}")
+                    print(f"   Features: {shp.get('feature_count')}")
+                    print(f"   Bbox: {shp.get('bbox')}")
+                    
+                    # Verify bbox format
+                    bbox = shp.get('bbox')
+                    if bbox and len(bbox) == 4:
+                        print(f"   ✅ Bbox format valid: [minLon, minLat, maxLon, maxLat]")
+                    else:
+                        print(f"   ❌ Bbox format invalid or missing!")
+            else:
+                print("\n⚠️ No shapefiles found")
+        else:
+            print(f"\n❌ API returned status code: {response.status_code}")
+            print(f"Response: {response.text}")
+            
+    except requests.exceptions.ConnectionError:
+        print("\n❌ Cannot connect to API server. Is it running on localhost:8000?")
+    except Exception as e:
+        print(f"\n❌ Error: {e}")
+
+if __name__ == "__main__":
+    test_shapefile_api()
diff --git a/test_shapefile_overlay.py b/test_shapefile_overlay.py
new file mode 100644
index 0000000..ee51744
--- /dev/null
+++ b/test_shapefile_overlay.py
@@ -0,0 +1,79 @@
+#!/usr/bin/env python
+"""
+Test script to verify shapefile overlay functionality
+"""
+import geopandas as gpd
+import numpy as np
+from pathlib import Path
+
+# Test shapefile path
+shapefile_path = "ChauThanh/HienTrang/ChauThanh_kiemke.shp"
+
+print("=" * 70)
+print("TESTING SHAPEFILE OVERLAY")
+print("=" * 70)
+
+# Check if file exists
+shp = Path(shapefile_path)
+print(f"\n1. Checking file existence:")
+print(f"   Path: {shp}")
+print(f"   Exists: {shp.exists()}")
+print(f"   Absolute: {shp.absolute()}")
+
+if shp.exists():
+    # Read shapefile
+    print(f"\n2. Reading shapefile...")
+    gdf = gpd.read_file(str(shp))
+    print(f"   Features: {len(gdf)}")
+    print(f"   CRS: {gdf.crs}")
+    print(f"   Bounds: {gdf.total_bounds}")
+    print(f"   Columns: {list(gdf.columns)}")
+    
+    # Check geometries
+    print(f"\n3. Checking geometries...")
+    valid_count = sum(1 for geom in gdf.geometry if geom is not None and geom.is_valid)
+    print(f"   Valid geometries: {valid_count} / {len(gdf)}")
+    
+    # Sample geometry bounds
+    if len(gdf) > 0:
+        sample_geom = gdf.geometry.iloc[0]
+        print(f"   Sample geometry type: {sample_geom.geom_type}")
+        print(f"   Sample geometry bounds: {sample_geom.bounds}")
+    
+    # Test reprojection to EPSG:4326
+    print(f"\n4. Testing reprojection to EPSG:4326...")
+    try:
+        gdf_4326 = gdf.to_crs("EPSG:4326")
+        print(f"   Success!")
+        print(f"   New bounds: {gdf_4326.total_bounds}")
+    except Exception as e:
+        print(f"   ERROR: {e}")
+    
+    # Test boundary extraction
+    print(f"\n5. Testing boundary extraction...")
+    boundaries = []
+    for geom in gdf.geometry:
+        if geom is not None and geom.is_valid:
+            boundary = geom.boundary
+            if boundary is not None:
+                boundaries.append(boundary)
+    print(f"   Extracted boundaries: {len(boundaries)}")
+    
+    # Test buffering
+    print(f"\n6. Testing buffer...")
+    buffer_size = 0.001  # degrees or meters depending on CRS
+    buffered = []
+    for boundary in boundaries[:10]:  # Test first 10
+        try:
+            buf = boundary.buffer(buffer_size)
+            buffered.append(buf)
+        except Exception as e:
+            print(f"   Buffer error: {e}")
+    print(f"   Successfully buffered: {len(buffered)} / 10")
+    
+else:
+    print("   ERROR: Shapefile not found!")
+
+print("\n" + "=" * 70)
+print("TEST COMPLETE")
+print("=" * 70)
diff --git a/test_shapefile_selection.html b/test_shapefile_selection.html
new file mode 100644
index 0000000..1c109b2
--- /dev/null
+++ b/test_shapefile_selection.html
@@ -0,0 +1,177 @@
+<!DOCTYPE html>
+<html lang="vi">
+<head>
+    <meta charset="UTF-8">
+    <title>Test Shapefile Selection</title>
+    <link rel="stylesheet" href="https://unpkg.com/leaflet@1.9.4/dist/leaflet.css" />
+    <style>
+        body { font-family: Arial, sans-serif; padding: 20px; }
+        #map { height: 400px; border: 2px solid #ccc; margin: 20px 0; }
+        .info-box { background: #f0f0f0; padding: 15px; margin: 10px 0; border-radius: 5px; }
+    </style>
+</head>
+<body>
+    <h1>🧪 Test Shapefile Auto-Select Bbox</h1>
+    
+    <div class="info-box">
+        <h3>Chọn Shapefile:</h3>
+        <select id="shapefileOverlay" onchange="onShapefileSelected(event)" style="width: 100%; padding: 10px; font-size: 14px;">
+            <option value="">-- Chọn shapefile --</option>
+        </select>
+    </div>
+    
+    <div id="map"></div>
+    
+    <div class="info-box">
+        <h3>Current Bbox:</h3>
+        <pre id="bboxInfo">Chưa chọn shapefile</pre>
+    </div>
+    
+    <div class="info-box">
+        <h3>Console Logs:</h3>
+        <pre id="console" style="max-height: 200px; overflow-y: auto; background: #000; color: #0f0; padding: 10px;"></pre>
+    </div>
+    
+    <script src="https://unpkg.com/leaflet@1.9.4/dist/leaflet.js"></script>
+    <script>
+        // Global variables
+        let map, drawnItems, selectedBbox = null;
+        
+        // Custom console.log to display in page
+        const originalLog = console.log;
+        console.log = function(...args) {
+            originalLog.apply(console, args);
+            const consoleEl = document.getElementById('console');
+            consoleEl.textContent += args.join(' ') + '\n';
+            consoleEl.scrollTop = consoleEl.scrollHeight;
+        };
+        
+        // Initialize map
+        function initMap() {
+            map = L.map('map').setView([10.0, 105.8], 10);
+            L.tileLayer('https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', {
+                attribution: '© OpenStreetMap contributors'
+            }).addTo(map);
+            
+            drawnItems = new L.FeatureGroup();
+            map.addLayer(drawnItems);
+            
+            console.log('✅ Map initialized');
+        }
+        
+        // Load shapefiles from API
+        async function loadShapefiles() {
+            try {
+                console.log('📡 Fetching shapefiles from API...');
+                const response = await fetch('http://localhost:8000/api/overlay/shapefiles');
+                const data = await response.json();
+                
+                const select = document.getElementById('shapefileOverlay');
+                select.innerHTML = '<option value="">-- Chọn shapefile --</option>';
+                
+                if (data.shapefiles && data.shapefiles.length > 0) {
+                    data.shapefiles.forEach(shp => {
+                        const option = document.createElement('option');
+                        option.value = shp.path;
+                        
+                        let label = `${shp.filename} - ${shp.feature_count} features`;
+                        if (shp.crs) {
+                            const crsCode = shp.crs.split(':').pop();
+                            label += ` | CRS: ${crsCode}`;
+                        }
+                        if (shp.bbox && shp.bbox.length === 4) {
+                            const [minLon, minLat, maxLon, maxLat] = shp.bbox;
+                            label += ` | [${minLon.toFixed(2)}, ${minLat.toFixed(2)}, ${maxLon.toFixed(2)}, ${maxLat.toFixed(2)}]`;
+                        }
+                        
+                        option.textContent = label;
+                        option.dataset.crs = shp.crs || '';
+                        option.dataset.bbox = JSON.stringify(shp.bbox || []);
+                        option.dataset.featureCount = shp.feature_count;
+                        
+                        select.appendChild(option);
+                    });
+                    
+                    console.log(`✅ Loaded ${data.shapefiles.length} shapefiles`);
+                } else {
+                    console.log('⚠️ No shapefiles found');
+                }
+            } catch (error) {
+                console.error('❌ Error loading shapefiles:', error);
+            }
+        }
+        
+        // Handle shapefile selection
+        function onShapefileSelected(event) {
+            console.log('🔔 Shapefile selection changed');
+            
+            const selectedOption = event.target.selectedOptions[0];
+            
+            if (!selectedOption || !selectedOption.value) {
+                console.log('ℹ️ No shapefile selected');
+                document.getElementById('bboxInfo').textContent = 'Chưa chọn shapefile';
+                return;
+            }
+            
+            const bboxData = selectedOption.dataset.bbox;
+            console.log('📦 Bbox data from option:', bboxData);
+            
+            if (!bboxData || bboxData === '[]') {
+                console.log('⚠️ No bbox data in selected option');
+                return;
+            }
+            
+            try {
+                const bbox = JSON.parse(bboxData);
+                console.log('📊 Parsed bbox:', bbox);
+                
+                if (bbox.length !== 4) {
+                    console.log('❌ Invalid bbox length:', bbox.length);
+                    return;
+                }
+                
+                const [minLon, minLat, maxLon, maxLat] = bbox;
+                
+                // Validate bbox
+                if (minLon < -180 || maxLon > 180 || minLat < -90 || maxLat > 90) {
+                    console.log('❌ Bbox out of valid range');
+                    return;
+                }
+                
+                console.log('✅ Valid bbox:', {minLon, minLat, maxLon, maxLat});
+                
+                // Update map
+                const bounds = [[minLat, minLon], [maxLat, maxLon]];
+                const rectangle = L.rectangle(bounds, {
+                    color: '#667eea',
+                    weight: 3,
+                    fillOpacity: 0.2
+                });
+                
+                drawnItems.clearLayers();
+                drawnItems.addLayer(rectangle);
+                map.fitBounds(bounds, { padding: [50, 50] });
+                
+                selectedBbox = {min_lon: minLon, min_lat: minLat, max_lon: maxLon, max_lat: maxLat};
+                
+                console.log('🗺️ Map updated with new bbox');
+                
+                // Update bbox info display
+                document.getElementById('bboxInfo').textContent = JSON.stringify(selectedBbox, null, 2);
+                
+                alert(`✅ Bbox updated!\n\nmin_lon: ${minLon.toFixed(4)}\nmin_lat: ${minLat.toFixed(4)}\nmax_lon: ${maxLon.toFixed(4)}\nmax_lat: ${maxLat.toFixed(4)}`);
+                
+            } catch (e) {
+                console.error('❌ Error:', e);
+            }
+        }
+        
+        // Initialize on load
+        window.onload = function() {
+            console.log('🚀 Page loaded, initializing...');
+            initMap();
+            loadShapefiles();
+        };
+    </script>
+</body>
+</html>