2

I am currently trying to read in a catalog of las files (statewide 10TB) and extract around 500 plots from a shapefile of central points. The shapefile reads in correctly, and I have checked it in ArcGIS to be sure. However, when I attempt to display the catalog of las tiles, they appear disjointed and fail to line up with any of the points from the shapefile (which should fall well within the spatial extent). It appears that the northing may be correct, but the easting could be off. I attempted to change the projection of the las files to match the shapefile, but this did not remedy the issue.

catalog:

class       : LAScatalog
extent      : 70000 , 3195000 , 940000 , 2380000 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
area        : 997518.4 kunits²
points      : 379.02 billion points
density     : 0.4 points/units²
num. files  : 41353 

Shapefile:

Object of class SpatialPointsDataFrame

Coordinates:
                min       max
coords.x1 -87.97372 -84.84951
coords.x2  37.80275  41.72992
Is projected: FALSE 
proj4string :
[+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
Number of points: 530
Data attributes:
      PLOT          LAT             LON        
 1005   :  1   Min.   :37.80   Min.   :-87.97  
 1024   :  1   1st Qu.:38.25   1st Qu.:-86.71  
 1031   :  1   Median :38.70   Median :-86.34  
 1035   :  1   Mean   :38.82   Mean   :-86.39  
 1042   :  1   3rd Qu.:39.22   3rd Qu.:-86.21  
 1046   :  1   Max.   :41.73   Max.   :-84.85  
 (Other):524   

Here is the code I used.

shp <- readOGR("D:/All LAS/Shapefile/plot4.shp")  
ctg <- catalog(folder = "D:/All LAS")  
plot(ctg)  
EABplots <- lasclip(ctg, shp, radius = 7.3)
0

1 Answer 1

4

Your points in your shapefile are recorded with longitude latitude.

                min       max
coords.x1 -87.97372 -84.84951
coords.x2  37.80275  41.72992

This is also confirmed by the CRS and the attribute named LON and LAT

+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

However your point clouds, is projected

extent : 70000, 3195000, 940000, 2380000 (xmin, xmax, ymin, ymax)

The coordinates definitively do no match. And the CRS of the LAScatalog is wrong. It was probably NA and I guess you set it manually and you set a wrong CRS to your point cloud.

You must reproject your shapefile in the coordinate system of your point cloud (assuming that all the files have the same CRS). If the dataset does not have a CRS recorded you must find it.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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