The code is designed to analyse point locations along each polyline.
It works when ttl_ln is small.
But when the whole polyline features are loaded (more than 300,000 polylines), Rgui will closed automatically.
I did some test by adding a bool test_tag. When test_tag is true, the Rgui closing happens, when false, then the 300,000 can be processed by only using subset() and repeating matrix creation in the loop. So it may not be due to using subset() or matrix().
Shut down comes when using scale_to_grid() defined in the Rcpp.
shp_file = 'some_ply.shp'
shp_1 = sf::st_read(shp_file)
shp_2 = sf::st_coordinates(shp_1)
ttl_ln = nrow(shp_1)
for (ln in c(1:ttl_ln)){
shp_3 = subset(shp_2, shp_2[,3] == ln)
# field 1 is coordinate X
# field 2 is coordinate Y
# field 3 is Line ID created by sf package
em_mat = matrix(0, nrow = 200, ncol=200)
if(test_tag){
# manipulations on points are followed
# a groupe of functions defined in a Rcpp file
# are used. They can be successfully compiled.
new_mat = scale_to_grid(shp_3, em_mat, t1, t2, t3, t4)
# t1, t2, t3, t4 are numbers with double type in Rcpp file.
}
}
Following is the code in Rcpp.
// [[Rcpp::export]]
NumericMatrix scale_to_grid(NumericMatrix shp,
NumericMatrix mat,
double o_x, double o_y,
double u_x, double u_y){
// mat : the matrix of grid
int ptr;
int nbr = shp.nrow();
int ind_x, ind_y;
double slope, new_x, new_y;
// go through each point, expect the last one (i.e. ptr < nbr -1)
for(ptr=0;ptr<nbr-1;ptr++){
ind_x = trunc((shp(ptr,0)-o_x)/u_x);
ind_y = trunc((shp(ptr,1)-o_y)/u_y);
mat(ind_x, ind_y)++;
if(abs(shp(ptr+1,0) - shp(ptr,0)) > u_x){
// the difference in x value is greater than the grid unit width
// there must be grid to fill with value '1'
slope = (shp(ptr+1,1) - shp(ptr,1)) / (shp(ptr+1,0) - shp(ptr,0));
// initialise
new_x = shp(ptr,0);
new_y = shp(ptr,1);
while(abs(new_x - shp(ptr+1,0)) > u_x){
if(shp(ptr+1,0) - shp(ptr,0) > 0){
new_x = new_x + u_x;
new_y = new_y + u_x * slope;
}else{
new_x = new_x - u_x;
new_y = new_y - u_x * slope;
}
ind_x = trunc((new_x-o_x)/u_x);
ind_y = trunc((new_y-o_y)/u_y);
mat(ind_x, ind_y)++;
}
}
}
return mat;}