I can get the postion with this code

pos = First[
  Flatten[#, 1] & /@ 
      "AdministrativeDivision", {"California", "UnitedStates"}][
     EntityProperty["AdministrativeDivision", "Polygon"]]]]

You can see the postion in the map


Mathematica graphics

I can get their elevation

eleData = 
   Flatten[#, 1] & /@ 
       "AdministrativeDivision", {"California", "UnitedStates"}][
      EntityProperty["AdministrativeDivision", "Polygon"]]]]]

Then I get the data

data = Flatten /@ Transpose[{pos, List /@ eleData}]

I can plot its discrete plot


Mathematica graphics

But how to connected those discrete points to get a smooth boundary?

  • $\begingroup$ If the points are returned in order, you could try ListLinePlot[pos]. The islands might cause a few problems. $\endgroup$ – QuantumDot Sep 17 '17 at 12:45
  • $\begingroup$ @QuantumDot ListLinePlot can plot 3D graphics? $\endgroup$ – yode Sep 17 '17 at 13:24
  • $\begingroup$ Some answers do essentially the same thing here, but with countries instead of a single state. Somewhat related: (60427) $\endgroup$ – Michael E2 Sep 17 '17 at 13:24
  • $\begingroup$ @MichaelE2 Thanks for the links,but in my case,I just want to get the boundary $\endgroup$ – yode Sep 17 '17 at 13:43
  • $\begingroup$ Yep, and that's what the answer in the question does. ("Wireframe": just change to FaceForm[None] or change Polygon to Line and append the first point(s) to the end(s).) $\endgroup$ – Michael E2 Sep 17 '17 at 13:52

Another approach is similar to Yode's answer, but to refine the boundary into smaller segments before calling GeoElevationData.

poly = EntityValue[
  Entity["AdministrativeDivision", {"California", "UnitedStates"}], 
] /. GeoPosition -> Identity;

Refine the boundary by imposing a maximum length:

bd = DiscretizeRegion[RegionBoundary[poly], MaxCellMeasure -> {1 -> .1}];

Now replace each 2D coordinate with it's 3D version. Here we make one bulk call to GeoElevationData to avoid the overhead of many server calls:

raw = GeoElevationData[GeoPosition[MeshCoordinates[bd]], Automatic, "GeoPosition"];
pts3D = First[raw][[All, {2, 1, 3}]];

Now we construct a MeshRegion (or equivalently we could use Graphics3D + GraphicsComplex):

ratio = Divide @@ Subtract @@@ RegionBounds[bd];
MeshRegion[pts3D, MeshCells[bd, 1], BoxRatios -> {ratio, 1, .1}]

enter image description here


Here's another (slightly hacky) approach. The idea is to take the BoundaryMeshRegion from GeoElevationData, take only the top, then take the topological boundary.

First the BoundaryMeshRegion:

cali = Entity["AdministrativeDivision", {"California", "UnitedStates"}];
reg = GeoElevationData[cali, Automatic, "Region"]

enter image description here

To get the top portion of this region, I select the primitives free of a coordinate that's at the bottom of the region:

floor = RegionBounds[reg][[3, 1]];
prims = Select[MeshPrimitives[reg, 2], FreeQ[floor]];

Visualize this:


enter image description here

I then tally the each line segment of each triangle and select the ones with multiplicity one. This will give me the boundary:

segs = Catenate[Partition[#, 2, 1, 1] & /@ prims[[All, 1]]];
boundary = Cases[Tally[Sort /@ segs], {_, 1}][[All, 1]];

To visualize, I exaggerate the z box ratio, assuring the xy ratios are to scale:

ratio = Divide @@ Most[Subtract @@@ RegionBounds[reg]];
Graphics3D[Line[segs], BoxRatios -> {1, ratio, .1}, Boxed -> False]

enter image description here

  • $\begingroup$ Ah, I was trying to separate out the islands with ConnectedMeshComponents[], but my machine hung. :o Nice solution! $\endgroup$ – J. M. will be back soon Sep 17 '17 at 18:43
  • $\begingroup$ It is a magic.... $\endgroup$ – yode Sep 17 '17 at 18:55

Here's one possible method:

california = Entity["AdministrativeDivision", {"California", "UnitedStates"}];

bounds = GeoBounds[california];
cheights = Reverse[QuantityMagnitude[GeoElevationData[Transpose[bounds],
                                                      GeoZoomLevel -> 4,
                                                      UnitSystem -> "Metric"]]];
crf = RegionMember[MapAt[Map[Reverse[#, 2] &, QuantityMagnitude[LatitudeLongitude[#]]] &, 
                         EntityValue[california, "Polygon"], 1]];

ListPlot3D[cheights, BoundaryStyle -> Thick, DataRange -> Reverse[bounds],
           Mesh -> None, PlotStyle -> None, RegionFunction -> (crf[{#1, #2}] &)]

3D boundary of California

Use Cases[] as usual if you need the actual Line[] objects.

  • $\begingroup$ I'm sorry I missed the elevation data information in my original post..That so sorry..But I have updated it. $\endgroup$ – yode Sep 17 '17 at 14:02
  • $\begingroup$ That's OK; I got what you meant, as you can see from the figure in this answer. $\endgroup$ – J. M. will be back soon Sep 17 '17 at 14:03
  • $\begingroup$ Yes,it is indeed. $\endgroup$ – yode Sep 17 '17 at 14:10
  • $\begingroup$ Could we have any fast solution?I like your result but I have to say to get that cheights will cost too many time.. $\endgroup$ – yode Sep 17 '17 at 14:40
  • 1
    $\begingroup$ EntityValue[ Entity["AdministrativeDivision", {"California", "UnitedStates"}], EntityProperty["AdministrativeDivision", "Polygon", {"ZoomLevel" -> 8}]] get better result.And more fast. $\endgroup$ – yode Sep 17 '17 at 16:39

As you geo-polygon has the 2D coordinates you need. However, the coordinate system ofGraphics and GeoGraphics are not the same. You need to Reverse the geo-polygon coordinate pairs to make them compatible in graphics. The GeoPosition head must also be removed. Finally a z-axis coordinate should be added (I used MapAt) for the polygon in Graphics3D.


geoPoly = 
 Entity["AdministrativeDivision", {"California", "UnitedStates"}][
  EntityProperty["AdministrativeDivision", "Polygon"]]


 MapAt[Append[0], {All, All, All}]@
  Reverse[geoPoly /. GeoPosition -> Identity, 4]

Mathematica graphics

Hope this helps.

  • $\begingroup$ I'm sorry I missed the elevation data information in my original post..That so sorry..But I have updated it. $\endgroup$ – yode Sep 17 '17 at 14:01
  • $\begingroup$ EntityValue[ Entity["AdministrativeDivision", {"California", "UnitedStates"}], EntityProperty["AdministrativeDivision", "Polygon", {"ZoomLevel" -> 9}]] get better result $\endgroup$ – yode Sep 17 '17 at 16:41

Considering the J.M and Edmund's answer,I figure out more faster method based on Entity

poly = EntityValue[
    Entity["AdministrativeDivision", {"California", "UnitedStates"}], 
     "Polygon", {"ZoomLevel" -> 6}]] /. GeoPosition -> Identity;
shape = MapAt[
   Append[Reverse[#], QuantityMagnitude[GeoElevationData[#]]] &, 
   poly, {All, All, All}];
Graphics3D[{FaceForm[], shape}, BoxRatios -> {1, 1, 1/3}]

enter image description here

  • 1
    $\begingroup$ So you don't want the real spatial configuration, but one where sea-level is assumed to be flat and lines of longitude are assumed to be parallel? (The non-flat coastline in the graphics is, I assume, due to numerical error in the data. That's not what I'm asking about. The coast viewed from the side should show the earth's curvature and not be roughly level.) $\endgroup$ – Michael E2 Sep 17 '17 at 17:14

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.