# Extracting a single layer from a SHP file

I have an ESRI ShapeFile (SHP) that I've loaded. It has a number of layers which I can see the names of by using:

Import["file.shp", "LayerNames"]


I get something like

{"Africover_Kenya_Spatial_Aggregation", "ke-landcover", "ke-forests", \ "ke-rangeland", "ke-agriculture", "ke-waterbodies", "ke-bareareas", \ "ke-urban", "ke-politicalboundary", "theme1"}

I would like to import just one layer ("ke-landcover" for example) from the data but I can't figure this out.

If I ask for the "Graphics" element, I think I get the layers individually, but I don't know how to verify whether this is true or not. I'm new to mathematica, and documentation for GIS applications is limited at best. Can anyone help me?

Thanks!

-
Import["file.shp", "Elements"] should tell you which elements are available for import. What does it say in your case? –  cormullion Feb 25 '13 at 15:07
{"Centering", "CentralScaleFactor", "CoordinateSystem", \ "CoordinateSystemInformation", "Data", "Datum", "Graphics", \ "GraphicsList", "GridOrigin", "InverseFlattening", "LayerNames", \ "LayerTypes", "LinearUnits", "Projection", "ProjectionName", \ "ReferenceModel", "SemimajorAxis", "SemiminorAxis", "SpatialRange", \ "StandardParallels"} –  ShaunLangley Feb 25 '13 at 15:15
Extracting "LayerNames" returns {"Africover_Kenya_Spatial_Aggregation", "ke-landcover", "ke-forests", \ "ke-rangeland", "ke-agriculture", "ke-waterbodies", "ke-bareareas", \ "ke-urban", "ke-politicalboundary", "theme1"} and I want just the second element of the data that is contained in this second layer. –  ShaunLangley Feb 25 '13 at 15:16
If I ask for the "Graphics" element, I think each layer is returned separately, but I'm not certain this is the case. I do not know how to ascertain whether or not this is true. –  ShaunLangley Feb 25 '13 at 15:18
Shapefiles do not have layers. Perhaps you are importing an entire directory or archive of shapefiles (a "SHP bundle")? In this case the layer names correspond to the individual shapefiles, which only have to be specifically named to be separately imported. –  whuber Feb 25 '13 at 16:12

Import the data like:

dat = First@Import["/tmp/Archive.zip", "Data"];


A good start after importing a .shp file is to find out what data it contains:

dat /. Rule[a_, b_] :> a
(* {"LayerName", "Geometry", "Labels", "LabeledData"} *)
"Labels" /. dat
(* {"USERLABEL", "CODE1", "CODE2", "CODE3", "LCCCODE", "LCCLEVEL", "HECTARES", "LC"} *)


Let's group the geometry data by the LC label and see what we get:

labels = "LC" /. ("LabeledData" /. dat);
geometry = "Geometry" /. dat;
labeled = {labels, geometry}\[Transpose];
gathered = GatherBy[labeled, First@# &];

gathered//Length (* 6 *)
Union@labels (* {"AG", "BA", "NVT", "NVW", "UR", "WAT"} *)

g = Graphics[MapIndexed[
{ColorData[1] @@ #2, Tooltip[#[[All, 2]], #[[1, 1]]]} &, gathered]
]


Using this and some guessing the labels probably correspond to:

"AG"  Cultivated Terrestrial & Aquatic
"BA"  Bare Areas
"NVT" Natural Vegetation - Terrestrial
"NVW" Natural Vegetation - Aquatic
"UR"
"WAT" Artificial & Natural Water


SHP files have a tendency towards completely cryptic labels

-
Nice spelunking –  belisarius Feb 25 '13 at 22:28
Awesome! Thanks! –  ShaunLangley Feb 25 '13 at 22:37
@ssch where did you learn about how to tease apart the data like that? Your solution was elegant, but I don't understand half of what you did. Do you have any recommendations for resources? Specifically I want to know about how to identify the different elements in an object and then how to interact with them individually like you did above. –  ShaunLangley Feb 26 '13 at 3:47
@ShaunLangley I don't know any good resources :( LabeledData is of the form {"label1"->{"label1 for geometry[[1]]", "label1 for geometry[[2]]", ...}, "label2"->{"label2 for geometry[[1]]", "label2 for geometry[[2]]", ...}, ...} So for instance the first element in geometry has the labels ("LabeledData" /. dat)[[All, 2, 1]] –  ssch Feb 26 '13 at 12:08

At least here this seems to get the whole file, select the second element of the second layer and draw it:

Quiet@Graphics["Geometry" /. Import["c:\usamap.zip", "Data"][[2, 2]]]


Please note that the file is a zip bundle

It has only one layer in it:

Import["c:\\Archive.zip", "LayerNames"]
(* "Africover_Kenya_Spatial_Aggregation" *)


It contains a map:

Import["c:\\Archive.zip"]


-
In my case, when I try and subset the element [[2,2]], it returns > Frame -> True –  ShaunLangley Feb 25 '13 at 16:20
@ShaunLangley So your comment "I want just the second element of the data that is contained in this second layer" seems misleading –  belisarius Feb 25 '13 at 16:23
I'm not sure how to phrase the question then. The data has multiple overlapping layers of data, identified with "LayerNames". I need to extract the data contained only in the second layer, but I don't know how to isolate that layer. –  ShaunLangley Feb 25 '13 at 16:29
@ShaunLangley Is Import["c:\usamap.zip", "Data"][[2]] what you want? –  belisarius Feb 25 '13 at 16:32
@ShaunLangley I guess it'll be easier if you upload your file somewhere –  belisarius Feb 25 '13 at 16:56