Gravimetry and Subsidence

class ResdataGrav:
    """Gravity simulation operations for reservoir monitoring."""
    
    def __init__(self, grid: Grid, porv_kw: ResdataKW):
        """
        Initialize gravity calculator.
        
        Args:
            grid (Grid): Reservoir grid
            porv_kw (ResdataKW): Pore volume keyword
        """
    
    def add_survey_GRAV(self, survey_name: str, x: float, y: float, z: float):
        """
        Add gravity survey point.
        
        Args:
            survey_name (str): Name of the survey
            x, y, z (float): Survey point coordinates
        """
    
    def add_survey_GRAVDIFF(self, survey_name: str, base_survey: str, 
                           x: float, y: float, z: float):
        """
        Add differential gravity survey point.
        
        Args:
            survey_name (str): Name of the survey
            base_survey (str): Reference survey name
            x, y, z (float): Survey point coordinates
        """
    
    def eval_grav(self, time: datetime = None) -> numpy.ndarray:
        """
        Evaluate gravity response.
        
        Args:
            time (datetime, optional): Evaluation time
            
        Returns:
            numpy.ndarray: Gravity values in milligals (mGal)
        """
    
    def load_PORV(self, porv_kw: ResdataKW):
        """Load pore volume data."""
    
    def load_RHO(self, rho_kw: ResdataKW):
        """Load density data."""
    
    def add_std_density(self, oil_density: float, gas_density: float, 
                       water_density: float):
        """
        Add standard fluid densities.
        
        Args:
            oil_density (float): Oil density (kg/m³)
            gas_density (float): Gas density (kg/m³)
            water_density (float): Water density (kg/m³)
        """

class ResdataSubsidence:
    """Subsidence simulation operations for geomechanical analysis."""
    
    def __init__(self, grid: Grid):
        """
        Initialize subsidence calculator.
        
        Args:
            grid (Grid): Reservoir grid
        """
    
    def add_survey_SUBSIDENCE(self, survey_name: str, x: float, y: float, z: float):
        """
        Add subsidence monitoring point.
        
        Args:
            survey_name (str): Name of the survey
            x, y, z (float): Monitoring point coordinates
        """
    
    def eval_subsidence(self, time: datetime = None) -> numpy.ndarray:
        """
        Evaluate subsidence response.
        
        Args:
            time (datetime, optional): Evaluation time
            
        Returns:
            numpy.ndarray: Subsidence values in meters
        """
    
    def load_PRESSURE(self, pressure_kw: ResdataKW):
        """Load pressure data for geomechanical calculations."""
    
    def load_PORV(self, porv_kw: ResdataKW):
        """Load pore volume data."""
    
    def add_node(self, x: float, y: float, z: float):
        """
        Add calculation node.
        
        Args:
            x, y, z (float): Node coordinates
        """

from resdata.gravimetry import ResdataGrav
from resdata.grid import Grid
from resdata.resfile import ResdataFile
import numpy as np
import matplotlib.pyplot as plt

# Load grid and simulation data
grid = Grid("SIMULATION.EGRID")
init_file = ResdataFile("SIMULATION.INIT")
restart_file = ResdataFile("SIMULATION.UNRST")

# Get pore volume
porv = init_file.get_kw("PORV")

# Initialize gravity calculator
grav_calc = ResdataGrav(grid, porv)

# Set standard fluid densities (kg/m³)
grav_calc.add_std_density(
    oil_density=850.0,    # Light oil
    gas_density=200.0,    # Gas at reservoir conditions
    water_density=1020.0  # Formation water
)

# Define gravity survey grid
survey_points = []
x_range = np.linspace(0, 4000, 21)  # 21 points over 4 km
y_range = np.linspace(0, 3000, 16)  # 16 points over 3 km
z_survey = 0.0  # Surface elevation

for i, y in enumerate(y_range):
    for j, x in enumerate(x_range):
        survey_name = f"GRAV_{i:02d}_{j:02d}"
        grav_calc.add_survey_GRAV(survey_name, x, y, z_survey)
        survey_points.append((x, y, survey_name))

print(f"Added {len(survey_points)} gravity survey points")

# Evaluate gravity response
gravity_response = grav_calc.eval_grav()
print(f"Gravity response shape: {gravity_response.shape}")
print(f"Gravity range: {gravity_response.min():.3f} to {gravity_response.max():.3f} mGal")

# Reshape for plotting
nx, ny = len(x_range), len(y_range)
gravity_grid = gravity_response.reshape(ny, nx)

# Plot gravity map
plt.figure(figsize=(12, 9))
contour = plt.contourf(x_range, y_range, gravity_grid, levels=20, cmap='RdBu_r')
plt.colorbar(contour, label='Gravity Anomaly (mGal)')
plt.xlabel('X (m)')
plt.ylabel('Y (m)')
plt.title('Gravity Survey - Base Case')
plt.axis('equal')
plt.grid(True, alpha=0.3)
plt.show()

from resdata.gravimetry import ResdataGrav
from resdata.grid import Grid
from resdata.resfile import ResdataFile
import numpy as np
import matplotlib.pyplot as plt

# Load data
grid = Grid("SIMULATION.EGRID")
init_file = ResdataFile("SIMULATION.INIT")
restart_file = ResdataFile("SIMULATION.UNRST")

# Get base data
porv = init_file.get_kw("PORV")

# Initialize gravity calculator
grav_calc = ResdataGrav(grid, porv)
grav_calc.add_std_density(850.0, 200.0, 1020.0)

# Add survey points over reservoir
survey_points = []
for i in range(5):
    for j in range(5):
        x = 1000 + j * 500  # 500m spacing
        y = 1000 + i * 500
        z = 0.0
        survey_name = f"GRAV_{i}_{j}"
        grav_calc.add_survey_GRAV(survey_name, x, y, z)
        survey_points.append((x, y, survey_name))

# Evaluate gravity at different time steps
time_steps = [0, 365, 730, 1095]  # Days: initial, 1yr, 2yr, 3yr
gravity_maps = {}

for day in time_steps:
    # Load density data for this time step
    # (This would typically involve calculating densities from saturation data)
    
    # For demonstration, evaluate gravity response
    gravity_response = grav_calc.eval_grav()
    gravity_maps[day] = gravity_response.reshape(5, 5)
    
    print(f"Day {day}: Gravity range {gravity_response.min():.3f} to "
          f"{gravity_response.max():.3f} mGal")

# Plot time-lapse gravity changes
fig, axes = plt.subplots(2, 2, figsize=(15, 12))
axes = axes.flatten()

x_coords = [1000 + j * 500 for j in range(5)]
y_coords = [1000 + i * 500 for i in range(5)]

for idx, day in enumerate(time_steps):
    ax = axes[idx]
    im = ax.contourf(x_coords, y_coords, gravity_maps[day], 
                     levels=20, cmap='RdBu_r')
    ax.set_title(f'Day {day} ({day/365:.1f} years)')
    ax.set_xlabel('X (m)')
    ax.set_ylabel('Y (m)')
    ax.axis('equal')
    fig.colorbar(im, ax=ax, label='Gravity (mGal)')

plt.tight_layout()
plt.show()

# Calculate gravity changes
base_gravity = gravity_maps[0]
for day in time_steps[1:]:
    gravity_change = gravity_maps[day] - base_gravity
    max_change = np.max(np.abs(gravity_change))
    print(f"Maximum gravity change at day {day}: {max_change:.3f} mGal")

from resdata.gravimetry import ResdataSubsidence
from resdata.grid import Grid
from resdata.resfile import ResdataFile
import numpy as np
import matplotlib.pyplot as plt

# Load grid and data
grid = Grid("SIMULATION.EGRID")
restart_file = ResdataFile("SIMULATION.UNRST")

# Initialize subsidence calculator
subsidence_calc = ResdataSubsidence(grid)

# Add monitoring points at surface
monitoring_points = []
for i in range(11):
    for j in range(11):
        x = j * 400  # 400m spacing over 4km
        y = i * 300  # 300m spacing over 3km
        z = 0.0      # Surface level
        
        subsidence_calc.add_survey_SUBSIDENCE(f"SUB_{i}_{j}", x, y, z)
        monitoring_points.append((x, y))

print(f"Added {len(monitoring_points)} subsidence monitoring points")

# Load pressure data
pressure = restart_file.get_kw("PRESSURE", 0)  # Initial pressure
subsidence_calc.load_PRESSURE(pressure)

# Get pore volume if available
try:
    porv = restart_file.get_kw("PORV")
    subsidence_calc.load_PORV(porv)
except:
    print("PORV not found in restart file")

# Evaluate subsidence
subsidence_response = subsidence_calc.eval_subsidence()
print(f"Subsidence range: {subsidence_response.min():.4f} to "
      f"{subsidence_response.max():.4f} m")

# Reshape for plotting
subsidence_grid = subsidence_response.reshape(11, 11)

# Plot subsidence map
x_coords = [j * 400 for j in range(11)]
y_coords = [i * 300 for i in range(11)]

plt.figure(figsize=(12, 9))
contour = plt.contourf(x_coords, y_coords, subsidence_grid, 
                      levels=20, cmap='RdYlBu')
plt.colorbar(contour, label='Subsidence (m)')
plt.xlabel('X (m)')
plt.ylabel('Y (m)')
plt.title('Surface Subsidence Map')
plt.axis('equal')
plt.grid(True, alpha=0.3)

# Add contour lines
contour_lines = plt.contour(x_coords, y_coords, subsidence_grid, 
                           levels=10, colors='black', alpha=0.4, linewidths=0.5)
plt.clabel(contour_lines, inline=True, fontsize=8, fmt='%.3f')

plt.show()

# Identify maximum subsidence
max_subsidence_idx = np.unravel_index(np.argmax(subsidence_grid), subsidence_grid.shape)
max_x = x_coords[max_subsidence_idx[1]]
max_y = y_coords[max_subsidence_idx[0]]
max_value = subsidence_grid[max_subsidence_idx]

print(f"Maximum subsidence: {max_value:.4f} m at ({max_x}, {max_y})")

from resdata.gravimetry import ResdataGrav, ResdataSubsidence
from resdata.grid import Grid
from resdata.resfile import ResdataFile
import numpy as np
import matplotlib.pyplot as plt

# Load data
grid = Grid("SIMULATION.EGRID")
init_file = ResdataFile("SIMULATION.INIT")
restart_file = ResdataFile("SIMULATION.UNRST")

# Get required data
porv = init_file.get_kw("PORV")
pressure = restart_file.get_kw("PRESSURE", 0)

# Initialize both calculators
grav_calc = ResdataGrav(grid, porv)
grav_calc.add_std_density(850.0, 200.0, 1020.0)

subsidence_calc = ResdataSubsidence(grid)
subsidence_calc.load_PRESSURE(pressure)
subsidence_calc.load_PORV(porv)

# Define common survey grid
survey_grid = []
for i in range(7):
    for j in range(7):
        x = 500 + j * 500  # 500m spacing
        y = 500 + i * 500
        z = 0.0
        
        # Add to both calculators
        grav_calc.add_survey_GRAV(f"POINT_{i}_{j}", x, y, z)
        subsidence_calc.add_survey_SUBSIDENCE(f"POINT_{i}_{j}", x, y, z)
        
        survey_grid.append((x, y))

# Evaluate both responses
gravity_response = grav_calc.eval_grav()
subsidence_response = subsidence_calc.eval_subsidence()

# Reshape for plotting
gravity_grid = gravity_response.reshape(7, 7)
subsidence_grid = subsidence_response.reshape(7, 7)

# Plot combined analysis
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(18, 5))

x_coords = [500 + j * 500 for j in range(7)]
y_coords = [500 + i * 500 for i in range(7)]

# Gravity map
im1 = ax1.contourf(x_coords, y_coords, gravity_grid, levels=15, cmap='RdBu_r')
ax1.set_title('Gravity Anomaly')
ax1.set_xlabel('X (m)')
ax1.set_ylabel('Y (m)')
fig.colorbar(im1, ax=ax1, label='Gravity (mGal)')
ax1.axis('equal')

# Subsidence map
im2 = ax2.contourf(x_coords, y_coords, subsidence_grid, levels=15, cmap='RdYlBu')
ax2.set_title('Surface Subsidence')
ax2.set_xlabel('X (m)')
ax2.set_ylabel('Y (m)')
fig.colorbar(im2, ax=ax2, label='Subsidence (m)')
ax2.axis('equal')

# Correlation plot
ax3.scatter(gravity_response, subsidence_response * 1000, alpha=0.7)  # Convert to mm
ax3.set_xlabel('Gravity Anomaly (mGal)')
ax3.set_ylabel('Subsidence (mm)')
ax3.set_title('Gravity vs Subsidence Correlation')
ax3.grid(True, alpha=0.3)

# Calculate correlation coefficient
correlation = np.corrcoef(gravity_response, subsidence_response)[0, 1]
ax3.text(0.05, 0.95, f'Correlation: {correlation:.3f}', 
         transform=ax3.transAxes, verticalalignment='top',
         bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))

plt.tight_layout()
plt.show()

print(f"Gravity-Subsidence Analysis Summary:")
print(f"  Gravity range: {gravity_response.min():.3f} to {gravity_response.max():.3f} mGal")
print(f"  Subsidence range: {subsidence_response.min()*1000:.2f} to {subsidence_response.max()*1000:.2f} mm")
print(f"  Correlation coefficient: {correlation:.3f}")

from resdata.gravimetry import ResdataGrav
from resdata.grid import Grid
from resdata.resfile import ResdataFile
import numpy as np
import matplotlib.pyplot as plt

# Load data
grid = Grid("SIMULATION.EGRID")
init_file = ResdataFile("SIMULATION.INIT")
porv = init_file.get_kw("PORV")

# Initialize gravity calculator
grav_calc = ResdataGrav(grid, porv)
grav_calc.add_std_density(850.0, 200.0, 1020.0)

# Test different survey designs
survey_designs = {
    'coarse': {'spacing': 1000, 'description': '1km spacing'},
    'medium': {'spacing': 500, 'description': '500m spacing'},
    'fine': {'spacing': 250, 'description': '250m spacing'}
}

results = {}

for design_name, design in survey_designs.items():
    spacing = design['spacing']
    
    # Create new gravity calculator for this design
    grav_test = ResdataGrav(grid, porv)
    grav_test.add_std_density(850.0, 200.0, 1020.0)
    
    # Add survey points
    survey_points = []
    x_range = np.arange(0, 4001, spacing)
    y_range = np.arange(0, 3001, spacing)
    
    for i, y in enumerate(y_range):
        for j, x in enumerate(x_range):
            grav_test.add_survey_GRAV(f"P_{i}_{j}", x, y, 0.0)
            survey_points.append((x, y))
    
    # Evaluate gravity
    gravity_response = grav_test.eval_grav()
    
    results[design_name] = {
        'points': len(survey_points),
        'spacing': spacing,
        'gravity_range': gravity_response.max() - gravity_response.min(),
        'gravity_std': np.std(gravity_response),
        'response': gravity_response
    }

# Compare survey designs
print("Survey Design Comparison:")
print(f"{'Design':<10} {'Points':<8} {'Spacing':<10} {'Range':<12} {'Std Dev':<10}")
print("-" * 60)

for name, result in results.items():
    print(f"{name:<10} {result['points']:<8} {result['spacing']:<10} "
          f"{result['gravity_range']:<12.3f} {result['gravity_std']:<10.3f}")

# Plot survey design comparison
fig, axes = plt.subplots(1, 3, figsize=(18, 5))

for idx, (name, result) in enumerate(results.items()):
    ax = axes[idx]
    
    spacing = result['spacing']
    x_range = np.arange(0, 4001, spacing)
    y_range = np.arange(0, 3001, spacing)
    
    gravity_grid = result['response'].reshape(len(y_range), len(x_range))
    
    im = ax.contourf(x_range, y_range, gravity_grid, levels=15, cmap='RdBu_r')
    ax.set_title(f"{name.title()} Grid ({result['points']} points)")
    ax.set_xlabel('X (m)')
    ax.set_ylabel('Y (m)')
    fig.colorbar(im, ax=ax, label='Gravity (mGal)')
    ax.axis('equal')

plt.tight_layout()
plt.show()

# Cost-benefit analysis
print("\nCost-Benefit Analysis:")
base_cost_per_point = 1000  # Cost per gravity measurement
for name, result in results.items():
    total_cost = result['points'] * base_cost_per_point
    resolution = result['gravity_std']
    efficiency = resolution / (total_cost / 1000)  # Resolution per k$
    
    print(f"{name:<10}: Cost ${total_cost:,}, "
          f"Resolution {resolution:.3f} mGal, "
          f"Efficiency {efficiency:.6f} mGal/k$")

# Gravity and subsidence response arrays
GravityResponse = numpy.ndarray      # Gravity values in milligals
SubsidenceResponse = numpy.ndarray   # Subsidence values in meters

# Survey point coordinates
SurveyPoint = tuple[float, float, float]  # (x, y, z)
SurveyGrid = list[SurveyPoint]

# Fluid density parameters
FluidDensities = dict[str, float]    # Keys: 'oil', 'gas', 'water' (kg/m³)

# Geophysical survey metadata
SurveyMetadata = dict[str, any]      # Survey configuration and parameters