2

I'm working with a global shapefile in R that is latlong projected. I'd like to reproject it into the appropriate UTM projection for individual zones. But, I'd like to chop it up into those zones before reprojection, instead of having multiple global reprojections into UTM most of which have meaningless coordinates. Also, computational efficiency. So, how would I split a latlong shapefile into areas corresponding to UTM zones before the reprojection (and a little overhang is fine)?

8
  • I am sure you have good reasons for this, but just in case, the world really doesn't need more data in UTM. Is there not a better choice of projection/s for your application? – mdsumner Feb 4 '17 at 6:29
  • I agree. I'm working with Landsat data to see if part of a scene falls along a coastline or not. I have a global coastline shapefile, in latlong, and all scenes are in UTM. Computationally, it's less expensive to reproject the coastline into UTM zones and then use that to use gIntersect or somesuch on rather than reproject each of the thousands of rasters I'm working with. – jebyrnes Feb 5 '17 at 16:18
  • They are in UTM? Gosh. Is that standard delivery for Landsat or some sub project? I had no idea, keen to explore can you point to some sources? Usually sat data has some rawer form or is swath data or some projection masked by longlat arrays. – mdsumner Feb 5 '17 at 21:25
  • Yeah, that's the standard. Like you, i assumed otherwise - WGS84, and then...had some problems that made me take a second look. That's been...an experience. – jebyrnes Feb 6 '17 at 19:53
  • So really you want a map source that will always plot this in context correctly? This slicing gores approach won't work all the time. Still personally seeking best ways to do this kind of thing for myself. – mdsumner Feb 6 '17 at 21:52
3

You can use crop from raster to crop to a defined extent, in this case your UTM zone approximate boundary.

So I have a world map shapefile:

> w
class       : SpatialPolygons 
features    : 1 
extent      : -179.9997, 180, -89.9999, 83.2479  (xmin, xmax, ymin, ymax)
coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 

And some xy coordinates of a box to crop to:

> xy
$x
[1] -4.386177 28.181688

$y
[1]  78.54324 -86.03500

You might want to tweak the latitude limits depending on how near the poles your world shapefile goes. Set the longitude (x) to your UTM zone.

Then crop.

> library(raster)
> wc = crop(w, extent(SpatialPoints(data.frame(xy))))
> plot(w)
> plot(wc, add=TRUE,col="red",lwd=2)

enter image description here

Note you better not have any odd geometry problems in your world shapefile...

2

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.