Sorry - here is essentially what I do. There are no error messages since it runs, but the result is incorrect because all points are shown as being a single point within a single polygon. The coordinates are different and that is likely the issue but I can’t figure out how to correct it. The bounding boxes show this:
shapefile_util@bbox
min max
x -374446.5 540037.5
y -604489.9 450020.7
v_comp@bbox
min max
lon -124.17578 -115.55391
lat 32.54582 40.80654library(tidyverse)## ── Attaching packages ────────────────────────────────────── tidyverse 1.2.1 ──## ✔ ggplot2 3.1.0 ✔ purrr 0.3.0
## ✔ tibble 2.0.1 ✔ dplyr 0.7.8
## ✔ tidyr 0.8.2 ✔ stringr 1.4.0
## ✔ readr 1.3.1 ✔ forcats 0.3.0## ── Conflicts ───────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()library(sp)
library(spatialEco)
library(maptools)## Checking rgeos availability: TRUE# set paths
dir_data <- "~/Dropbox/CA util" # change this for your machine
vouchers_csv <- file.path(dir_data, "vouchers_subset.csv")
utilities_shp <- file.path(
dir_data,
"UtilityShapefile/California_Electric_Utility_Service_Territory_.shp")
#import voucher data
vouchers_df <- read_csv(vouchers_csv)## Parsed with column specification:
## cols(
## .default = col_character(),
## ID = col_double(),
## PaymentDate = col_logical(),
## Amount = col_double(),
## ModelYear = col_double(),
## EngineYear = col_double(),
## IsActive = col_logical(),
## lon = col_double(),
## lat = col_double(),
## AssemblyDistrict = col_double(),
## SenateDistrict = col_double(),
## CensusTract = col_double()
## )## See spec(...) for full column specifications.sp, maptools functions#convert points to spatialpointsdataframe
vouchers_pt <- SpatialPointsDataFrame(
vouchers_df[, 15:16],
vouchers_df[, 1:22],
proj4string = CRS("+proj=longlat +datum=WGS84"))
vouchers_pt@bbox## min max
## lon -122.78356 -115.56913
## lat 32.72263 40.57319#Import shapefile of counties
utilities_ply <- readShapeSpatial(
utilities_shp,
proj4string = CRS( "+proj=longlat +datum=WGS84"))## Warning: readShapeSpatial is deprecated; use rgdal::readOGR or sf::st_read## Warning: readShapePoly is deprecated; use rgdal::readOGR or sf::st_readutilities_ply@bbox## min max
## x -374446.5 540037.5
## y -604489.9 450020.7#Match data with utility polygons
vouchers_in_utilities_pt <- point.in.poly(vouchers_pt, utilities_ply)## although coordinates are longitude/latitude, st_intersects assumes that they are planarhead(vouchers_in_utilities_pt)## coordinates ID VoucherID CreatedOn PaymentDate
## 1 (-122.4935, 37.95006) 11658 14-8262-ARB 11/15/2018 8:00 NA
## 2 (-122.4935, 37.95006) 11659 14-8263-ARB 11/15/2018 8:00 NA
## 3 (-122.4935, 37.95006) 11660 14-8264-ARB 11/15/2018 8:00 NA
## 4 (-122.4935, 37.95006) 11661 14-8265-ARB 11/15/2018 8:00 NA
## 5 (-122.4935, 37.95006) 11662 14-8266-ARB 11/15/2018 8:00 NA
## 6 (-122.4935, 37.95006) 11663 14-8267-ARB 11/15/2018 8:00 NA
## Amount Dlr_Company Manufacturer VehicleType ModelYear EngineYear
## 1 80000 Chanje Energy, Inc. Chanje Delivery 2018 2018
## 2 80000 Chanje Energy, Inc. Chanje Delivery 2018 2018
## 3 80000 Chanje Energy, Inc. Chanje Delivery 2018 2018
## 4 80000 Chanje Energy, Inc. Chanje Delivery 2018 2018
## 5 80000 Chanje Energy, Inc. Chanje Delivery 2018 2018
## 6 80000 Chanje Energy, Inc. Chanje Delivery 2018 2018
## Description EngineType
## 1 Chanje Truck V8100, 4x2, 59.8 sq ft. frontal area Electric
## 2 Chanje Truck V8100, 4x2, 59.8 sq ft. frontal area Electric
## 3 Chanje Truck V8100, 4x2, 59.8 sq ft. frontal area Electric
## 4 Chanje Truck V8100, 4x2, 59.8 sq ft. frontal area Electric
## 5 Chanje Truck V8100, 4x2, 59.8 sq ft. frontal area Electric
## 6 Chanje Truck V8100, 4x2, 59.8 sq ft. frontal area Electric
## GrossVehicleWeight IsActive lon lat AssemblyDistrict
## 1 16,535 TRUE -122.4935 37.95006 10
## 2 16,535 TRUE -122.4935 37.95006 10
## 3 16,535 TRUE -122.4935 37.95006 10
## 4 16,535 TRUE -122.4935 37.95006 10
## 5 16,535 TRUE -122.4935 37.95006 10
## 6 16,535 TRUE -122.4935 37.95006 10
## AssemblyDistName SenateDistrict SenateDistName County CensusTract
## 1 MARIN 2 NORCO Marin County 6041112202
## 2 MARIN 2 NORCO Marin County 6041112202
## 3 MARIN 2 NORCO Marin County 6041112202
## 4 MARIN 2 NORCO Marin County 6041112202
## 5 MARIN 2 NORCO Marin County 6041112202
## 6 MARIN 2 NORCO Marin County 6041112202
## OBJECTID Shape_Leng Name Acronym Category
## 1 26 0.2374779 Pacific Gas & Electric Company PG&E IOU
## 2 26 0.2374779 Pacific Gas & Electric Company PG&E IOU
## 3 26 0.2374779 Pacific Gas & Electric Company PG&E IOU
## 4 26 0.2374779 Pacific Gas & Electric Company PG&E IOU
## 5 26 0.2374779 Pacific Gas & Electric Company PG&E IOU
## 6 26 0.2374779 Pacific Gas & Electric Company PG&E IOU
## Shape__Are Shape__Len
## 1 172701145927 6364469
## 2 172701145927 6364469
## 3 172701145927 6364469
## 4 172701145927 6364469
## 5 172701145927 6364469
## 6 172701145927 6364469maptools::readShapeSpatial()Per the documentation accessed from the R Console via ?readShapeSpatial or online at readShapeSpatial function | R Documentation, you’re using the deprecated function maptools::readShapeSpatial() that does not read in the spatial projection information. Per documentation you should use rgdal::readOGR() or better yet sf::st_read():
The use of this function is deprecated and it is not being maintained. Use rgdal::readOGR() or sf::st_read() instead - both of these read the coordinate reference system from the input file, while this deprecated function does not. For writing, use rgdal::writeOGR() or sf::st_write() instead.
Most significantly, this function doesn’t properly read the coordinate reference system and manually setting it doesn’t transform it.
rgdal::ogrInfo(utilities_shp)## Source: "/Users/bbest/Dropbox/CA util/UtilityShapefile/California_Electric_Utility_Service_Territory_.shp", layer: "California_Electric_Utility_Service_Territory_"
## Driver: ESRI Shapefile; number of rows: 59
## Feature type: wkbPolygon with 2 dimensions
## Extent: (-374446.5 -604489.9) - (540037.5 450020.7)
## CRS: +proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs
## LDID: 0
## Number of fields: 7
## name type length typeName
## 1 OBJECTID 12 10 Integer64
## 2 Shape_Leng 2 24 Real
## 3 Name 4 80 String
## 4 Acronym 4 80 String
## 5 Category 4 80 String
## 6 Shape__Are 2 24 Real
## 7 Shape__Len 2 24 Realsf functionslibrary(sf)## Linking to GEOS 3.6.1, GDAL 2.1.3, PROJ 4.9.3vouchers_pt <- st_as_sf(vouchers_df, coords = c("lon", "lat"), crs = 4326)
vouchers_pt## Simple feature collection with 500 features and 20 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: -122.7836 ymin: 32.72263 xmax: -115.5691 ymax: 40.57319
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
## # A tibble: 500 x 21
## ID VoucherID CreatedOn PaymentDate Amount Dlr_Company Manufacturer
## <dbl> <chr> <chr> <lgl> <dbl> <chr> <chr>
## 1 11658 14-8262-… 11/15/20… NA 80000 Chanje Ene… Chanje
## 2 11659 14-8263-… 11/15/20… NA 80000 Chanje Ene… Chanje
## 3 11660 14-8264-… 11/15/20… NA 80000 Chanje Ene… Chanje
## 4 11661 14-8265-… 11/15/20… NA 80000 Chanje Ene… Chanje
## 5 11662 14-8266-… 11/15/20… NA 80000 Chanje Ene… Chanje
## 6 11663 14-8267-… 11/15/20… NA 80000 Chanje Ene… Chanje
## 7 11664 14-8268-… 11/15/20… NA 80000 Chanje Ene… Chanje
## 8 11665 14-8269-… 11/15/20… NA 80000 Chanje Ene… Chanje
## 9 11666 14-8270-… 11/15/20… NA 80000 Chanje Ene… Chanje
## 10 11667 14-8271-… 11/15/20… NA 80000 Chanje Ene… Chanje
## # … with 490 more rows, and 14 more variables: VehicleType <chr>,
## # ModelYear <dbl>, EngineYear <dbl>, Description <chr>,
## # EngineType <chr>, GrossVehicleWeight <chr>, IsActive <lgl>,
## # AssemblyDistrict <dbl>, AssemblyDistName <chr>, SenateDistrict <dbl>,
## # SenateDistName <chr>, County <chr>, CensusTract <dbl>, geometry <POINT
## # [°]>#Import shapefile of counties
utilities_ply <- read_sf(utilities_shp)
utilities_ply## Simple feature collection with 59 features and 7 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -374446.5 ymin: -604489.9 xmax: 540037.5 ymax: 450020.7
## epsg (SRID): NA
## proj4string: +proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs
## # A tibble: 59 x 8
## OBJECTID Shape_Leng Name Acronym Category Shape__Are Shape__Len
## <dbl> <dbl> <chr> <chr> <chr> <dbl> <dbl>
## 1 1 1.36 Lass… <NA> POU 5.47e 9 617963.
## 2 2 1.36 Paci… <NA> IOU 2.50e10 992989.
## 3 3 1.36 Surp… <NA> Co-op 6.38e 9 508510.
## 4 4 1.36 Paci… <NA> Co-op 3.71e 9 437655.
## 5 5 1.78 Trin… <NA> POU 8.31e 9 607666.
## 6 6 1.36 Lass… <NA> Combo 4.63e 8 181591.
## 7 7 1.36 Plum… <NA> Co-op 3.67e 9 389560.
## 8 8 0 City… <NA> POU 1.24e 7 28641.
## 9 9 0 Bigg… <NA> POU 1.62e 6 10533.
## 10 10 0 Redd… <NA> POU 1.59e 8 146491.
## # … with 49 more rows, and 1 more variable: geometry <MULTIPOLYGON [m]>utilities_ply <- read_sf(utilities_shp) %>%
st_transform(crs = 4326)
utilities_ply## Simple feature collection with 59 features and 7 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -124.4152 ymin: 32.53434 xmax: -114.1312 ymax: 42.00952
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
## # A tibble: 59 x 8
## OBJECTID Shape_Leng Name Acronym Category Shape__Are Shape__Len
## <dbl> <dbl> <chr> <chr> <chr> <dbl> <dbl>
## 1 1 1.36 Lass… <NA> POU 5.47e 9 617963.
## 2 2 1.36 Paci… <NA> IOU 2.50e10 992989.
## 3 3 1.36 Surp… <NA> Co-op 6.38e 9 508510.
## 4 4 1.36 Paci… <NA> Co-op 3.71e 9 437655.
## 5 5 1.78 Trin… <NA> POU 8.31e 9 607666.
## 6 6 1.36 Lass… <NA> Combo 4.63e 8 181591.
## 7 7 1.36 Plum… <NA> Co-op 3.67e 9 389560.
## 8 8 0 City… <NA> POU 1.24e 7 28641.
## 9 9 0 Bigg… <NA> POU 1.62e 6 10533.
## 10 10 0 Redd… <NA> POU 1.59e 8 146491.
## # … with 49 more rows, and 1 more variable: geometry <MULTIPOLYGON [°]>vouchers_in_utilities_pt <- st_join(
vouchers_pt, utilities_ply, join = st_intersects)## although coordinates are longitude/latitude, st_intersects assumes that they are planarvouchers_in_utilities_pt## Simple feature collection with 500 features and 27 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: -122.7836 ymin: 32.72263 xmax: -115.5691 ymax: 40.57319
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
## First 10 features:
## ID VoucherID CreatedOn PaymentDate Amount
## 1 11658 14-8262-ARB 11/15/2018 8:00 NA 80000
## 2 11659 14-8263-ARB 11/15/2018 8:00 NA 80000
## 3 11660 14-8264-ARB 11/15/2018 8:00 NA 80000
## 4 11661 14-8265-ARB 11/15/2018 8:00 NA 80000
## 5 11662 14-8266-ARB 11/15/2018 8:00 NA 80000
## 6 11663 14-8267-ARB 11/15/2018 8:00 NA 80000
## 7 11664 14-8268-ARB 11/15/2018 8:00 NA 80000
## 8 11665 14-8269-ARB 11/15/2018 8:00 NA 80000
## 9 11666 14-8270-ARB 11/15/2018 8:00 NA 80000
## 10 11667 14-8271-ARB 11/15/2018 8:00 NA 80000
## Dlr_Company Manufacturer VehicleType ModelYear EngineYear
## 1 Chanje Energy, Inc. Chanje Delivery 2018 2018
## 2 Chanje Energy, Inc. Chanje Delivery 2018 2018
## 3 Chanje Energy, Inc. Chanje Delivery 2018 2018
## 4 Chanje Energy, Inc. Chanje Delivery 2018 2018
## 5 Chanje Energy, Inc. Chanje Delivery 2018 2018
## 6 Chanje Energy, Inc. Chanje Delivery 2018 2018
## 7 Chanje Energy, Inc. Chanje Delivery 2018 2018
## 8 Chanje Energy, Inc. Chanje Delivery 2018 2018
## 9 Chanje Energy, Inc. Chanje Delivery 2018 2018
## 10 Chanje Energy, Inc. Chanje Delivery 2018 2018
## Description EngineType
## 1 Chanje Truck V8100, 4x2, 59.8 sq ft. frontal area Electric
## 2 Chanje Truck V8100, 4x2, 59.8 sq ft. frontal area Electric
## 3 Chanje Truck V8100, 4x2, 59.8 sq ft. frontal area Electric
## 4 Chanje Truck V8100, 4x2, 59.8 sq ft. frontal area Electric
## 5 Chanje Truck V8100, 4x2, 59.8 sq ft. frontal area Electric
## 6 Chanje Truck V8100, 4x2, 59.8 sq ft. frontal area Electric
## 7 Chanje Truck V8100, 4x2, 59.8 sq ft. frontal area Electric
## 8 Chanje Truck V8100, 4x2, 59.8 sq ft. frontal area Electric
## 9 Chanje Truck V8100, 4x2, 59.8 sq ft. frontal area Electric
## 10 Chanje Truck V8100, 4x2, 59.8 sq ft. frontal area Electric
## GrossVehicleWeight IsActive AssemblyDistrict AssemblyDistName
## 1 16,535 TRUE 10 MARIN
## 2 16,535 TRUE 10 MARIN
## 3 16,535 TRUE 10 MARIN
## 4 16,535 TRUE 10 MARIN
## 5 16,535 TRUE 10 MARIN
## 6 16,535 TRUE 10 MARIN
## 7 16,535 TRUE 10 MARIN
## 8 16,535 TRUE 10 MARIN
## 9 16,535 TRUE 10 MARIN
## 10 16,535 TRUE 17 ESF
## SenateDistrict SenateDistName County CensusTract OBJECTID
## 1 2 NORCO Marin County 6041112202 26
## 2 2 NORCO Marin County 6041112202 26
## 3 2 NORCO Marin County 6041112202 26
## 4 2 NORCO Marin County 6041112202 26
## 5 2 NORCO Marin County 6041112202 26
## 6 2 NORCO Marin County 6041112202 26
## 7 2 NORCO Marin County 6041112202 26
## 8 2 NORCO Marin County 6041112202 26
## 9 2 NORCO Marin County 6041112202 26
## 10 11 SF San Francisco County 6075980900 26
## Shape_Leng Name Acronym Category Shape__Are
## 1 0.2374779 Pacific Gas & Electric Company PG&E IOU 172701145927
## 2 0.2374779 Pacific Gas & Electric Company PG&E IOU 172701145927
## 3 0.2374779 Pacific Gas & Electric Company PG&E IOU 172701145927
## 4 0.2374779 Pacific Gas & Electric Company PG&E IOU 172701145927
## 5 0.2374779 Pacific Gas & Electric Company PG&E IOU 172701145927
## 6 0.2374779 Pacific Gas & Electric Company PG&E IOU 172701145927
## 7 0.2374779 Pacific Gas & Electric Company PG&E IOU 172701145927
## 8 0.2374779 Pacific Gas & Electric Company PG&E IOU 172701145927
## 9 0.2374779 Pacific Gas & Electric Company PG&E IOU 172701145927
## 10 0.2374779 Pacific Gas & Electric Company PG&E IOU 172701145927
## Shape__Len geometry
## 1 6364469 POINT (-122.4935 37.95006)
## 2 6364469 POINT (-122.4935 37.95006)
## 3 6364469 POINT (-122.4935 37.95006)
## 4 6364469 POINT (-122.4935 37.95006)
## 5 6364469 POINT (-122.4935 37.95006)
## 6 6364469 POINT (-122.4935 37.95006)
## 7 6364469 POINT (-122.4935 37.95006)
## 8 6364469 POINT (-122.4935 37.95006)
## 9 6364469 POINT (-122.4935 37.95006)
## 10 6364469 POINT (-122.3947 37.74809)Read documentation on functions
Try making a reproducible example when enlisting help from others: https://reprex.tidyverse.org/