0
# Geometry Operations
1
2
Comprehensive 2D and 3D geometric operations including point sets, regions, polylines, surfaces, and spatial analysis tools for reservoir characterization and visualization.
3
4
## Capabilities
5
6
### Point Set Operations
7
8
Management and analysis of geometric point collections with spatial operations.
9
10
```python { .api }
11
class GeoPointset:
12
"""Set of geometric points with spatial operations."""
13
14
def __init__(self):
15
"""Create empty point set."""
16
17
def add_point(self, x: float, y: float, z: float = 0.0):
18
"""
19
Add point to set.
20
21
Args:
22
x, y (float): 2D coordinates
23
z (float, optional): Z-coordinate for 3D points
24
"""
25
26
def size(self) -> int:
27
"""Get number of points."""
28
29
def get_point(self, index: int) -> tuple:
30
"""
31
Get point by index.
32
33
Args:
34
index (int): Point index
35
36
Returns:
37
tuple: (x, y, z) coordinates
38
"""
39
40
def bounding_box(self) -> tuple:
41
"""
42
Get bounding box of all points.
43
44
Returns:
45
tuple: (xmin, ymin, xmax, ymax, zmin, zmax)
46
"""
47
```
48
49
### Geographic Regions
50
51
Geographic region operations for spatial selection and containment testing.
52
53
```python { .api }
54
class GeoRegion:
55
"""Geographic region operations and spatial queries."""
56
57
def __init__(self):
58
"""Create empty geographic region."""
59
60
def contains_point(self, x: float, y: float) -> bool:
61
"""
62
Check if point is inside region.
63
64
Args:
65
x, y (float): Point coordinates
66
67
Returns:
68
bool: True if point is inside region
69
"""
70
71
def select_inside(self, pointset: GeoPointset) -> list:
72
"""
73
Select points inside region.
74
75
Args:
76
pointset (GeoPointset): Point set to filter
77
78
Returns:
79
list: Indices of points inside region
80
"""
81
82
def select_outside(self, pointset: GeoPointset) -> list:
83
"""Select points outside region."""
84
```
85
86
### Polyline Operations
87
88
C-based high-performance polyline operations for complex geometric calculations.
89
90
```python { .api }
91
class CPolyline:
92
"""C-based polyline operations for high performance."""
93
94
def __init__(self):
95
"""Create empty polyline."""
96
97
def add_point(self, x: float, y: float):
98
"""Add point to polyline."""
99
100
def size(self) -> int:
101
"""Get number of points."""
102
103
def get_point(self, index: int) -> tuple:
104
"""Get point by index."""
105
106
def segment_length(self, segment: int) -> float:
107
"""
108
Get length of specific segment.
109
110
Args:
111
segment (int): Segment index
112
113
Returns:
114
float: Segment length
115
"""
116
117
def get_length(self) -> float:
118
"""Get total polyline length."""
119
120
def close(self):
121
"""Close polyline by connecting last point to first."""
122
123
def extend_to_edge(self, xmin: float, ymin: float, xmax: float, ymax: float):
124
"""Extend polyline to bounding box edges."""
125
```
126
127
### Polyline Collections
128
129
Management of multiple polylines for complex geometric structures.
130
131
```python { .api }
132
class CPolylineCollection:
133
"""Collection of polylines for complex structures."""
134
135
def __init__(self):
136
"""Create empty polyline collection."""
137
138
def add_polyline(self, polyline: CPolyline):
139
"""Add polyline to collection."""
140
141
def size(self) -> int:
142
"""Get number of polylines."""
143
144
def get_polyline(self, index: int) -> CPolyline:
145
"""Get polyline by index."""
146
```
147
148
### Python Polyline Implementation
149
150
Pure Python polyline implementation with additional convenience methods.
151
152
```python { .api }
153
class Polyline:
154
"""Python polyline implementation with convenience methods."""
155
156
def __init__(self):
157
"""Create empty polyline."""
158
159
def add_point(self, x: float, y: float):
160
"""Add point to polyline."""
161
162
def size(self) -> int:
163
"""Get number of points."""
164
165
def get_points(self) -> list:
166
"""
167
Get all points as list.
168
169
Returns:
170
list: List of (x, y) tuples
171
"""
172
173
def unzip(self) -> tuple:
174
"""
175
Separate coordinates into x and y arrays.
176
177
Returns:
178
tuple: (x_array, y_array)
179
"""
180
181
def extend_to_bbox(self, bbox: tuple):
182
"""Extend polyline to bounding box."""
183
```
184
185
### Surface Operations
186
187
3D surface representation and manipulation for reservoir characterization.
188
189
```python { .api }
190
class Surface:
191
"""3D surface operations and analysis."""
192
193
def __init__(self, nx: int, ny: int, xmin: float, ymin: float,
194
dx: float, dy: float):
195
"""
196
Create surface grid.
197
198
Args:
199
nx, ny (int): Grid dimensions
200
xmin, ymin (float): Origin coordinates
201
dx, dy (float): Grid spacing
202
"""
203
204
def get_nx(self) -> int:
205
"""Get X-dimension size."""
206
207
def get_ny(self) -> int:
208
"""Get Y-dimension size."""
209
210
def get_value(self, x: float, y: float) -> float:
211
"""
212
Get surface value at coordinates.
213
214
Args:
215
x, y (float): World coordinates
216
217
Returns:
218
float: Surface value (typically elevation or property)
219
"""
220
221
def set_value(self, i: int, j: int, value: float):
222
"""
223
Set surface value at grid indices.
224
225
Args:
226
i, j (int): Grid indices
227
value (float): Surface value
228
"""
229
230
def add_surface(self, other_surface):
231
"""Add another surface to this one."""
232
233
def copy(self) -> Surface:
234
"""Create copy of surface."""
235
236
def shift(self, offset: float):
237
"""Shift all surface values by offset."""
238
239
def scale(self, factor: float):
240
"""Scale all surface values by factor."""
241
242
def write(self, filename: str):
243
"""Write surface to file."""
244
```
245
246
### Geometric Tools
247
248
Static utility functions for geometric calculations and analysis.
249
250
```python { .api }
251
class GeometryTools:
252
"""Static geometric utility functions."""
253
254
@staticmethod
255
def distance(x1: float, y1: float, x2: float, y2: float) -> float:
256
"""
257
Calculate distance between two points.
258
259
Args:
260
x1, y1 (float): First point coordinates
261
x2, y2 (float): Second point coordinates
262
263
Returns:
264
float: Euclidean distance
265
"""
266
267
@staticmethod
268
def area_polygon(points: list) -> float:
269
"""
270
Calculate polygon area.
271
272
Args:
273
points (list): List of (x, y) coordinate tuples
274
275
Returns:
276
float: Polygon area
277
"""
278
279
@staticmethod
280
def convex_hull(points: list) -> list:
281
"""
282
Calculate convex hull of points.
283
284
Args:
285
points (list): List of (x, y) coordinate tuples
286
287
Returns:
288
list: Convex hull points
289
"""
290
```
291
292
### XYZ File I/O
293
294
Input/output operations for XYZ coordinate data files.
295
296
```python { .api }
297
class XYZIo:
298
"""XYZ file I/O operations."""
299
300
@staticmethod
301
def save(filename: str, points: list):
302
"""
303
Save points to XYZ file.
304
305
Args:
306
filename (str): Output file path
307
points (list): List of (x, y, z) tuples
308
"""
309
310
@staticmethod
311
def load(filename: str) -> list:
312
"""
313
Load points from XYZ file.
314
315
Args:
316
filename (str): Input file path
317
318
Returns:
319
list: List of (x, y, z) tuples
320
"""
321
```
322
323
## Usage Examples
324
325
### Surface Analysis
326
327
```python
328
from resdata.geometry import Surface, GeometryTools
329
import numpy as np
330
import matplotlib.pyplot as plt
331
332
# Create surface representing structure map
333
nx, ny = 50, 40
334
xmin, ymin = 0.0, 0.0
335
dx, dy = 100.0, 100.0 # 100m grid spacing
336
337
surface = Surface(nx, ny, xmin, ymin, dx, dy)
338
339
# Create synthetic structural surface (anticline)
340
for j in range(ny):
341
for i in range(nx):
342
x = xmin + i * dx
343
y = ymin + j * dy
344
345
# Anticline structure
346
center_x, center_y = 2500.0, 2000.0
347
dist_from_center = GeometryTools.distance(x, y, center_x, center_y)
348
elevation = 1000.0 - (dist_from_center / 50.0) ** 2 / 100.0
349
350
surface.set_value(i, j, elevation)
351
352
print(f"Surface dimensions: {surface.get_nx()} x {surface.get_ny()}")
353
354
# Sample surface values
355
sample_points = [(1000, 1500), (2500, 2000), (4000, 3000)]
356
for x, y in sample_points:
357
elevation = surface.get_value(x, y)
358
print(f"Elevation at ({x}, {y}): {elevation:.1f} m")
359
360
# Save surface to file
361
surface.write("structure_map.surf")
362
```
363
364
### Polyline Operations
365
366
```python
367
from resdata.geometry import CPolyline, Polyline, GeometryTools
368
369
# Create well trajectory as polyline
370
trajectory = CPolyline()
371
372
# Add waypoints for deviated well
373
waypoints = [
374
(1000, 2000), # Surface location
375
(1050, 2000), # Start of deviation
376
(1200, 2050), # Mid-point
377
(1500, 2200), # Target location
378
]
379
380
for x, y in waypoints:
381
trajectory.add_point(x, y)
382
383
print(f"Trajectory points: {trajectory.size()}")
384
print(f"Total length: {trajectory.get_length():.1f} m")
385
386
# Calculate segment lengths
387
for i in range(trajectory.size() - 1):
388
seg_length = trajectory.segment_length(i)
389
print(f"Segment {i+1}: {seg_length:.1f} m")
390
391
# Create Python polyline for easier manipulation
392
py_trajectory = Polyline()
393
for x, y in waypoints:
394
py_trajectory.add_point(x, y)
395
396
# Extract coordinates for plotting
397
points = py_trajectory.get_points()
398
x_coords, y_coords = py_trajectory.unzip()
399
400
# Plot trajectory
401
plt.figure(figsize=(10, 8))
402
plt.plot(x_coords, y_coords, 'bo-', linewidth=2, markersize=6)
403
plt.xlabel('X (m)')
404
plt.ylabel('Y (m)')
405
plt.title('Well Trajectory')
406
plt.grid(True)
407
plt.axis('equal')
408
409
# Add labels for waypoints
410
for i, (x, y) in enumerate(points):
411
plt.annotate(f'P{i+1}', (x, y), xytext=(5, 5), textcoords='offset points')
412
413
plt.show()
414
```
415
416
### Point Set Analysis
417
418
```python
419
from resdata.geometry import GeoPointset, GeoRegion, GeometryTools
420
import numpy as np
421
422
# Create point set representing well locations
423
wells = GeoPointset()
424
425
# Add well locations
426
well_locations = [
427
(1000, 1500, -2000), # (x, y, z)
428
(1200, 1800, -2100),
429
(1500, 1600, -1950),
430
(1800, 2000, -2200),
431
(2000, 1900, -2150),
432
]
433
434
for x, y, z in well_locations:
435
wells.add_point(x, y, z)
436
437
print(f"Number of wells: {wells.size()}")
438
439
# Get bounding box
440
bbox = wells.bounding_box()
441
print(f"Bounding box: X[{bbox[0]:.0f}, {bbox[2]:.0f}], "
442
f"Y[{bbox[1]:.0f}, {bbox[3]:.0f}], Z[{bbox[4]:.0f}, {bbox[5]:.0f}]")
443
444
# Analyze well spacing
445
distances = []
446
for i in range(wells.size()):
447
for j in range(i + 1, wells.size()):
448
p1 = wells.get_point(i)
449
p2 = wells.get_point(j)
450
451
# 2D distance (ignore depth)
452
dist = GeometryTools.distance(p1[0], p1[1], p2[0], p2[1])
453
distances.append(dist)
454
455
if distances:
456
print(f"Well spacing statistics:")
457
print(f" Minimum: {min(distances):.0f} m")
458
print(f" Maximum: {max(distances):.0f} m")
459
print(f" Average: {np.mean(distances):.0f} m")
460
```
461
462
### Polygon Area Calculations
463
464
```python
465
from resdata.geometry import GeometryTools
466
import matplotlib.pyplot as plt
467
468
# Define reservoir boundary polygon
469
reservoir_boundary = [
470
(0, 0),
471
(2000, 0),
472
(2500, 1000),
473
(2000, 2000),
474
(1000, 2200),
475
(0, 1500),
476
(0, 0) # Close polygon
477
]
478
479
# Calculate area
480
area = GeometryTools.area_polygon(reservoir_boundary)
481
print(f"Reservoir area: {area/1e6:.2f} km²")
482
483
# Find convex hull
484
convex_hull = GeometryTools.convex_hull(reservoir_boundary)
485
hull_area = GeometryTools.area_polygon(convex_hull)
486
print(f"Convex hull area: {hull_area/1e6:.2f} km²")
487
print(f"Concavity factor: {area/hull_area:.3f}")
488
489
# Plot polygon and convex hull
490
fig, ax = plt.subplots(figsize=(10, 8))
491
492
# Plot original boundary
493
boundary_x = [p[0] for p in reservoir_boundary]
494
boundary_y = [p[1] for p in reservoir_boundary]
495
ax.plot(boundary_x, boundary_y, 'b-', linewidth=2, label='Reservoir Boundary')
496
ax.fill(boundary_x, boundary_y, alpha=0.3, color='blue')
497
498
# Plot convex hull
499
hull_x = [p[0] for p in convex_hull] + [convex_hull[0][0]] # Close hull
500
hull_y = [p[1] for p in convex_hull] + [convex_hull[0][1]]
501
ax.plot(hull_x, hull_y, 'r--', linewidth=2, label='Convex Hull')
502
503
ax.set_xlabel('X (m)')
504
ax.set_ylabel('Y (m)')
505
ax.set_title('Reservoir Boundary Analysis')
506
ax.legend()
507
ax.grid(True)
508
ax.axis('equal')
509
plt.show()
510
```
511
512
### XYZ Data Processing
513
514
```python
515
from resdata.geometry import XYZIo, Surface
516
import numpy as np
517
518
# Create synthetic elevation data
519
points = []
520
for i in range(100):
521
# Random points with elevation trend
522
x = np.random.uniform(0, 5000)
523
y = np.random.uniform(0, 4000)
524
z = 1000 + 0.0002 * x - 0.0001 * y + np.random.normal(0, 10)
525
points.append((x, y, z))
526
527
# Save to XYZ file
528
XYZIo.save("elevation_data.xyz", points)
529
530
# Load data back
531
loaded_points = XYZIo.load("elevation_data.xyz")
532
print(f"Loaded {len(loaded_points)} elevation points")
533
534
# Extract coordinates for analysis
535
x_coords = [p[0] for p in loaded_points]
536
y_coords = [p[1] for p in loaded_points]
537
z_coords = [p[2] for p in loaded_points]
538
539
print(f"Elevation statistics:")
540
print(f" Min: {min(z_coords):.1f} m")
541
print(f" Max: {max(z_coords):.1f} m")
542
print(f" Mean: {np.mean(z_coords):.1f} m")
543
print(f" Std: {np.std(z_coords):.1f} m")
544
545
# Plot elevation map
546
plt.figure(figsize=(12, 8))
547
scatter = plt.scatter(x_coords, y_coords, c=z_coords, cmap='terrain', s=20)
548
plt.colorbar(scatter, label='Elevation (m)')
549
plt.xlabel('X (m)')
550
plt.ylabel('Y (m)')
551
plt.title('Elevation Map from XYZ Data')
552
plt.axis('equal')
553
plt.grid(True, alpha=0.3)
554
plt.show()
555
```
556
557
## Types
558
559
```python { .api }
560
# Coordinate types
561
Point2D = tuple[float, float] # (x, y)
562
Point3D = tuple[float, float, float] # (x, y, z)
563
BoundingBox = tuple[float, float, float, float, float, float] # (xmin, ymin, xmax, ymax, zmin, zmax)
564
565
# Geometric collections
566
PointList = list[Point2D | Point3D]
567
PolygonPoints = list[Point2D]
568
ConvexHull = list[Point2D]
569
570
# Surface grid parameters
571
SurfaceGrid = tuple[int, int, float, float, float, float] # (nx, ny, xmin, ymin, dx, dy)
572
```