rcsvrasterterrasf-package

How to read a csv (raster) into R and intersect it with a shapefile?

I'm trying to put together a csv dataset(csv_1) and a polygon shapefile (cities_shp) in R. The csv is a raster of 100meters resolution, with a variable of interest named "variable_1" (Shown in the picture below as multiple points). The shapefile has polygons for each multiple cities across the world (shown in the picture below as a single polygon).

I'm trying to create an intersection of the csv file and the polygon shapefile, so I can identify the value of "variable_1" across the city. As you can see in the picture below, the csv points do not correspond geographically to the "cities_shp" locations.

#####################
# Inputs
#####################
## Libraries ##
library(terra)
library(sf)
library(dplyr)
library(mapview)
library(raster)

# Read City polygon
cities_shp = st_read("Input/.../GHS_STAT_UCDB2015MT_GLOBE_R2019A_V1_1.shp") 
cities_shp = st_make_valid(cities_shp)

# Read CSV file
csv_1 = read.csv("Input/.../csv_1.csv") 


#####################
# Graph 
#####################
csv_1_graph = st_as_sf(csv_1, coords = c("x", "y"), crs = 3395)
csv_1_graph = st_transform(csv_1_graph, crs = st_crs(cities_shp))
mapview(csv_1_graph) + mapview(cities_shp)

enter image description here

What I did next was the following steps:

  1. to try to rasterize the "csv_1" file and then,
  2. to convert this new raster to a shapefile,
  3. to use the shapefile from step 2, so I can intersect it with the "cities_shp".

I only was able to achieve the 2nd step, because I wasn't able to visualize this new shapefile using the appropriate coordinates. With any of the 2 options I use (below), I don't get how to correct the CRS so I can intersect it with the "cities_shp". Because when I try it, I get an error with any of the options I show below.

Does anyone know how to read this csv file and then intersect it with the shapefile?

Note: I'm showing you just one city in this case, but I have to do it for 10,000 across the world. So, any help is really appreciated!

#####################
# 1st try
#####################
aux_2 <- try(rast(csv_1))
terra::plot(aux_2)
terra::plot(aux_2$variable_1)
aux_3 <- raster(aux_2)
mapview(aux_3,zcol = "variable_1",legend = T)
# Warning message:
#   In rasterCheckAdjustProjection(x, method) :
#   supplied layer has no projection information and is shown without background map
aux_2_polygons = as.polygons(aux_2)   |> st_as_sf()
class(aux_2_polygons)
new_df <- aux_2_polygons %>% st_as_sf(sf_column_name = "geometry", crs = 4326)
mapview(new_df)
#writing: substituting ENGCRS["Undefined Cartesian SRS with unknown unit"] for missing CRS

#####################
# 2nd try
#####################
geoName <- csv_1
geoName_raster <- raster(xmn = min(geoName$x), xmx = max(geoName$x), 
                         ymn = min(geoName$y), ymx = max(geoName$y), 
                         res = 100,crs = "+proj=longlat +datum=WGS84")
geoName_rasterize <- rasterize(geoName[, c('x', 'y')], geoName_raster, 
                               geoName[, 'variable_1'], fun = sum)
plot(geoName_rasterize)
points_matrix <- rasterToPoints(geoName_rasterize)
points_df <- as.data.frame(points_matrix) 
points_sf <- st_as_sf(points_df, coords = c("x", "y"), crs = "+proj=longlat +datum=WGS84")
points_sf = st_transform(points_sf, crs = 4326) #4326
plot(points_sf)
mapview(points_sf)
Solution

  • You assigned the wrong CRS to your csv_1 data. It appears that the correct CRS is WGS84 Pseudo Mercator/EPSG:3857.

    However, without knowing of the source of your csv, or how it was created, this may not be correct. You will need to ask the creator of the data to confirm whether it is EPSG:3857 or not.

    Either way:

    library(sf)
library(mapview)

# Assign EPSG:3857 CRS
csv_1_graph <- st_as_sf(csv_1, coords = c("x", "y"), crs = 3857)

# Project to st_crs(cities_shp)
csv_1_graph <- st_transform(csv_1_graph, st_crs(cities_shp))

mapview(csv_1_graph) +
  mapview(cities_shp)

    1