Howto: Display 2D plot in 3D

In GIS field, sometimes we like to stack several 2D plots together and display them inside a 3D box. For Mathematica, we need to define a function to convert a 2D plot into a 3D graphic object. I will use a small Geotiff as an example, you can download it here if you like to try the code.

(*get the ElevationRange, then import data *)

Import["smalldem.tif", {"Geotiff", "ElevationRange"}]

data = Import["smalldem.tif", {"Geotiff", "Data"}];

Then we a create the contour plot with 100 feet contours

c1=ListContourPlot[data,MaxPlotPoints->30,Contours->Function[Range[650,850,100]],ColorFunction->”DarkTerrain”,PlotRange->{640,850}]

Here the function we need to convert 2D plot into 3D

to3d[plot_,height_,opacity_]:=Module[{newplot}, newplot = First@Graphics[plot];newplot=N@newplot /. {x_?AtomQ,y_?AtomQ}->{x,y, height} ;
newplot /. GraphicsComplex[xx__]->{Opacity[opacity], GraphicsComplex[xx]}];

This function has three parameters: 2D plot, height, and opacity

Let’s create two more contour plots with 50 and 20 feet contours respectively. Then we can stack them together by setting them in different heights.

Show[{Graphics3D[to3d[c1,30,0.75]]}, Graphics3D[to3d[c2,20,0.75]], Graphics3D[to3d[c3,10,0.75]], Lighting->"Neutral", BoxRatios->{1,1,0.8},Axes->True]

Then we like to stack the original geotiff at the very bottom. This time we need to convert the raster into 3D. I use the example you can find in Listplot3D (check the section of “Neat Examples”).

r1=ReliefPlot[data,ColorFunction->colorf,ImagePadding->None, Frame->False, ImageSize->{800,800}];

pic = Reverse[ImageData[r1]];

bg=ListPlot3D[Table[1,{x,1,800,5},{y,1,800,5}],Mesh->None, VertexColors->{pic[[1;;800;;5,1;;800;;5]]},DataRange->{{1,800},{1,800}}, Lighting->"Neutral"]

The final product:

Here is the complete notebook.

Display large DEM with MaxPlotPoints

MaxPlotPoints is a very useful option when works with large DEM.

dem=Import[tif,{"GeoTIFF","Data"}];
Dimensions[dem]
{2825, 5075}

Try ReliefPlot[dem], it takes a very long processing time in my not-up-to-date computer, on the other hand, with limited monitor resolutions, it doesn’t make sense to draw every point of this DEM.

Let’s first try:

ReliefPlot[dem, MaxPlotPoints –> 25]

 

Not quite good, the pixel is in rectangle shape rather than square, then with the small adjustment:

ReliefPlot[dem, MaxPlotPoints -> {50, 25}]

 

ReliefPlot[dem, MaxPlotPoints -> {1000, 500},
ColorFunction -> "GreenBrownTerrain"]

Click to see the large image, please ignore the left part, it is a quality issue from the DEM data source.

The same rule goes with ListPlot3D, if you want to get a 3D version.

ListPlot3D[dem, MaxPlotPoints -> {600, 300}, Mesh -> None,
Axes -> None, ColorFunction -> "GreenBrownTerrain",
BoxRatios -> {2, 1, 0.2}, ImageSize -> 1200]

Import water data from National Weather Service

Besides real-time water data from USGS in the previous post, we like to import more data from National Weather Service, it has the forecast data on water stage.

The data is distributed in XML format.

NWSxml=Import[http://www.crh.noaa.gov/ahps2/xml/nori3_hgirg.xml];

(* function to extract data from xml *)

extractdata[xml_,group_]:=Module[{dataxml,groupdata},dataxml=Cases[NWSxml,XMLElement[group,_,_],Infinity];
groupdata=Transpose[Cases[dataxml,XMLElement["datum",{},{___,XMLElement[#,_,{x_}],___}]:>x,Infinity]&/@{"valid","primary","secondary"}];
groupdata /. {x_,y_,z_}:> {StringReplace[StringDrop[x,-9],{"T"->" "}],ToExpression[y],ToExpression[z]*1000}];

(* extract data*)

observeddata = extractdata[NWSxml, "observed"];
forecastdata = extractdata[NWSxml, "forecast"];

(* produce the graph *)

Flood stage is at 11.0 feet, it seems the white river may be in the flood stage in next several days.

Import USGS national real-time water data

The built-in real-time weather data in Mathematica 7.0 is wonderful. You may be interested in importing more related data.

USGS has web-interface for Daily Streamflow Conditions

Here is the simple steps from importing Gauge height and Discharge data

USGS 03353800 White Lick Creek at Mooresville, IN

(*import html as plain text *)

rawdata=Import[“http://waterdata.usgs.gov/nwis/uv?cb_00060=on&cb_00065=on&format=html&period=7&site_no=03353800”];

(* extract table information *)

rawdata=StringReplace[rawdata,","->""];
result=StringCases[rawdata,x:DatePattern[{"Month","/","Day","/","Year"," ","Hour",":","Minute"}]~~Whitespace~~y:NumberString~~Whitespace~~z:NumberString:>{x,ToExpression[y],ToExpression[z]}];

(*draw graph*)

DateListPlot[result[[All,1;;2]],PlotStyle->{Blue, PointSize[Medium]},PlotLabel->"Site: 03353800",Joined->True,Filling->Bottom,FrameLabel->{"Time","Gage Height (feet)"},DateFunction:>(DateList[{#,{"Month","/","Day","/","Year"," ","Hour",":","Minute"}}]&),DateTicksFormat->{"MonthShort","/","Day"," ","Hour",":","Minute"}]

By the way, the code is tested in Mathematica 6.0.

Create the map with a image background

Here is something you can do with Mathematica 7.0 Image Processing

Here is a good picture of  American Flag

Then create a mask for it and add them together:

mask= Graphics[{Black, CountryData["UnitedStates", {"Polygon", "Mercator"}]}, ImageSize -> {500, 335}]

flagImage = ImageAdd[img, mask]

 

Then create another new image with the boundary:

bgk = Graphics[{FaceForm[], EdgeForm[{Black}], CountryData["UnitedStates", {"Polygon", "Mercator"}]}, ImageSize -> {500, 335}]

ImageMultiply[bgk, flagImage]