Custom Algorithms
Developing custom algorithms and specialized functions.
This page documents the custom algorithms example.
1"""
2Advanced Example 3: Custom Algorithms and Specialized Functions
3===============================================================
4
5This example demonstrates:
6- Creating custom Earth Engine algorithms
7- Implementing specialized mathematical functions
8- Building reusable code libraries
9- Performance optimization techniques
10- Advanced computational methods
11
12Prerequisites:
13- Strong programming background
14- Understanding of Earth Engine architecture
15- Knowledge of mathematical algorithms
16- Experience with optimization techniques
17"""
18
19import ee
20import math
21import numpy as np
22from typing import List, Dict, Tuple, Union, Optional
23
24class CustomAlgorithms:
25 """Library of custom Earth Engine algorithms and specialized functions."""
26
27 def __init__(self, project_id: str):
28 """Initialize custom algorithms library."""
29 self.project_id = project_id
30 self.initialize_ee()
31
32 def initialize_ee(self):
33 """Initialize Earth Engine with error handling."""
34 try:
35 ee.Initialize(project=self.project_id)
36 print("✓ Earth Engine initialized for custom algorithms")
37 except Exception as e:
38 print(f"✗ Failed to initialize Earth Engine: {e}")
39 raise
40
41 def adaptive_threshold_segmentation(self, image: ee.Image,
42 window_size: int = 15,
43 threshold_factor: float = 0.1) -> ee.Image:
44 """
45 Implement adaptive threshold segmentation algorithm.
46
47 Args:
48 image: Input single-band image
49 window_size: Size of local window for threshold calculation
50 threshold_factor: Factor for threshold adjustment
51
52 Returns:
53 ee.Image: Binary segmented image
54 """
55 print(f"🔍 Applying adaptive threshold segmentation (window: {window_size})")
56
57 # Create kernel for local statistics
58 kernel = ee.Kernel.square(radius=window_size//2, units='pixels')
59
60 # Calculate local mean and standard deviation
61 local_mean = image.reduceNeighborhood(
62 reducer=ee.Reducer.mean(),
63 kernel=kernel
64 )
65
66 local_stddev = image.reduceNeighborhood(
67 reducer=ee.Reducer.stdDev(),
68 kernel=kernel
69 )
70
71 # Adaptive threshold = local_mean + threshold_factor * local_stddev
72 adaptive_threshold = local_mean.add(
73 local_stddev.multiply(threshold_factor)
74 )
75
76 # Apply threshold
77 segmented = image.gt(adaptive_threshold)
78
79 return segmented.rename('adaptive_segmentation')
80
81 def texture_analysis_glcm(self, image: ee.Image,
82 window_size: int = 7,
83 angles: List[int] = [0, 45, 90, 135]) -> ee.Image:
84 """
85 Calculate texture features using Gray-Level Co-occurrence Matrix (GLCM).
86
87 Args:
88 image: Input single-band image
89 window_size: Size of analysis window
90 angles: List of angles for GLCM calculation
91
92 Returns:
93 ee.Image: Multi-band image with texture features
94 """
95 print(f"📐 Calculating GLCM texture features (window: {window_size})")
96
97 # Normalize image to 0-255 range for GLCM
98 normalized = image.unitScale(0, 255).uint8()
99
100 texture_bands = []
101
102 for angle in angles:
103 # Calculate GLCM for each angle
104 glcm = normalized.glcmTexture(size=window_size, kernel=None, average=True)
105
106 # Extract specific texture measures
107 contrast = glcm.select(f'.*_contrast')
108 dissimilarity = glcm.select(f'.*_diss')
109 homogeneity = glcm.select(f'.*_idm')
110 energy = glcm.select(f'.*_asm')
111 entropy = glcm.select(f'.*_ent')
112
113 # Rename bands with angle suffix
114 contrast = contrast.rename(f'contrast_{angle}')
115 dissimilarity = dissimilarity.rename(f'dissimilarity_{angle}')
116 homogeneity = homogeneity.rename(f'homogeneity_{angle}')
117 energy = energy.rename(f'energy_{angle}')
118 entropy = entropy.rename(f'entropy_{angle}')
119
120 texture_bands.extend([contrast, dissimilarity, homogeneity, energy, entropy])
121
122 # Combine all texture bands
123 texture_image = ee.Image.cat(texture_bands)
124
125 return texture_image
126
127 def morphological_operations(self, binary_image: ee.Image,
128 operation: str = 'opening',
129 kernel_size: int = 3,
130 iterations: int = 1) -> ee.Image:
131 """
132 Implement morphological operations for binary images.
133
134 Args:
135 binary_image: Input binary image
136 operation: Type of operation ('erosion', 'dilation', 'opening', 'closing')
137 kernel_size: Size of morphological kernel
138 iterations: Number of iterations to apply
139
140 Returns:
141 ee.Image: Processed binary image
142 """
143 print(f"🔧 Applying morphological {operation} (kernel: {kernel_size}, iterations: {iterations})")
144
145 # Create morphological kernel
146 kernel = ee.Kernel.square(radius=kernel_size//2, units='pixels')
147
148 def erosion(img):
149 """Morphological erosion."""
150 return img.reduceNeighborhood(
151 reducer=ee.Reducer.min(),
152 kernel=kernel
153 )
154
155 def dilation(img):
156 """Morphological dilation."""
157 return img.reduceNeighborhood(
158 reducer=ee.Reducer.max(),
159 kernel=kernel
160 )
161
162 def opening(img):
163 """Morphological opening (erosion followed by dilation)."""
164 eroded = erosion(img)
165 return dilation(eroded)
166
167 def closing(img):
168 """Morphological closing (dilation followed by erosion)."""
169 dilated = dilation(img)
170 return erosion(dilated)
171
172 # Define operations
173 operations = {
174 'erosion': erosion,
175 'dilation': dilation,
176 'opening': opening,
177 'closing': closing
178 }
179
180 if operation not in operations:
181 raise ValueError(f"Unknown operation: {operation}")
182
183 # Apply operation for specified iterations
184 result = binary_image
185 for i in range(iterations):
186 result = operations[operation](result)
187
188 return result.rename(f'{operation}_result')
189
190 def spectral_angle_mapper(self, image: ee.Image,
191 reference_spectra: Dict[str, List[float]]) -> ee.Image:
192 """
193 Implement Spectral Angle Mapper (SAM) classification.
194
195 Args:
196 image: Multi-band input image
197 reference_spectra: Dictionary of class names and their spectral signatures
198
199 Returns:
200 ee.Image: Classification result with spectral angle distances
201 """
202 print(f"📊 Applying Spectral Angle Mapper with {len(reference_spectra)} classes")
203
204 # Get band names
205 band_names = image.bandNames()
206 num_bands = band_names.length()
207
208 classification_bands = []
209
210 for class_name, spectrum in reference_spectra.items():
211 # Convert reference spectrum to Earth Engine image
212 ref_spectrum = ee.Image.constant(spectrum).rename(band_names)
213
214 # Calculate spectral angle
215 # SAM = arccos(sum(pixel * reference) / (||pixel|| * ||reference||))
216
217 # Dot product
218 dot_product = image.multiply(ref_spectrum).reduce(ee.Reducer.sum())
219
220 # Magnitudes
221 image_magnitude = image.pow(2).reduce(ee.Reducer.sum()).sqrt()
222 ref_magnitude = ref_spectrum.pow(2).reduce(ee.Reducer.sum()).sqrt()
223
224 # Cosine of angle
225 cos_angle = dot_product.divide(image_magnitude.multiply(ref_magnitude))
226
227 # Spectral angle in radians
228 spectral_angle = cos_angle.acos()
229
230 # Add to classification bands
231 classification_bands.append(spectral_angle.rename(f'sam_{class_name}'))
232
233 # Combine all SAM bands
234 sam_image = ee.Image.cat(classification_bands)
235
236 # Find class with minimum spectral angle
237 class_assignment = sam_image.reduce(ee.Reducer.min(len(reference_spectra) + 1))
238
239 return sam_image.addBands(class_assignment.rename('sam_classification'))
240
241 def principal_component_analysis(self, image: ee.Image,
242 region: ee.Geometry,
243 scale: int = 30) -> Dict[str, ee.Image]:
244 """
245 Implement Principal Component Analysis (PCA) transformation.
246
247 Args:
248 image: Multi-band input image
249 region: Region for calculating statistics
250 scale: Scale for calculations
251
252 Returns:
253 Dict containing PCA components and statistics
254 """
255 print("📈 Performing Principal Component Analysis")
256
257 # Get band names
258 band_names = image.bandNames()
259
260 # Calculate mean values
261 mean_dict = image.reduceRegion(
262 reducer=ee.Reducer.mean(),
263 geometry=region,
264 scale=scale,
265 maxPixels=1e9
266 )
267
268 # Create mean image
269 means = ee.Image.constant(mean_dict.values(band_names))
270
271 # Center the data (subtract means)
272 centered = image.subtract(means)
273
274 # Calculate covariance matrix
275 # This is a simplified implementation - full PCA requires eigenvalue decomposition
276 covariance = centered.toArray().reduce(
277 reducer=ee.Reducer.covariance(),
278 axes=[0, 1]
279 )
280
281 # For demonstration, create simplified PC transformation
282 # In practice, you would need to compute eigenvectors
283
284 # Simplified PC1 (weighted average emphasizing variance)
285 pc1_weights = [0.4, 0.3, 0.2, 0.1] # Example weights
286 pc1 = image.expression(
287 'b1*w1 + b2*w2 + b3*w3 + b4*w4',
288 {
289 'b1': image.select(0),
290 'b2': image.select(1),
291 'b3': image.select(2) if image.bandNames().size().gt(2) else image.select(0),
292 'b4': image.select(3) if image.bandNames().size().gt(3) else image.select(0),
293 'w1': pc1_weights[0],
294 'w2': pc1_weights[1],
295 'w3': pc1_weights[2],
296 'w4': pc1_weights[3]
297 }
298 ).rename('PC1')
299
300 # Simplified PC2 (orthogonal to PC1)
301 pc2_weights = [0.1, -0.2, 0.3, 0.4]
302 pc2 = image.expression(
303 'b1*w1 + b2*w2 + b3*w3 + b4*w4',
304 {
305 'b1': image.select(0),
306 'b2': image.select(1),
307 'b3': image.select(2) if image.bandNames().size().gt(2) else image.select(0),
308 'b4': image.select(3) if image.bandNames().size().gt(3) else image.select(0),
309 'w1': pc2_weights[0],
310 'w2': pc2_weights[1],
311 'w3': pc2_weights[2],
312 'w4': pc2_weights[3]
313 }
314 ).rename('PC2')
315
316 pca_image = ee.Image.cat([pc1, pc2])
317
318 return {
319 'pca_image': pca_image,
320 'mean_values': mean_dict,
321 'centered_image': centered
322 }
323
324 def edge_detection_sobel(self, image: ee.Image) -> ee.Image:
325 """
326 Implement Sobel edge detection algorithm.
327
328 Args:
329 image: Input single-band image
330
331 Returns:
332 ee.Image: Edge magnitude and direction
333 """
334 print("🔍 Applying Sobel edge detection")
335
336 # Sobel kernels
337 sobel_x = ee.Kernel.fixed(3, 3, [
338 [-1, 0, 1],
339 [-2, 0, 2],
340 [-1, 0, 1]
341 ])
342
343 sobel_y = ee.Kernel.fixed(3, 3, [
344 [-1, -2, -1],
345 [ 0, 0, 0],
346 [ 1, 2, 1]
347 ])
348
349 # Apply Sobel filters
350 gradient_x = image.convolve(sobel_x).rename('gradient_x')
351 gradient_y = image.convolve(sobel_y).rename('gradient_y')
352
353 # Calculate magnitude and direction
354 magnitude = gradient_x.pow(2).add(gradient_y.pow(2)).sqrt().rename('edge_magnitude')
355 direction = gradient_y.atan2(gradient_x).rename('edge_direction')
356
357 return ee.Image.cat([gradient_x, gradient_y, magnitude, direction])
358
359 def watershed_segmentation(self, image: ee.Image,
360 markers: ee.Image,
361 connectivity: int = 8) -> ee.Image:
362 """
363 Simplified watershed segmentation implementation.
364
365 Args:
366 image: Input grayscale image
367 markers: Seed points for watershed
368 connectivity: Pixel connectivity (4 or 8)
369
370 Returns:
371 ee.Image: Segmented image
372 """
373 print(f"💧 Applying watershed segmentation (connectivity: {connectivity})")
374
375 # This is a simplified implementation
376 # Full watershed requires iterative region growing
377
378 # Create distance transform from markers
379 distance = markers.distance(ee.Kernel.euclidean(radius=100))
380
381 # Apply watershed-like segmentation using distance and image gradient
382 gradient = self.edge_detection_sobel(image).select('edge_magnitude')
383
384 # Combine distance and gradient information
385 watershed_function = distance.multiply(-1).add(gradient.multiply(0.1))
386
387 # Create regions using connected components
388 # This is a simplified approach - true watershed is more complex
389 regions = watershed_function.connectedComponents(
390 connectedness=ee.Kernel.plus() if connectivity == 4 else ee.Kernel.square(1),
391 maxSize=1000
392 )
393
394 return regions.rename('watershed_segments')
395
396 def fractal_dimension_calculation(self, binary_image: ee.Image,
397 scales: List[int] = [1, 2, 4, 8, 16]) -> ee.Image:
398 """
399 Calculate fractal dimension using box-counting method.
400
401 Args:
402 binary_image: Input binary image
403 scales: List of box sizes for calculation
404
405 Returns:
406 ee.Image: Fractal dimension map
407 """
408 print(f"📐 Calculating fractal dimension with scales: {scales}")
409
410 box_counts = []
411
412 for scale in scales:
413 # Create kernel for box counting
414 kernel = ee.Kernel.square(radius=scale, units='pixels')
415
416 # Count boxes containing features
417 box_count = binary_image.reduceNeighborhood(
418 reducer=ee.Reducer.sum(),
419 kernel=kernel
420 ).gt(0) # Convert to binary (box contains feature or not)
421
422 # Sum boxes in larger neighborhoods to get local fractal properties
423 local_count = box_count.reduceNeighborhood(
424 reducer=ee.Reducer.sum(),
425 kernel=ee.Kernel.square(radius=scale*2, units='pixels')
426 )
427
428 box_counts.append(local_count.rename(f'count_{scale}'))
429
430 # Combine all scales
431 count_image = ee.Image.cat(box_counts)
432
433 # Calculate fractal dimension using linear regression
434 # This is a simplified implementation
435 # True fractal dimension requires log-log regression across scales
436
437 # For demonstration, use ratio of counts at different scales
438 if len(scales) >= 2:
439 fractal_approx = count_image.select(f'count_{scales[0]}').divide(
440 count_image.select(f'count_{scales[-1]}').add(1)
441 ).log().rename('fractal_dimension')
442 else:
443 fractal_approx = count_image.select(0).rename('fractal_dimension')
444
445 return count_image.addBands(fractal_approx)
446
447 def multi_scale_analysis(self, image: ee.Image,
448 scales: List[int] = [1, 2, 4, 8, 16, 32],
449 operations: List[str] = ['mean', 'stdDev', 'variance']) -> ee.Image:
450 """
451 Perform multi-scale analysis using different scales and operations.
452
453 Args:
454 image: Input image
455 scales: List of analysis scales
456 operations: List of operations to perform
457
458 Returns:
459 ee.Image: Multi-scale feature image
460 """
461 print(f"🔍 Multi-scale analysis with {len(scales)} scales and {len(operations)} operations")
462
463 scale_bands = []
464
465 for scale in scales:
466 kernel = ee.Kernel.square(radius=scale, units='pixels')
467
468 for operation in operations:
469 if operation == 'mean':
470 result = image.reduceNeighborhood(
471 reducer=ee.Reducer.mean(),
472 kernel=kernel
473 )
474 elif operation == 'stdDev':
475 result = image.reduceNeighborhood(
476 reducer=ee.Reducer.stdDev(),
477 kernel=kernel
478 )
479 elif operation == 'variance':
480 result = image.reduceNeighborhood(
481 reducer=ee.Reducer.variance(),
482 kernel=kernel
483 )
484 elif operation == 'entropy':
485 # Approximate entropy using histogram
486 result = image.reduceNeighborhood(
487 reducer=ee.Reducer.entropy(),
488 kernel=kernel
489 )
490 else:
491 continue
492
493 scale_bands.append(result.rename(f'{operation}_scale_{scale}'))
494
495 return ee.Image.cat(scale_bands)
496
497 def optimize_algorithm_performance(self) -> Dict[str, str]:
498 """
499 Provide performance optimization strategies for custom algorithms.
500
501 Returns:
502 Dict of optimization strategies
503 """
504 strategies = {
505 'vectorization': 'Use ee.Image operations instead of loops when possible',
506 'kernel_optimization': 'Choose appropriate kernel sizes for operations',
507 'band_selection': 'Process only necessary bands to reduce computation',
508 'scale_matching': 'Match processing scale to data resolution',
509 'chunked_processing': 'Break large operations into smaller chunks',
510 'memory_management': 'Use appropriate data types and avoid large intermediate results',
511 'caching': 'Cache frequently used intermediate results',
512 'parallel_processing': 'Design algorithms to leverage Earth Engine\'s parallelism'
513 }
514
515 print("⚡ Performance Optimization Strategies:")
516 for strategy, description in strategies.items():
517 print(f" • {strategy}: {description}")
518
519 return strategies
520
521def main():
522 """Main function demonstrating custom algorithms."""
523
524 # Initialize custom algorithms library
525 algorithms = CustomAlgorithms('your-project-id')
526
527 print("="*80)
528 print("🧮 CUSTOM ALGORITHMS AND SPECIALIZED FUNCTIONS")
529 print("="*80)
530
531 # Load test image
532 image = (ee.Image('LANDSAT/LC08/C02/T1_L2/LC08_044034_20140318')
533 .select(['SR_B2', 'SR_B3', 'SR_B4', 'SR_B5'])
534 .multiply(0.0000275).add(-0.2))
535
536 test_region = ee.Geometry.Rectangle([-122.5, 37.5, -122.0, 38.0])
537
538 print(f"Using test image: Landsat 8")
539 print(f"Test region: San Francisco Bay Area")
540
541 # Example 1: Adaptive threshold segmentation
542 print("\n1️⃣ Adaptive Threshold Segmentation")
543 print("-" * 40)
544
545 # Use NDVI for segmentation
546 ndvi = image.normalizedDifference(['SR_B5', 'SR_B4'])
547 segmented = algorithms.adaptive_threshold_segmentation(
548 ndvi, window_size=15, threshold_factor=0.1
549 )
550 print("✓ Adaptive segmentation completed")
551
552 # Example 2: Texture analysis
553 print("\n2️⃣ Texture Analysis (GLCM)")
554 print("-" * 30)
555
556 # Use NIR band for texture analysis
557 nir_band = image.select('SR_B5')
558 texture_features = algorithms.texture_analysis_glcm(
559 nir_band, window_size=7, angles=[0, 45, 90, 135]
560 )
561 print(f"✓ Texture analysis completed with {texture_features.bandNames().size().getInfo()} features")
562
563 # Example 3: Morphological operations
564 print("\n3️⃣ Morphological Operations")
565 print("-" * 35)
566
567 # Create binary image from NDVI threshold
568 binary_image = ndvi.gt(0.3)
569
570 # Apply different morphological operations
571 operations = ['erosion', 'dilation', 'opening', 'closing']
572 morph_results = {}
573
574 for operation in operations:
575 result = algorithms.morphological_operations(
576 binary_image, operation=operation, kernel_size=3, iterations=1
577 )
578 morph_results[operation] = result
579 print(f"✓ {operation} completed")
580
581 # Example 4: Spectral Angle Mapper
582 print("\n4️⃣ Spectral Angle Mapper")
583 print("-" * 30)
584
585 # Define reference spectra (example values)
586 reference_spectra = {
587 'water': [0.05, 0.08, 0.06, 0.04],
588 'vegetation': [0.04, 0.08, 0.15, 0.35],
589 'urban': [0.15, 0.18, 0.22, 0.25],
590 'soil': [0.12, 0.15, 0.18, 0.20]
591 }
592
593 sam_result = algorithms.spectral_angle_mapper(image, reference_spectra)
594 print(f"✓ SAM classification with {len(reference_spectra)} classes")
595
596 # Example 5: Principal Component Analysis
597 print("\n5️⃣ Principal Component Analysis")
598 print("-" * 35)
599
600 pca_results = algorithms.principal_component_analysis(
601 image, test_region, scale=30
602 )
603 print("✓ PCA transformation completed")
604 print(f"✓ PC bands: {pca_results['pca_image'].bandNames().getInfo()}")
605
606 # Example 6: Edge detection
607 print("\n6️⃣ Sobel Edge Detection")
608 print("-" * 28)
609
610 # Use first band for edge detection
611 edges = algorithms.edge_detection_sobel(image.select(0))
612 print(f"✓ Edge detection completed with {edges.bandNames().size().getInfo()} output bands")
613
614 # Example 7: Multi-scale analysis
615 print("\n7️⃣ Multi-Scale Analysis")
616 print("-" * 28)
617
618 multiscale_features = algorithms.multi_scale_analysis(
619 ndvi,
620 scales=[1, 2, 4, 8, 16],
621 operations=['mean', 'stdDev', 'variance']
622 )
623 num_features = multiscale_features.bandNames().size().getInfo()
624 print(f"✓ Multi-scale analysis with {num_features} features")
625
626 # Example 8: Fractal dimension
627 print("\n8️⃣ Fractal Dimension Calculation")
628 print("-" * 35)
629
630 fractal_result = algorithms.fractal_dimension_calculation(
631 binary_image, scales=[1, 2, 4, 8, 16]
632 )
633 print("✓ Fractal dimension calculated")
634
635 # Example 9: Performance optimization
636 print("\n9️⃣ Performance Optimization")
637 print("-" * 32)
638
639 optimization_strategies = algorithms.optimize_algorithm_performance()
640
641 # Summary
642 print("\n" + "="*80)
643 print("📊 CUSTOM ALGORITHMS SUMMARY")
644 print("="*80)
645
646 print("\n🎯 Algorithms Demonstrated:")
647 print("• Adaptive threshold segmentation")
648 print("• Gray-Level Co-occurrence Matrix (GLCM) texture analysis")
649 print("• Morphological operations (erosion, dilation, opening, closing)")
650 print("• Spectral Angle Mapper (SAM) classification")
651 print("• Principal Component Analysis (PCA)")
652 print("• Sobel edge detection")
653 print("• Multi-scale spatial analysis")
654 print("• Fractal dimension calculation")
655
656 print("\n📈 Results Generated:")
657 print(f"• Segmented regions: {segmented.bandNames().size().getInfo()} bands")
658 print(f"• Texture features: {texture_features.bandNames().size().getInfo()} bands")
659 print(f"• Morphological results: {len(morph_results)} operations")
660 print(f"• SAM classification: {sam_result.bandNames().size().getInfo()} bands")
661 print(f"• Edge features: {edges.bandNames().size().getInfo()} bands")
662 print(f"• Multi-scale features: {num_features} bands")
663
664 print("\n🏆 Advanced Techniques Applied:")
665 print("• Custom mathematical implementations")
666 print("• Neighborhood analysis and convolution")
667 print("• Multi-dimensional data processing")
668 print("• Performance optimization strategies")
669 print("• Object-oriented algorithm design")
670 print("• Reusable function libraries")
671
672 print("\n✅ Custom Algorithms Development Complete!")
673
674if __name__ == "__main__":
675 main()