I converted a DWG file to an ESRI Shapefile using FME Workbench. Now I wish to create an elevation raster. I'm not sure how to extract the elevation info.


# https://www.dropbox.com/s/n3nvssl3swa8ryn/Massings.shp?dl=0
shapename <- read_sf('Massings.shp')
plot(shapename, axes=TRUE)

# Data has no assigned CRS. Set it to what website says it is.
st_crs(shapename) <- 2019

# Now can transform CRS to lat/long OK
shp_4326 <- st_transform(shapename, crs=4326)
plot(shp_4326, axes=TRUE)

# Create empty raster
ext <- extent(c(xmin = -79.395, xmax = -79.385, ymin = 43.638, ymax = 43.644))
gridsize <- 5e-5 
r <- raster(ext, res = gridsize)

# Fill with elevation data from shapefile ???
# dwntwn_raster <- rasterize(shp_crs, r, field = ???)

Shapefile 24.8MB

  • FIrst, you need a .shx to go with that shapefile - this can be generated using the ogrinfo command line tool. Next, your data has lots of overlapping polygons, with different Z coordinates. How are you expecting to extract an elevation from that? – Spacedman Apr 18 at 16:22
  • The elevation data in the shapefile is in the Z dimension. Every vertex of every polygon in your data has a Z coordinate, and they don't have to be the same, so your polygons could be sloping or even non-planar. You can extract the coordinates for every vertex of every polygon with st_coordinates. You could then get the mean or max or min Z coordinate by computing on the coordinates over the polygon grouping variable. Note there's some empty geometries that will throw an error. But the multiple overlapping geometry question remains. – Spacedman Apr 18 at 16:24

Your Answer

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

Browse other questions tagged or ask your own question.