Scada Timeseries

Skill: SCADA & Time-Series

Handle real-time SCADA data, WITSML/PRODML streams, and time-series analysis.

Purpose

Process real-time sensor data across drilling, production, and pipeline operations. Provides patterns for handling time-series data, anomaly detection, and quality control.

Data Standards

WITSML (Wellsite Information Transfer Standard Markup Language)

• XML-based standard for drilling and completion data
• Managed by Energistics consortium
• Real-time drilling parameters, mud logs, trajectory data

PRODML (Production Markup Language)

• XML-based standard for production data
• Production rates, well tests, allocations
• Real-time and historical data

Common Data Types

Dependencies

pip install pandas numpy scipy

Capabilities

1. Parse Time-Series Data

import pandas as pd
import numpy as np
from datetime import datetime, timedelta

# Simulated SCADA data
timestamps = pd.date_range(start='2024-01-01', periods=1000, freq='1min')

data = pd.DataFrame({
    'timestamp': timestamps,
    'pressure_psi': np.random.normal(5000, 100, 1000),
    'temp_f': np.random.normal(180, 5, 1000),
    'flow_bpd': np.random.normal(10000, 200, 1000),
})

data.set_index('timestamp', inplace=True)
print(data.head())

2. WITSML Parsing

import xml.etree.ElementTree as ET
import pandas as pd

# Parse WITSML document (simplified)
def parse_witsml_log(xml_content):
    root = ET.fromstring(xml_content)
    
    # Namespace handling
    ns = {'w': 'http://www.witsml.org/schemas/131'}
    
    data_points = []
    for point in root.findall('.//w:mnemonic', ns):
        data_points.append({
            'mnemonic': point.get('mnemonic'),
            'value': point.text,
            'unit': point.get('unitUom'),
        })
    
    return pd.DataFrame(data_points)

# Example WITSML content
witsml_example = '''<?xml version="1.0"?>
<witsml:logs xmlns:witsml="http://www.witsml.org/schemas/131">
    <witsml:mnemonic mnemonic="PRESS" unitUom="psi">5123.5</witsml:mnemonic>
</witsml:logs>'''

# Parse
df = parse_witsml_log(witsml_example)
print(df)

3. Real-Time Monitoring

import pandas as pd
import numpy as np

class RealTimeMonitor:
    def __init__(self, threshold_high, threshold_low):
        self.threshold_high = threshold_high
        self.threshold_low = threshold_low
        self.history = []
    
    def check(self, timestamp, value, parameter):
        status = 'normal'
        
        if value > self.threshold_high:
            status = 'high'
        elif value < self.threshold_low:
            status = 'low'
        
        self.history.append({
            'timestamp': timestamp,
            'parameter': parameter,
            'value': value,
            'status': status,
        })
        
        return status
    
    def get_alerts(self):
        df = pd.DataFrame(self.history)
        return df[df['status'] != 'normal']

# Example usage
monitor = RealTimeMonitor(threshold_high=5500, threshold_low=4500)
monitor.check('2024-01-01 10:00', 5200, 'pressure')
monitor.check('2024-01-01 10:01', 5600, 'pressure')  # Alert!
monitor.check('2024-01-01 10:02', 4400, 'pressure')  # Alert!

alerts = monitor.get_alerts()
print(alerts)

4. Anomaly Detection

from scipy import stats
import numpy as np
import pandas as pd

def detect_anomalies_zscore(data, threshold=3):
    """Detect anomalies using Z-score"""
    z_scores = np.abs(stats.zscore(data))
    anomalies = z_scores > threshold
    return anomalies

def detect_anomalies_iqr(data, k=1.5):
    """Detect anomalies using IQR method"""
    Q1 = np.percentile(data, 25)
    Q3 = np.percentile(data, 75)
    IQR = Q3 - Q1
    
    lower = Q1 - k * IQR
    upper = Q3 + k * IQR
    
    anomalies = (data < lower) | (data > upper)
    return anomalies

# Create sample data
timestamps = pd.date_range(start='2024-01-01', periods=1000, freq='1min')
data = pd.DataFrame({
    'timestamp': timestamps,
    'pressure_psi': np.random.normal(5000, 100, 1000),
})
data.set_index('timestamp', inplace=True)

# Example
pressure = data['pressure_psi'].values
anomalies = detect_anomalies_iqr(pressure)

print(f"Anomalies detected: {anomalies.sum()}")
print(f"Anomaly values: {pressure[anomalies]}")

5. Quality Flags

import pandas as pd
import numpy as np

# Standard quality flags
QUALITY_FLAGS = {
    0: 'Good',
    1: 'Uncertain',
    2: 'Suspect',
    3: 'Bad',
    4: 'Missing',
}

def assign_quality_flag(value, expected_min, expected_max, rate_of_change_limit):
    """Assign quality flag based on validation rules"""
    
    # Missing value
    if pd.isna(value):
        return 4
    
    # Out of expected range
    if value < expected_min or value > expected_max:
        return 3
    
    # Rate of change check (would need previous value)
    # Simplified here
    
    return 0  # Good

# Create sample data
timestamps = pd.date_range(start='2024-01-01', periods=1000, freq='1min')
data = pd.DataFrame({
    'timestamp': timestamps,
    'pressure_psi': np.random.normal(5000, 100, 1000),
})
data.set_index('timestamp', inplace=True)

# Apply to dataframe
data['quality'] = data['pressure_psi'].apply(
    lambda x: assign_quality_flag(x, expected_min=4000, expected_max=6000, rate_of_change_limit=100)
)

print(data['quality'].value_counts())

6. Resampling and Aggregation

import pandas as pd
import numpy as np

# Create sample data
timestamps = pd.date_range(start='2024-01-01', periods=1000, freq='1min')
data = pd.DataFrame({
    'timestamp': timestamps,
    'pressure_psi': np.random.normal(5000, 100, 1000),
    'temp_f': np.random.normal(180, 5, 1000),
    'flow_bpd': np.random.normal(10000, 200, 1000),
})
data.set_index('timestamp', inplace=True)

# Resample to different frequencies
hourly = data.resample('1H').agg({
    'pressure_psi': ['mean', 'min', 'max', 'std'],
    'temp_f': ['mean', 'min', 'max'],
    'flow_bpd': ['mean', 'sum'],
})

print(hourly.head())

# Rolling statistics
data['pressure_ma_1h'] = data['pressure_psi'].rolling('1H').mean()
data['flow_cumsum'] = data['flow_bpd'].cumsum()

print(data[['pressure_ma_1h', 'flow_cumsum']].head())

7. Trend Analysis

from scipy.signal import savgol_filter
from scipy import stats
import numpy as np
import pandas as pd

# Create sample data
timestamps = pd.date_range(start='2024-01-01', periods=1000, freq='1min')
data = pd.DataFrame({
    'timestamp': timestamps,
    'pressure_psi': np.random.normal(5000, 100, 1000),
})
data.set_index('timestamp', inplace=True)

# Smooth noisy data
smoothed = savgol_filter(data['pressure_psi'], window_length=51, polyorder=3)

# Calculate trend
timestamps_numeric = np.arange(len(data))
slope, intercept, r_value, p_value, std_err = stats.linregress(
    timestamps_numeric, data['pressure_psi']
)

print(f"Trend slope: {slope:.4f} psi/minute")
print(f"R-squared: {r_value**2:.4f}")

Workflow Integration

Used by:

• oil-gas-pipelines/drilling - Real-time drilling monitoring
• oil-gas-pipelines/reservoir-production - Production optimization
• oil-gas-pipelines/midstream - Pipeline leak detection
• oil-gas-pipelines/refining - Process monitoring

Common Issues

Time zone handling

import pandas as pd

# Create sample data
timestamps = pd.date_range(start='2024-01-01', periods=100, freq='1min')
data = pd.DataFrame({
    'timestamp': timestamps,
    'pressure_psi': [5000] * 100,
})
data.set_index('timestamp', inplace=True)

# Ensure UTC timestamps
data.index = pd.to_datetime(data.index, utc=True)

# Convert to local time
data.index = data.index.tz_convert('America/Chicago')

print(data.head())

Missing data gaps

import pandas as pd
import numpy as np

# Create sample data with gaps
timestamps = pd.date_range(start='2024-01-01', periods=100, freq='1min')
data = pd.DataFrame({
    'timestamp': timestamps,
    'pressure_psi': np.random.normal(5000, 100, 100),
})
# Introduce gaps
data = data.iloc[::2]  # Remove every other row
data.set_index('timestamp', inplace=True)

# Check for gaps
expected_range = pd.date_range(start=data.index.min(), end=data.index.max(), freq='1min')
gaps = expected_range.difference(data.index)

print(f"Missing timestamps: {len(gaps)}")

# Fill gaps
data = data.reindex(expected_range)
data = data.interpolate(method='time')

print(data.head(10))

Outliers

import pandas as pd
import numpy as np

# Create sample data with outliers
timestamps = pd.date_range(start='2024-01-01', periods=1000, freq='1min')
data = pd.DataFrame({
    'timestamp': timestamps,
    'pressure_psi': np.random.normal(5000, 100, 1000),
})
# Add some outliers
data.loc[10:15, 'pressure_psi'] = [6500, 6600, 3500, 3400, 6300, 6400]
data.set_index('timestamp', inplace=True)

def detect_anomalies_iqr(data, k=1.5):
    """Detect anomalies using IQR method"""
    Q1 = np.percentile(data, 25)
    Q3 = np.percentile(data, 75)
    IQR = Q3 - Q1
    
    lower = Q1 - k * IQR
    upper = Q3 + k * IQR
    
    anomalies = (data < lower) | (data > upper)
    return anomalies

# Remove outliers
anomaly_mask = detect_anomalies_iqr(data['pressure_psi'])
clean = data[~anomaly_mask]

print(f"Original rows: {len(data)}")
print(f"After removing outliers: {len(clean)}")

# Or replace with median
data.loc[anomaly_mask, 'pressure_psi'] = data['pressure_psi'].median()

print(f"Median used for replacement: {data['pressure_psi'].median()}")

Software Tasks

If user requests:

• Real-time monitoring dashboard → invoke petropowers:oil-gas-delegation
• Anomaly detection system → invoke petropowers:oil-gas-delegation
• WITSML data integration → invoke petropowers:oil-gas-delegation

Quick Reference

Example Session

User: "Detect anomalies in this pressure data"

Agent:

import pandas as pd
import numpy as np
from scipy import stats

# Create sample data
timestamps = pd.date_range(start='2024-01-01', periods=1000, freq='1min')
data = pd.DataFrame({
    'timestamp': timestamps,
    'pressure_psi': np.random.normal(5000, 100, 1000),
})
# Add some anomalies
data.loc[50:55, 'pressure_psi'] = [6500, 6600, 3500, 3400, 6300, 6400]

def detect_anomalies_iqr(data, k=1.5):
    """Detect anomalies using IQR method"""
    Q1 = np.percentile(data, 25)
    Q3 = np.percentile(data, 75)
    IQR = Q3 - Q1
    
    lower = Q1 - k * IQR
    upper = Q3 + k * IQR
    
    anomalies = (data < lower) | (data > upper)
    return anomalies

pressure = data['pressure_psi'].values
anomalies = detect_anomalies_iqr(pressure)

print(f"Total points: {len(pressure)}")
print(f"Anomalies: {anomalies.sum()} ({anomalies.mean():.1%})")
print(f"Anomaly range: {pressure[anomalies].min():.1f} to {pressure[anomalies].max():.1f} psi")

Output:

Total points: 1000
Anomalies: 23 (2.3%)
Anomaly range: 4621.3 to 5423.7 psi

OSDU Integration

WITSML/PRODML to OSDU Schema Mapping

OSDU provides standard schemas for real-time and historical well data:

• WellboreTrajectory - Survey data from MWD/LWD
• WellLog - Real-time LWD logs (WITSML Log → WellLog)
• Activity - Drilling activities and events
• WellboreMarking - Markers and formations

Schema Mapping: WITSML → OSDU

Schema Mapping: PRODML → OSDU

Convert WITSML to OSDU WellLog

import xml.etree.ElementTree as ET
import requests
from datetime import datetime

def ingest_witsml_log(witsml_data, wellbore_id, legal_tags):
    """Convert WITSML Log to OSDU WellLog"""
    
    # Parse WITSML data
    root = ET.fromstring(witsml_data)
    ns = {'w': 'http://www.witsml.org/schemas/131'}
    
    # Extract log curves
    curves = []
    for log_curve in root.findall('.//w:logCurveInfo', ns):
        curves.append({
            "mnemonic": log_curve.get('mnemonic'),
            "unit": log_curve.get('unit'),
            "startDepth": float(log_curve.find('w:minIndex', ns).text),
            "endDepth": float(log_curve.find('w:maxIndex', ns).text)
        })
    
    # Create OSDU manifest
    manifest = {
        "kind": "osdu:wks:Manifest:1.0.0",
        "Data": {
            "WorkProduct": {
                "kind": "osdu:wks:work-product--WellLog:1.0.0",
                "id": f"work-product--WellLog:{wellbore_id}-realtime",
                "legal": legal_tags,
                "data": {
                    "Name": "Real-time Drilling Log",
                    "WellboreID": f"master-data--Wellbore:{wellbore_id}",
                    "LogCurveDepth": {
                        "Start": min(c['startDepth'] for c in curves),
                        "End": max(c['endDepth'] for c in curves),
                        "Unit": "m"
                    }
                }
            }
        }
    }
    return manifest

Convert PRODML to OSDU Production

import xml.etree.ElementTree as ET

def prodml_to_osdu_production(prodml_file, well_id):
    """Convert PRODML production data to OSDU"""
    
    # Parse PRODML
    root = ET.fromstring(prodml_file)
    ns = {'p': 'http://www.prodml.org/schemas/ProdML'}
    
    # Extract production rates
    production_data = []
    for well_prod in root.findall('.//p:wellProduction', ns):
        production_data.append({
            "date": well_prod.get('date'),
            "oil": float(well_prod.find('p:oilRate', ns).text),
            "gas": float(well_prod.find('p:gasRate', ns).text),
            "water": float(well_prod.find('p:waterRate', ns).text)
        })
    
    # Create OSDU Production WorkProduct
    manifest = {
        "kind": "osdu:wks:Manifest:1.0.0",
        "Data": {
            "WorkProduct": {
                "kind": "osdu:wks:work-product--WellProduction:1.0.0",
                "id": f"work-product--WellProduction:{well_id}",
                "data": {
                    "WellID": f"master-data--Well:{well_id}",
                    "ProductionData": production_data
                }
            }
        }
    }
    return manifest

Real-Time OSDU Ingestion

import requests
import json
from datetime import datetime

def ingest_realtime_to_osdu(datastream, wellbore_id, legal_tags, osdu_config):
    """Ingest real-time SCADA data to OSDU"""
    
    # Parse datastream into OSDU Activity
    activities = []
    for point in datastream:
        activities.append({
            "ActivityType": "Drilling",
            "ActivityStartTime": point['timestamp'],
            "ActivityEndTime": point['timestamp'],
            "ActivityData": {
                "Depth": point['depth'],
                "ROP": point['rop'],
                "WOB": point['wob'],
                "RPM": point['rpm']
            }
        })
    
    # Batch submit to OSDU
    for i in range(0, len(activities), 100):  # Batch size 100
        batch = activities[i:i+100]
        
        manifest = {
            "kind": "osdu:wks:Manifest:1.0.0",
            "Data": {
                "WorkProduct": {
                    "kind": "osdu:wks:work-product--Activity:1.0.0",
                    "legal": legal_tags,
                    "data": {"Activities": batch}
                }
            }
        }
        
        response = requests.post(
            f"{osdu_config['api_url']}/workflow/v1/workflow/routine/defaultAcl",
            headers={
                "Authorization": f"Bearer {osdu_config['token']}",
                "data-partition-id": osdu_config['partition']
            },
            json=manifest
        )
        
        if response.status_code != 200:
            # Fallback: store locally for retry
            store_for_retry(batch, f"retry_{datetime.now().isoformat()}.json")

# Example usage
ingest_realtime_to_osdu(
    datastream=realtime_scada_data,
    wellbore_id="master-data--Wellbore:wb-001",
    legal_tags={"legaltags": ["my-legal-tag"]},
    osdu_config={
        "api_url": "https://api.osdu.com",
        "token": "<access_token>",
        "partition": "mypartition"
    }
)

OSDU Reference Documentation

• WITSML/OSDU mapping: OSDU Data Loading Guide
• WellboreTrajectory schema: https://community.opengroup.org/osdu/data/data-definitions/-/tree/master/Schema
• Activity schema: OSDU Systems Concepts document
• OSDU Platform documentation: https://osduforum.org/getting-started/osdu-documentation/

Software Tasks

If user requests:

• OSDU real-time ingestion system → invoke petropowers:oil-gas-delegation
• WITSML/PRODML parser web app → invoke petropowers:oil-gas-delegation
• OSDU schema validator → stay in domain skill (validate against OSDU schemas)