# IMPORT LIBRARIES
import pygmt
import pandas as pd
import numpy as np
import sys

# Show versions
print("Python version:", sys.version)
print("PyGMT version:", pygmt.__version__)
print("Pandas version:", pd.__version__)
print("NumPy version:", np.__version__)

Python version: 3.11.13 | packaged by conda-forge | (main, Jun  4 2025, 14:39:58) [MSC v.1943 64 bit (AMD64)]
PyGMT version: v0.17.0
Pandas version: 2.3.3
NumPy version: 2.3.3

# Read data (Pennsylvania). Assumption is that all unconventional wells are fracked.
pa_data = pd.read_csv("PA_Unc_WellsSummary.csv")
pa_lats = pa_data["LATITUDE"]
pa_lons = pa_data["LONGITUDE"]

# Read data (New York). Assumption is that all horizontal wells are fracked.
ny_data = pd.read_csv("NY_wellspublic.csv", low_memory=False)
horizontal_ny = ny_data[ny_data["Slant"].str.contains("Horizontal", case=False, na=False)]
ny_lats_horizontal = horizontal_ny["Surface_latitude"]
ny_lons_horizontal = horizontal_ny["Surface_Longitude"]


# FIRST FIGURE
fig = pygmt.Figure()

# Define region (centered on New York and Pennsylvania)
region = [-82, -72, 39, 45]

# Base map with local projection (Lambert Conformal)
fig.basemap(region=region, projection="L-77/42/33/45/12c", frame=["a", "+tPA & NY Fracked Wells (Points)"])
fig.coast(land="green", water="lightblue", shorelines="1/0.5p", borders=["1/0.5p,black", "2/0.5p,black" ])

# Plot all filtered wells 
fig.plot(x=pa_lons, y=pa_lats, style="c0.1c", fill="red", pen="black", label="PA")
fig.plot(x=ny_lons_horizontal, y=ny_lats_horizontal, style="c0.1c", fill="blue", pen="black", label="NY")

# Add legend to distibguish wells in each state
fig.legend(position="JBR+jBR+o0.2c")


# SECOND FIGURE
fig.shift_origin(xshift="w+3c")

# Base map with local projection (Lambert Conformal)
fig.basemap(region=region, projection="L-77/42/33/45/12c", frame=["ag", "+tPA & NY Fracked Wells (Density)"])

# Combine PA and NY wells
all_lons = np.concatenate([pa_lons.values, ny_lons_horizontal.values])
all_lats = np.concatenate([pa_lats.values, ny_lats_horizontal.values])

# Remove any NaNs
mask = ~np.isnan(all_lons) & ~np.isnan(all_lats)
all_lons = all_lons[mask]
all_lats = all_lats[mask]

# The xyz2grd function expects z values, so we'll create a dummy array of ones with ones_like.
z = np.ones_like(all_lons)
z = np.ones_like(all_lats)

# Create grid
grid = pygmt.xyz2grd(
    x=all_lons,
    y=all_lats,
    z=z,
    region=region,
    spacing=0.1,
    duplicate="n",      # 'n' tells it to count points per cell, cellsize dictated by spacing param
    registration="g"    # gridline registration
)

# Add grid to map
fig.grdimage(grid)

# Add coastlines/borders
fig.coast(shorelines="1/0.5p", borders=["1/0.5p,black", "2/0.5p,black" ])

# Add colorbar
fig.colorbar(frame=["x+l# of wells in bin"], position="JMR+o0.5c/0c+w8c")

# Show figure
#fig.show()
#fig.savefig("frac_map.jpg")

#PA Unconventional wells data source: https://maps.carnegiemnh.org/index.php/projects/unconventional-wells/download/
#NY Horizontal wells data source: https://app.box.com/s/i7w2tm3tlp4fqmoe3pzelyjx5gbjj2lk/folder/138487259453