4
$\begingroup$

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!

$\endgroup$
9
  • $\begingroup$ Import["file.shp", "Elements"] should tell you which elements are available for import. What does it say in your case? $\endgroup$
    – cormullion
    Commented Feb 25, 2013 at 15:07
  • $\begingroup$ {"Centering", "CentralScaleFactor", "CoordinateSystem", \ "CoordinateSystemInformation", "Data", "Datum", "Graphics", \ "GraphicsList", "GridOrigin", "InverseFlattening", "LayerNames", \ "LayerTypes", "LinearUnits", "Projection", "ProjectionName", \ "ReferenceModel", "SemimajorAxis", "SemiminorAxis", "SpatialRange", \ "StandardParallels"} $\endgroup$ Commented Feb 25, 2013 at 15:15
  • $\begingroup$ 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. $\endgroup$ Commented Feb 25, 2013 at 15:16
  • $\begingroup$ 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. $\endgroup$ Commented Feb 25, 2013 at 15:18
  • $\begingroup$ 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. $\endgroup$
    – whuber
    Commented Feb 25, 2013 at 16:12

2 Answers 2

6
$\begingroup$

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

$\endgroup$
4
  • $\begingroup$ Nice spelunking $\endgroup$ Commented Feb 25, 2013 at 22:28
  • $\begingroup$ Awesome! Thanks! $\endgroup$ Commented Feb 25, 2013 at 22:37
  • $\begingroup$ @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. $\endgroup$ Commented Feb 26, 2013 at 3:47
  • $\begingroup$ @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]] $\endgroup$
    – ssch
    Commented Feb 26, 2013 at 12:08
3
$\begingroup$

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

$\endgroup$
8
  • $\begingroup$ In my case, when I try and subset the element [[2,2]], it returns > Frame -> True $\endgroup$ Commented Feb 25, 2013 at 16:20
  • $\begingroup$ @ShaunLangley So your comment "I want just the second element of the data that is contained in this second layer" seems misleading $\endgroup$ Commented Feb 25, 2013 at 16:23
  • $\begingroup$ 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. $\endgroup$ Commented Feb 25, 2013 at 16:29
  • $\begingroup$ @ShaunLangley Is Import["c:\usamap.zip", "Data"][[2]] what you want? $\endgroup$ Commented Feb 25, 2013 at 16:32
  • 1
    $\begingroup$ @ShaunLangley I guess it'll be easier if you upload your file somewhere $\endgroup$ Commented Feb 25, 2013 at 16:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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