tylermorganwall / MusaMasterclass

My masterclass hosted by the Penn MUSA department on 3D mapping and visualization in R
https://upenn.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=ae6642b2-c2be-4ad0-846c-aafb00e00129
127 stars 51 forks source link

EPSG code use when cropping #1

Closed r-leyshon closed 4 years ago

r-leyshon commented 4 years ago

Hi Tyler,

Great class!

My issue is with the custom function lat_long_to_other:

lat_long_to_other = function(lat,long, epsg_code) {
  data = data.frame(long=long, lat=lat)
  coordinates(data) <- ~ long+lat
  #Input--lat/long
  proj4string(data) <- CRS("+init=epsg:XXXX")
  #Convert to coordinate system specified by EPSG code
  xy = data.frame(spTransform(data, CRS(paste0("+init=epsg:", epsg_code))))
  colnames(xy) = c("x","y")
  return(unlist(xy))
}

bottomleft = lat_long_to_other(51.4452463, -3.1812256, 1314)
topright = lat_long_to_other(51.4667095, -3.163987, 1314)

As you can see, I've found the espg code for my chosen location, but get an error when attempting to run. I Believe the error is probably caused by not knowing the correct value to pass to CRS("+init=epsg:XXXX")

I can see that this is a different value from the Miami-Dade espg used in the masterclass markdown. Do you have any tips for establishing the correct value for this argument?

Here is the error: Error in CRS("+init=epsg:1314") : no options found in 'init' file

tylermorganwall commented 4 years ago

Thanks! The lat_long_to_other() function is meant for the simple case of converting WGS84 lat/long coordinates (which is is EPSG code 4326) to some other coordinate system system with an EPSG code, but it's not designed to handle more complex coordinate conversions.

I recommend reading the following resources to see if they help with your problem.

Geocomputation in R: https://geocompr.robinlovelace.net/reproj-geo-data.html

Intro to GIS and Spatial Analysis: https://mgimond.github.io/Spatial/coordinate-systems-in-r.html

tylermorganwall commented 4 years ago

Question: Is that verbatim the code you're using? The code in the masterclass is the following:

lat_long_to_other = function(lat,long, epsg_code) {
  data = data.frame(long=long, lat=lat)
  coordinates(data) <- ~ long+lat
  #Input--lat/long
  proj4string(data) <- CRS("+init=epsg:4326")
  #Convert to coordinate system specified by EPSG code
  xy = data.frame(spTransform(data, CRS(paste0("+init=epsg:", epsg_code))))
  colnames(xy) = c("x","y")
  return(unlist(xy))
}

As you can see, the EPSG code for lat/long data is entered in as 4326 and not XXXX. The EPSG code you specify in the function argument should be your destination CRS, and the 4326 number need not be altered.

r-leyshon commented 4 years ago

Good question - The lidar data I'm using has the following metadata: Spatial Reference System: OSGB36 epsg = 1314 I think this is probably the issue - knowing which code to input into the CRS call. I tried 1314 and it returned the above error. I've managed to avoid the issue by masking the map with a spatial polygon with the required extent. But it would be great to figure out this issue, as I see myself doing a lot of this type of work in the future.

tylermorganwall commented 4 years ago

If you're using lat/long data, you shouldn't need to change the proj4string(data) <- CRS("+init=epsg:4326") line to anything else (unless you're using a different reference ellipsoid than WGS84). Your destination CRS should be specified in the epsg_code argument in the function call itself, and you shouldn't need to change the inside of the function.

If you're converting from something other than lat/long to your lidar coordinate system, then you will have to do the conversion yourself (this function will not work).