Hallie Sheikh
Hallie Sheikh

Reputation: 431

How to get a subset from a shapefile with the correct outline?

I have a Hindu Kush Himalayan region shapefile that can be found in http://rds.icimod.org/Home/DataDetail?metadataId=3924 , and I have lat long of Nepal that have taken from this website (http://rds.icimod.org/Home/DataDetail?metadataId=19590&searchlist=True).

The geographical extent of Nepal is Geographic Extent East: 88.19456,Geographic Extent West:80.0522,Geographic Extent North:30.42472,Geographic Extent South: 26.36836

Now I am trying to subset Nepal from Hindu Kush Himalayan shapefile. This is my code:

mountains<-readOGR("outline.shp") #hindukushhimalayanshapefile
sub <- crop(mountains, extent( 80.0522, 88.18456, 26.36836, 30.42472))
plot(sub)

But the plot sub (that is of Nepal) is not shown properly with a proper outline. There are straight lines showing in the top. How can I get a proper subset of Nepal with the right outline. Am I putting the extent wrong? Help would be appreciated

sub plot

Upvotes: 1

Views: 220

Answers (1)

drJones
drJones

Reputation: 1233

You cannot subset country boundaries from that shapefile as it does not contain information on country boundaries. You would have to use a shapefile such as this one which contains information on country boundaries: http://rds.icimod.org/Home/DataDetail?metadataId=1218, or subset using a separate Nepal boundary shapefile instead of an extent:

require(maptools)
require(rgdal)

hkh_shp=readOGR("/Downloads/data/outline.shp")

data(wrld_simpl) ##this shapefile is quite coarse, you could substitute another one
nepal_shp=wrld_simpl[which(wrld_simpl$NAME=="Nepal"),]

##CRS are similar but not identical, so need to transform
nepal_shp=spTransform(nepal_shp,crs(hkh_shp))

plot(hkh_shp)
lines(nepal_shp,col="red")

enter image description here

##crop
hkh_sub_shp=crop(hkh_shp,nepal_shp)
plot(hkh_sub_shp) ##note, will look better with higher resolution shapefile

enter image description here

Upvotes: 1

Related Questions