library(sf)## Warning: package 'sf' was built under R version 4.0.3## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1library(dplyr)## Warning: package 'dplyr' was built under R version 4.0.3##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(downloader)## Warning: package 'downloader' was built under R version 4.0.3library(tigris)## Warning: package 'tigris' was built under R version 4.0.3## To enable
## caching of data, set `options(tigris_use_cache = TRUE)` in your R script or .Rprofile.#Download zip file of wells
#downloader::download('ftp://164.64.106.6/Public/OCD/OCD%20GIS%20Data/Shapefiles/WellGIS.zip', destfile = 'WellGIS.zip', mode = 'wb')
#Unzip it to your directory
#unzip('WellGIS.zip')
#Delete zip file
#unlink('WellGIS.zip')
#Read into a sf dataframe
sf2 <- sf::read_sf('/Users/john/Documents/WellsGIS.shp')
sf2## Simple feature collection with 125475 features and 29 fields (with 7 geometries empty)
## geometry type: POINT
## dimension: XY
## bbox: xmin: -108.9999 ymin: 32.00003 xmax: -103.0238 ymax: 36.99971
## geographic CRS: NAD83
## # A tibble: 125,475 x 30
## id name type_code type status_cod status ogrid ogrid_name district_c
## <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <chr> <dbl>
## 1 30-0~ EUNI~ O Oil A Active 5380 XTO ENERG~ 1
## 2 30-0~ UTE ~ O Oil A Active 148219 KEYSTONE ~ 3
## 3 30-0~ EUNI~ I Inje~ P Plugg~ 5380 XTO ENERG~ 1
## 4 30-0~ HATC~ O Oil A Active 371618 RELIABLE ~ 3
## 5 30-0~ RINC~ G Gas A Active 372286 ENDURING ~ 3
## 6 30-0~ YAGE~ G Gas A Active 372286 ENDURING ~ 3
## 7 30-0~ NEW ~ G Gas A Active 5380 XTO ENERG~ 1
## 8 30-0~ ARNO~ O Oil P Plugg~ 5380 XTO ENERG~ 1
## 9 30-0~ NORT~ I Inje~ A Active 21355 SOUTHWEST~ 2
## 10 30-0~ ARNO~ O Oil A Active 21355 SOUTHWEST~ 2
## # ... with 125,465 more rows, and 21 more variables: district <chr>,
## # county_cod <dbl>, county <chr>, ulstr <chr>, latitude <dbl>,
## # longitude <dbl>, projection <chr>, directiona <chr>, details <chr>,
## # files <chr>, status2 <chr>, symbology <chr>, year_spudd <chr>,
## # spud_date <date>, measured_v <dbl>, true_verti <dbl>, pool_id_li <chr>,
## # effective_ <date>, plug_date <date>, GlobalID <chr>, geometry <POINT [°]>#Filter counties
sf2 <- sf2 %>% filter(county == 'Eddy'|county == 'Lea')
#Create new reservoir name column; I only really care about Wolfcamp/Bone Spring
sf2$reservoir[grepl('WOLFC', sf2$pool_id_li)] <- 'Wolfcamp'## Warning: Unknown or uninitialised column: `reservoir`.sf2$reservoir[grepl('BONE S', sf2$pool_id_li)] <- 'Bone Spring'
#Clean up Section/Twn/Rng column
sf2$ulstr[nchar(sf2$ulstr) == 12] <- substr(sf2$ulstr[nchar(sf2$ulstr)==12], 3, 12)
#Filter to non P&A'd wells
sf2 <- sf2 %>% filter(status == 'Active'|status == 'New'|grepl('emporary', status))
#Create new operator column
sf2$operator <- sf2$ogrid_name
sf2$operator[grepl('CHEVRON', sf2$operator)] <- 'CHEVRON'
sf2$operator[grepl('Chevron', sf2$operator)] <- 'CHEVRON'
sf2$operator[grepl('DEVON', sf2$operator)] <- 'DEVON'
sf2$operator[grepl('OXY', sf2$operator)] <- 'OXY'
sf2$operator[grepl('OCCIDENTAL', sf2$operator)] <- 'OXY'
sf2$operator[grepl('XTO', sf2$operator)] <- 'XTO'
sf2$operator[grepl('CIMAREX', sf2$operator)] <- 'CIMAREX'
sf2$operator[grepl('COG O', sf2$operator)] <- 'CONCHO'
sf2$operator[grepl('COG P', sf2$operator)] <- 'CONCHO'
sf2$operator[grepl('EOG', sf2$operator)] <- 'EOG'
sf2$operator[grepl('MATADOR', sf2$operator)] <- 'MATADOR'
sf2$operator[grepl('Spur', sf2$operator)] <- 'SPUR'
sf2$operator[grepl('TAP ROCK', sf2$operator)] <- 'TAP ROCK'
sf2$operator[grepl('CONOCO', sf2$operator)] <- 'CONOCO'
sf2$operator[grepl('MARATHON', sf2$operator)] <- 'MARATHON'library(sf)
library(dplyr)
library(downloader)
library(tigris)
library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.0.3## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --## v ggplot2 3.3.2 v purrr 0.3.4
## v tibble 3.0.4 v stringr 1.4.0
## v tidyr 1.1.2 v forcats 0.5.0
## v readr 1.4.0## Warning: package 'ggplot2' was built under R version 4.0.3## Warning: package 'tibble' was built under R version 4.0.3## Warning: package 'tidyr' was built under R version 4.0.3## Warning: package 'readr' was built under R version 4.0.3## Warning: package 'purrr' was built under R version 4.0.3## Warning: package 'stringr' was built under R version 4.0.3## Warning: package 'forcats' was built under R version 4.0.3## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()downloader::download('https://www.nm.blm.gov/shapeFiles/state_wide/SURFACE_OWN.zip', destfile = 'SURFACE_OWN.zip', mode = 'wb')
unzip('SURFACE_OWN.zip')
unlink('SURFACE_OWN.zip')
#Read in sf dataframe
surface <- read_sf('Surface_Ownership.shp')
#Grab New Mexico county data
counties <- tigris::counties(state = 'NM')##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========= | 14%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|============== | 21%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================ | 24%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================= | 34%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|============================== | 44%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 46%
|
|================================= | 47%
|
|================================= | 48%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 49%
|
|=================================== | 50%
|
|=================================== | 51%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|===================================== | 54%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|========================================== | 61%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================ | 64%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|=================================================== | 74%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 79%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|========================================================== | 84%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================= | 94%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%counties <- st_as_sf(counties)
#Transform both to WGS84
surface <- surface %>% st_transform(crs = 4326)
counties <- counties %>% st_transform(crs = 4326)
#Filter to Eddy/Lea
counties <- counties %>% filter(NAME == 'Eddy'|NAME == 'Lea') %>% dplyr::select(county = NAME)
#Join county file and remove any na's
surface <- surface %>% st_join(counties)## although coordinates are longitude/latitude, st_intersects assumes that they are planarsurface <- surface %>% filter(!is.na(county))
surface <- surface %>% filter(!duplicated(globalid))library(leaflet)## Warning: package 'leaflet' was built under R version 4.0.5leaflet(options = leafletOptions(zoomControl = FALSE)) %>%
addProviderTiles(providers$OpenStreetMap) %>%
addPolygons(data = counties, fillColor = 'transparent', color = 'black', weight = 2) %>%
addPolygons(data = surface %>% filter(own == 'BLM'), fillColor = 'red', fillOpacity = 0.25, color = 'black', weight = 1)%>%
addMiniMap(position = 'topleft')wcList <- sf2 %>% filter(reservoir == 'Wolfcamp') %>%
st_transform(4326) %>% st_join(surface) %>%
filter(!duplicated(ulstr)) %>%
mutate(own = replace(own, own != 'BLM', 'Other')) %>%
data.frame() %>% group_by(own) %>% summarise(count = n()*640) %>%
filter(!is.na(own))## although coordinates are longitude/latitude, st_intersects assumes that they are planar## `summarise()` ungrouping output (override with `.groups` argument)boneSp <- sf2 %>% filter(reservoir == 'Bone Spring') %>%
st_transform(4326) %>% st_join(surface) %>%
filter(!duplicated(ulstr)) %>%
mutate(own = replace(own, own != 'BLM', 'Other')) %>%
data.frame() %>% group_by(own) %>% summarise(count = n()*640) %>%
filter(!is.na(own))## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## `summarise()` ungrouping output (override with `.groups` argument)library("data.table")##
## Attaching package: 'data.table'## The following object is masked from 'package:purrr':
##
## transpose## The following objects are masked from 'package:dplyr':
##
## between, first, lastlibrary("sp")## Warning: package 'sp' was built under R version 4.0.3library("rgdal")## Warning: package 'rgdal' was built under R version 4.0.3## rgdal: version: 1.5-18, (SVN revision 1082)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: C:/Users/john/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: C:/Users/john/Documents/R/win-library/4.0/rgdal/proj
## Linking to sp version:1.4-4
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.library("KernSmooth")## KernSmooth 2.23 loaded
## Copyright M. P. Wand 1997-2009library("raster")## Warning: package 'raster' was built under R version 4.0.5##
## Attaching package: 'raster'## The following object is masked from 'package:data.table':
##
## shift## The following object is masked from 'package:tidyr':
##
## extract## The following object is masked from 'package:dplyr':
##
## selectlibrary("ggplot2")
library("RColorBrewer")## Warning: package 'RColorBrewer' was built under R version 4.0.3#Convert sf dataframe to data table after filtering to just Wolfcamp wells
dat<- data.table(sf2 %>% filter(reservoir == 'Wolfcamp'))
#Create Kernel Density Output
kde <- bkde2D(dat[ , list(longitude, latitude)],
bandwidth=c(.0045, .0068), gridsize = c(500,500))
#Convert kde to a Raster for plot
KernelDensityRaster <- raster(list(x=kde$x1 ,y=kde$x2 ,z = kde$fhat))
KernelDensityRaster@data@values[which(KernelDensityRaster@data@values < 1)] <- NA
#Color scale for the raster object (11 bins)
palRaster <- colorBin(rev(brewer.pal(n=11, name='Spectral')),bins = 10,
domain = KernelDensityRaster@data@values, na.color = "transparent")
#Create the leaflet Plot
leaflet(options = leafletOptions(zoomControl = FALSE)) %>% addTiles() %>%
addRasterImage(KernelDensityRaster,
colors = palRaster,
opacity = .8) %>%
addLegend(pal = palRaster,
values = KernelDensityRaster@data@values,
title = "Well Density<br>Wolfcamp")%>%
addPolygons(data = counties, fillColor = 'transparent', color = 'black', weight = 2) %>% addPolygons(data = surface %>% filter(own == 'BLM'), fillColor = 'red', fillOpacity = 0.15, color = 'black', weight = 1)## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in CRS definition## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in CRS definitionlibrary("data.table")
library("sp")
library("rgdal")
library("KernSmooth")
library("raster")
library("ggplot2")
library("RColorBrewer")
#Convert sf dataframe to data table after filtering to just Bone Spring wells
dat<- data.table(sf2 %>% filter(reservoir == 'Bone Spring'))
#Create Kernel Density Output
kde <- bkde2D(dat[ , list(longitude, latitude)],
bandwidth=c(.0045, .0068), gridsize = c(500,500))
#Convert kde to a Raster for plot
KernelDensityRaster <- raster(list(x=kde$x1 ,y=kde$x2 ,z = kde$fhat))
KernelDensityRaster@data@values[which(KernelDensityRaster@data@values < 1)] <- NA
#Color scale for the raster object (11 bins)
palRaster <- colorBin(rev(brewer.pal(n=11, name='Spectral')),bins = 10,
domain = KernelDensityRaster@data@values, na.color = "transparent")
#Create the leaflet Plot
leaflet(options = leafletOptions(zoomControl = FALSE)) %>% addTiles() %>%
addRasterImage(KernelDensityRaster,
colors = palRaster,
opacity = .8) %>%
addLegend(pal = palRaster,
values = KernelDensityRaster@data@values,
title = "Well Density<br>BoneSpring")%>%
addPolygons(data = counties, fillColor = 'transparent', color = 'black', weight = 2) %>% addPolygons(data = surface %>% filter(own == 'BLM'), fillColor = 'red', fillOpacity = 0.15, color = 'black', weight = 1)## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in CRS definition## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in CRS definitionlibrary(tidyr)
#Filter to south of Carlsbad
wc1 <- sf2 %>% filter(latitude < 32.42)
#Find sections and operators in the data
operatorList <- wc1 %>% data.frame() %>%
group_by(ulstr, operator) %>% summarise(wells = n()) %>%
ungroup() %>% arrange(operator)## `summarise()` regrouping output by 'ulstr' (override with `.groups` argument)#Join the geometry object to the operatorList dataframe
locs <- sf2 %>% filter(!duplicated(ulstr))
locs <- locs[,c('ulstr', 'geometry')]
operatorList <- operatorList %>% left_join(locs)## Joining, by = "ulstr"#Transform to WGS84
operatorList <- st_as_sf(operatorList) %>% st_transform(crs = 4326)
#Join BLM grid to dataframe
operatorList <- st_as_sf(operatorList) %>% st_join(surface)## although coordinates are longitude/latitude, st_intersects assumes that they are planar## although coordinates are longitude/latitude, st_intersects assumes that they are planar#Remove duplicates
operatorList <- operatorList %>% filter(!duplicated(paste0(operator, ulstr)))
#Find operators with at least 40000 total acres and get split of BLM/Other
operatorList <- operatorList %>%
mutate(own = replace(own, own != 'BLM', 'Other')) %>%
data.frame() %>% group_by(operator, own) %>% summarise(count = n()*640) %>%
filter(!is.na(own)) %>% group_by(operator) %>% filter(sum(count) >= 40000) %>% ungroup() %>%
spread(own, count) %>%
mutate(blmFrac = BLM/(BLM+Other), total = BLM+Other)## `summarise()` regrouping output by 'operator' (override with `.groups` argument)Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.
##Midland Side Research
library(tidyverse)
library(rvest)## Warning: package 'rvest' was built under R version 4.0.3## Loading required package: xml2## Warning: package 'xml2' was built under R version 4.0.3##
## Attaching package: 'rvest'## The following object is masked from 'package:purrr':
##
## pluck## The following object is masked from 'package:readr':
##
## guess_encodinglibrary(stringr)
library(ggmap)## Warning: package 'ggmap' was built under R version 4.0.5## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.## Please cite ggmap if you use it! See citation("ggmap") for details.library(data.table)
library(sp)
library(rgdal)
library(KernSmooth)
library(raster)
library(ggplot2)
library(RColorBrewer)
library(sf)
library(dplyr)
library(downloader)
library(tigris)
library(leaflet)
sf2 <- read_sf("/Users/john/Downloads/Verdicos/Exports/WellReportingPoint.shp")
#sf2
target <- c("Producing", "DUC")
target2 <- c("WOLFCAMP", "SPRABERRY")
sf3 <- sf2 %>% filter(WellStatus %in% target) %>% filter(CountyName == 'Midland'|CountyName == 'Martin'|CountyName == 'Andrews'|CountyName == 'Howard'|CountyName == 'Glasscock'|CountyName == 'Reagan'|CountyName == 'Upton'|CountyName == 'Ector') %>% filter(Reservoir %in% target2)
counties <- tigris::counties(state = 'TX')
counties <- st_as_sf(counties)
#RRMIV <- read_sf("/Users/john/Downloads/Verdicos/Exports/RockRiverIV.shp")
#RRMIV <- RRMIV %>% st_transform(crs = 4326)
sf3 <- sf3 %>% st_transform(crs = 4326)
counties <- counties %>% st_transform(crs = 4326)
counties <- counties %>% filter(NAME == 'Midland'|NAME == 'Martin'|NAME == 'Andrews'|NAME == 'Howard'|NAME == 'Glasscock'|NAME == 'Reagan'|NAME == 'Upton'|NAME == 'Ector')
##WOLFCAMP ANALYSIS
dat<- data.table(sf3 %>% filter(Reservoir == 'WOLFCAMP'))
kde <- bkde2D(dat[ , list(POINT_X, POINT_Y)],
bandwidth=c(.01, .01), gridsize = c(1000,1000))
#Convert kde to a Raster for plot
KernelDensityRaster <- raster(list(x=kde$x1 ,y=kde$x2 ,z = kde$fhat))
KernelDensityRaster@data@values[which(KernelDensityRaster@data@values < 1)] <- NA
#Color scale for the raster object (11 bins)
palRaster <- colorBin(rev(brewer.pal(n=11, name='Spectral')),bins = 10,
domain = KernelDensityRaster@data@values, na.color = "transparent")
#Create the leaflet Plot
leaflet(options = leafletOptions(zoomControl = FALSE)) %>% addTiles() %>%
addRasterImage(KernelDensityRaster,
colors = palRaster,
opacity = .8) %>%
addLegend(pal = palRaster,
values = KernelDensityRaster@data@values,
title = "Well Density<br>Wolfcamp")%>%
addPolygons(data = counties, fillColor = 'transparent', color = 'black', weight = 2) #%>%## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in CRS definition## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in CRS definition #addPolygons(data = RRMIV, fillColor = 'transparent', color = 'orange', weight = 3)
##SPRABERRY ANALYSIS
dat<- data.table(sf3 %>% filter(Reservoir == 'SPRABERRY'))
kde <- bkde2D(dat[ , list(POINT_X, POINT_Y)],
bandwidth=c(.01, .01), gridsize = c(1000,1000))
#Convert kde to a Raster for plot
KernelDensityRaster <- raster(list(x=kde$x1 ,y=kde$x2 ,z = kde$fhat))
KernelDensityRaster@data@values[which(KernelDensityRaster@data@values < 1)] <- NA
#Color scale for the raster object (11 bins)
palRaster <- colorBin(rev(brewer.pal(n=11, name='Spectral')),bins = 10,
domain = KernelDensityRaster@data@values, na.color = "transparent")
#Create the leaflet Plot
leaflet(options = leafletOptions(zoomControl = FALSE)) %>% addTiles() %>%
addRasterImage(KernelDensityRaster,
colors = palRaster,
opacity = .8) %>%
addLegend(pal = palRaster,
values = KernelDensityRaster@data@values,
title = "Well Density<br>Spraberry")%>%
addPolygons(data = counties, fillColor = 'transparent', color = 'black', weight = 2) #>% addPolygons(data = RRMIV, fillColor = 'transparent', color = 'orange', weight = 3)## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in CRS definition## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in CRS definition##Delaware Side Research
sf2 <- read_sf("/Users/john/Downloads/Verdicos/Exports/WellReportingPoint.shp")
sf2## Simple feature collection with 483788 features and 38 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: -112.7188 ymin: 0 xmax: 0 ymax: 48.99818
## geographic CRS: WGS 84
## # A tibble: 483,788 x 39
## API DLink RegNo LeaseName WellName WellNumber WellNameNo CountyName
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 0504~ "htt~ 0504~ GEISKE GEISKE 42B-26-692 GEISKE #4~ Garfield
## 2 0504~ "htt~ 0504~ HENRY HENRY PA 414-3 HENRY #PA~ Garfield
## 3 0504~ "htt~ 0504~ SAVAGE SAVAGE RWF 544-26 SAVAGE #R~ Garfield
## 4 0504~ "htt~ 0504~ MCPHERSON MCPHERS~ A12 MCPHERSON~ Garfield
## 5 0504~ "htt~ 0504~ AP AP 321-18-695 AP #321-1~ Garfield
## 6 0504~ "htt~ 0504~ WILLIAMS WILLIAMS GM 421-32 WILLIAMS ~ Garfield
## 7 0504~ "htt~ 0504~ WILLIAMS WILLIAMS GM 311-32 WILLIAMS ~ Garfield
## 8 0504~ "htt~ 0504~ ROBINSON ROBINSON C12 ROBINSON ~ Garfield
## 9 0504~ "htt~ 0504~ MCPHERSON MCPHERS~ A9 MCPHERSON~ Garfield
## 10 0504~ "htt~ 0504~ AP AP 322-18-695 AP #322-1~ Garfield
## # ... with 483,778 more rows, and 31 more variables: District <chr>,
## # StateAbbr <chr>, Operator <chr>, Field <chr>, Reservoir <chr>,
## # WellStatus <chr>, WellType <chr>, BoreDir <chr>, PmtDepth <chr>,
## # PmtAppr <dbl>, PmtApprDt <date>, PmtExp <dbl>, PmtExpDt <date>,
## # OilGathr <chr>, GasGathr <chr>, CumOil <dbl>, CumGas <dbl>,
## # FrstProdDt <date>, FrstProd <dbl>, LastProdDt <date>, LastProd <dbl>,
## # MeasDepth <chr>, VertDepth <chr>, SpudDt <date>, Spud <dbl>,
## # FrstCompDt <date>, FrstComp <dbl>, LastUpFrm <chr>, POINT_X <dbl>,
## # POINT_Y <dbl>, geometry <POINT [°]>target <- c("Producing", "DUC")
target2 <- c("WOLFCAMP", "BONE SPRING")
sf3 <- sf2 %>% filter(WellStatus %in% target) %>% filter(CountyName == 'Loving'|CountyName == 'Ward'|CountyName == 'Winkler'|CountyName == 'Pecos') %>% filter(Reservoir %in% target2)
counties <- tigris::counties(state = 'TX')
counties <- st_as_sf(counties)
#RRMIV <- read_sf("/Users/john/Downloads/Verdicos/Exports/RockRiverIV.shp")
#RRMIV <- RRMIV %>% st_transform(crs = 4326)
sf3 <- sf3 %>% st_transform(crs = 4326)
counties <- counties %>% st_transform(crs = 4326)
counties <- counties %>% filter(NAME == 'Loving'|NAME == 'Ward'|NAME == 'Winkler'|NAME == 'Pecos')
##WOLFCAMP ANALYSIS
dat<- data.table(sf3 %>% filter(Reservoir == 'WOLFCAMP'))
kde <- bkde2D(dat[ , list(POINT_X, POINT_Y)],
bandwidth=c(.01, .01), gridsize = c(1000,1000))## Warning in bkde2D(dat[, list(POINT_X, POINT_Y)], bandwidth = c(0.01, 0.01), :
## Binning grid too coarse for current (small) bandwidth: consider increasing
## 'gridsize'#Convert kde to a Raster for plot
KernelDensityRaster <- raster(list(x=kde$x1 ,y=kde$x2 ,z = kde$fhat))
KernelDensityRaster@data@values[which(KernelDensityRaster@data@values < 1)] <- NA
#Color scale for the raster object (11 bins)
palRaster <- colorBin(rev(brewer.pal(n=11, name='Spectral')),bins = 10,
domain = KernelDensityRaster@data@values, na.color = "transparent")
#Create the leaflet Plot
leaflet(options = leafletOptions(zoomControl = FALSE)) %>% addTiles() %>%
addRasterImage(KernelDensityRaster,
colors = palRaster,
opacity = .8) %>%
addLegend(pal = palRaster,
values = KernelDensityRaster@data@values,
title = "Well Density<br>Wolfcamp")%>%
addPolygons(data = counties, fillColor = 'transparent', color = 'black', weight = 2) #%>%## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in CRS definition## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in CRS definition # addPolygons(data = RRMIV, fillColor = 'transparent', color = 'orange', weight = 3)
##SPRABERRY ANALYSIS
dat<- data.table(sf3 %>% filter(Reservoir == 'BONE SPRING'))
kde <- bkde2D(dat[ , list(POINT_X, POINT_Y)],
bandwidth=c(.01, .01), gridsize = c(1000,1000))
#Convert kde to a Raster for plot
KernelDensityRaster <- raster(list(x=kde$x1 ,y=kde$x2 ,z = kde$fhat))
KernelDensityRaster@data@values[which(KernelDensityRaster@data@values < 1)] <- NA
#Color scale for the raster object (11 bins)
palRaster <- colorBin(rev(brewer.pal(n=11, name='Spectral')),bins = 10,
domain = KernelDensityRaster@data@values, na.color = "transparent")
#Create the leaflet Plot
leaflet(options = leafletOptions(zoomControl = FALSE)) %>% addTiles() %>%
addRasterImage(KernelDensityRaster,
colors = palRaster,
opacity = .8) %>%
addLegend(pal = palRaster,
values = KernelDensityRaster@data@values,
title = "Well Density<br>Spraberry")%>%
addPolygons(data = counties, fillColor = 'transparent', color = 'black', weight = 2)#%>%## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in CRS definition## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in CRS definition #addPolygons(data = RRMIV, fillColor = 'transparent', color = 'orange', weight = 3)