3

I'm trying to clip a catalog using lidR (v2.0.1) package. I've tried with lasclipCircle and lasclip. It gives me the error: no point in ROI, but when I'm viewing the shapefile using mapview (I've exported the X and Y coordinates from this shapefile) and the catalog, they are at the same place.

Here's my code.

page1 <- catalog("O:/LiDAR/NB_LiDAR_16/page1")
Buffer_16 <- readOGR(dsn = "O:/2018_Regen/Derived GIS", layer = "Buffer_16")
X_16 <- Buffer_16@data[["DDLONG"]]
Y_16 <- Buffer_16@data[["DDLAT"]]
page1_points <- lasclipCircle(page1, xcenter = X_16,ycenter = Y_16, radius = 10)

or

page1_points <- lasclip(page1, Buffer_16)

I get the following error:

Error : No point found for within region of interest 2.

I got the same message for ROI 1 and all the other ROI as well.

Summary of Buffer_16

Object of class SpatialPolygonsDataFrame
Coordinates:
      min     max
x 2457220 2613228
y 7490508 7662516
Is projected: TRUE 
proj4string :
[+init=epsg:2953 +proj=sterea +lat_0=46.5 +lon_0=-66.5 +k=0.999912 +x_0=2500000 +y_0=7500000 +ellps=GRS80
+towgs84=0,0,0,0,0,0,0 +units=m +no_defs]
Data attributes:
    OBJECTID         MAPNO         PLOTID          X             Y           DDLAT           DDLONG       SCHEDULE
 Min.   : 62.0   4721   : 11   23239  :  1   2541224: 11   7600512:  9   Min.   :46.41   Min.   :-67.07   1:179   
 1st Qu.:257.5   4921   : 11   26157  :  1   2539224: 10   7554512:  8   1st Qu.:47.09   1st Qu.:-66.46           
 Median :454.0   5230   : 10   26159  :  1   2543224:  9   7604512:  8   Median :47.35   Median :-66.14           
 Mean   :415.1   4425   :  9   26160  :  1   2505224:  8   7646512:  8   Mean   :47.34   Mean   :-66.12           
 3rd Qu.:578.5   4714   :  9   26369  :  1   2501224:  7   7602512:  7   3rd Qu.:47.61   3rd Qu.:-65.88           
 Max.   :790.0   5523   :  9   26370  :  1   2507224:  7   7648512:  7   Max.   :47.96   Max.   :-65.00           
                 (Other):120   (Other):173   (Other):127   (Other):132                                            
 ISPSP     OBJECTID_1   objectid_2 PHOTO_YEAR  LABEL     NUMBERING   Shape_Leng       Shape_Area       
 0:131   Min.   :5    Min.   :2    2016:179   2016:179   0:179     Min.   :749732   Min.   :1.876e+10  
 1: 48   1st Qu.:5    1st Qu.:2                                    1st Qu.:749732   1st Qu.:1.876e+10  
         Median :5    Median :2                                    Median :749732   Median :1.876e+10  
         Mean   :5    Mean   :2                                    Mean   :749732   Mean   :1.876e+10  
         3rd Qu.:5    3rd Qu.:2                                    3rd Qu.:749732   3rd Qu.:1.876e+10  
         Max.   :5    Max.   :2                                    Max.   :749732   Max.   :1.876e+10`

Summary of page1:

class       : LAScatalog
extent      : 2457000 , 2590000 , 7542000 , 7647000 (xmin, xmax, ymin, ymax)
coord. ref. : +init=epsg:2953 +proj=sterea +lat_0=46.5 +lon_0=-66.5 +k=0.999912 +x_0=2500000 +y_0=7500000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
area        : 43 km²
points      : 844.56 million points
density     : 19.6 points/m²
num. files  : 43 
Summary of the processing options:
  - Catalog will be processed by file respecting the original tiling pattern
  - Catalog will be processed using 1 core(s).
Summary of the output options:
  - Outputs will be returned in R objects.
Summary of the input options:
  - readLAS will be called with the following select option: * 
  - readLAS will be called with the following filter option:    `  
  • 2
    Could you please confirm that all of your data are in the same projected coordinate system? – Aaron Feb 27 at 19:33
  • Please provide more details e.g. a summary of Buffer_16 and a summary of page1. Aaron is right but it may also be because the region of interest 2 (plot 2) belongs in a region with no data (data removed by the provider or any other reason). Plot number 1 (at least) has been extracted so coordinates seem to match – JRR Feb 28 at 4:53
  • Btw, the message should not be an error, not even a warning. It should only be a message to highlight the fact that you've extracted a 0 points region. – JRR Feb 28 at 8:20
  • JRR: I got the same message for ROI 1 and all the other ROI. Here's the projected coordinate system of Buffer_16 : [+proj=sterea +lat_0=46.5 +lon_0=-66.5 +k=0.999912 +x_0=2500000 +y_0=7500000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs] And for page1: +init=epsg:2953 +proj=sterea +lat_0=46.5 +lon_0=-66.5 +k=0.999912 +x_0=2500000 +y_0=7500000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs – Stéphanie Landry Feb 28 at 13:33
  • Please show us a summary of your objects. – JRR Feb 28 at 20:39
1

When doing that:

X_16 <- Buffer_16@data[["DDLONG"]]
Y_16 <- Buffer_16@data[["DDLAT"]]
page1_points <- lasclipCircle(page1, xcenter = X_16,ycenter = Y_16, radius = 10)

you are trying to extract long/lat coordinates from projected coordinates. So basically you tried to extract a disc at (47, -66) from data at (2457000, 7490000). No data is the expected behavior. You must recompute DDLONG DDLAT in the good coordinate system. Maybe the attributes X and Y of your shapefile are what you actually need. But in your data they are factors that must be converted to numbers first.

X_16 <- as.numeric(as.character(Buffer_16@data[["X"]]))
Y_16 <- as.numeric(as.character(Buffer_16@data[["Y"]]))
page1_points <- lasclipCircle(page1, X_16, Y_16, 10)

If many of your coordinates don't actually belong into your tiles you will extract many empty point clouds. In this case filter the empty ones

not.empty <- Filter(Negate(is.empty), page1_points)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.