Tell me more ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.
up vote 3 down vote favorite
1

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!

share|improve this question
Import["file.shp", "Elements"] should tell you which elements are available for import. What does it say in your case? – cormullion Feb 25 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 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 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 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 at 16:12
show 4 more comments

2 Answers

up vote 3 down vote accepted

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]
    ]

kenya

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

share|improve this answer
Nice spelunking – belisarius Feb 25 at 22:28
Awesome! Thanks! – ShaunLangley Feb 25 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 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 at 12:08
up vote 3 down vote

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]]]

Mathematica graphics

Please note that the file is a zip bundle

Edit: using your uploaded data file

It has only one layer in it:

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

It contains a map:

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

enter image description here

share|improve this answer
In my case, when I try and subset the element [[2,2]], it returns > Frame -> True – ShaunLangley Feb 25 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 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 at 16:29
@ShaunLangley Is Import["c:\usamap.zip", "Data"][[2]] what you want? – belisarius Feb 25 at 16:32
1  
@ShaunLangley I guess it'll be easier if you upload your file somewhere – belisarius Feb 25 at 16:56
show 3 more comments

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.