Подсчет количества точек в растровой сетке

Question

Я пытаюсь, чтобы растровая сетка распознала количество точек в ней и сосчитала их. Я создал растр, и у меня есть точки GPS. Я почти уверен, что и сетка, и точки, которые у меня есть, находятся в одном и том же CRS, но когда я пытаюсь построить их вместе, точки не появляются. Это то, что я сделал

Я создал растр и у меня есть точки GPS. Я почти уверен, что и сетка, и точки, которые у меня есть, находятся в одном и том же CRS, но когда я пытаюсь построить их вместе, точки не появляются. Это то, что я сделал

### 1. Create Raster###
#set the origin of raster
ori <- SpatialPoints(cbind(113.58353,-25.80607), proj4string = CRS("+init=epsg:4326"))

#convert the projection of ori and use EPSG: 3857 (Spherical Mercator)
ori_t <- spTransform(ori, CRSobj = CRS("+init=epsg:32749"))

#Determine the extent of the grid
# The origin has been rounded to the nearest 100
x_ori <- round(coordinates(ori_t)[1, 1]/100) * 100
y_ori <- round(coordinates(ori_t)[1, 2]/100) * 100

# Define how many cells for x and y axis
x_cell <- 13
y_cell <- 13

# Define the resolution to be 2000 meters
cell_size <- 2000

# Create the extent
ext <- extent(x_ori, x_ori + (x_cell * cell_size), 
              y_ori, y_ori + (y_cell * cell_size))
# Initialize a raster layer
ras <- raster(ext)

# Set the resolution to be
res(ras) <- c(cell_size, cell_size)
ras[] <- 0

# Project the raster
projection(ras) <- CRS("+init=epsg:32749")

### 2. Plotting the coordinates in the same CRS###
GPS <- read.csv("wow.1994.2017.gps.csv", sep=";")
# example of what the data looks like
head(GPS)
  project       Date alliance letterCode longitude  latitude
1       1 15.05.1994        5        AJA  113.6651 -25.71582
2       1 26.08.1994        5        AJA  113.6058 -25.66950
3       1 09.09.1994        5        AJA  113.6858 -25.70700
4       1 04.08.1995        5        AJA  113.6126 -25.68989
5       1 04.08.1995        5        AJA  113.6070 -25.68986
6       1 12.08.1995        5        AJA  113.6345 -25.64736

xy_GPS = GPS[c("longitude", "latitude")]
coordinates(xy_GPS)=c("longitude","latitude")
GPS_sp<-SpatialPointsDataFrame(xy_GPS, GPS) 


### 3. Tried to visualize them together ###
tab <- table(cellFromXY(ras, GPS_sp))
ras[as.numeric(names(tab))] <- tab
plot(ras)
points(GPS_sp, pch=20)

Answer 1

Ваш растр в UTM с extent(ras) = c (759000, 785000, 7143200, 7169200), в то время как ваши GPS-точки находятся в широте. R не будет автоматически проецировать очки для вас. Вам нужно присвоить своим GPS-точкам CRS, а затем спроецировать точки в тот же CRS, что и ras:

xy_GPS = GPS[,c("longitude", "latitude")]
coordinates(xy_GPS) = ~longitude+latitude
proj4string(xy_GPS) = CRS("+init=epsg:4326")
xy_GPS = spTransform(xy_GPS, crs(ras))

plot(ras)
points(xy_GPS)