A halfway home for wayward thoughts.

  • RSS 12 Angry Men

    • 2012 in review December 31, 2012
      The WordPress.com stats helper monkeys prepared a 2012 annual report for this blog. Here’s an excerpt: 19,000 people fit into the new Barclays Center to see Jay-Z perform. This blog was viewed about 77,000 times in 2012. If it were a concert at the Barclays Center, it would take about 4 sold-out performances for that […]
      Angry Political Optimist
  • Latitude

Planetary mapping with QGIS Part Three – Creating Shapefiles From Planetary Datasets

Posted by Jeff on September 16, 2011

In Part One, I covered how to set up the QGIS workspace to correspond to a custom ellipsoid (Mercury in this case), rather than a terrestrial one.

Part Two dealt with creating polygons and using the GDAL/OGR command line to transform the resulting shapefiles into raw coordinate lists.

This time, I’ll cover the process of creating shapefiles from raw data and importing them into QGIS.

Let’s say I have a file of tab-delimited values of   longitude    latitude    elevation
131.1758   -11.6734   1.5
131.1759   -11.6722   1.2
131.1794   -11.4167   0.4
130.3352   -6.20119   -0.3

To get this into meter map-coordinates, use the PROJ.4 utility “proj” to convert them, basically reversing the process from Part Two.

proj +proj=eqc +lat_ts=0 +a=2440000 +b=2440000 +no_defs test_lon_lat_elev.tab | sed 's/\t/,/g' > test_lon_lat_elev_proj.csv

Where test_lon_lat_elev.tab is the input file listed above and test_lon_lat_elev_proj.csv is the comma-delimited output file. Note that I used the sed utility to convert the tabs to commas.

The resulting file should look like this:


Now we need to add the header row with field descriptions so that the file becomes this:


As far as I can tell, PROJ.4 prefers a tab-delimited input file. If you try to use a comma-delimited one, it fails to properly parse the long/lat/alt values.
OGR, on the other hand, needs comma-delimited files as input, since it uses the concept of “drivers” to parse the file, as will be shown next.

Now that the input is properly formatted, we can use OGR to create the shapefiles.  OGR uses the concept of “virtual formats” to describe the data contained in the input file.  These are basically XML-styles that must be prepared for each file to be converted.  A .csv file will use the following VRT file:

 <OGRVRTLayer name="test_lon_lat_elev_proj">
 <GeometryField encoding="PointFromColumns" x="Longitude" y="Latitude"/>

In the above .vrt file, test_lon_lat_elev_proj.csv is the file produced above.
OGRVRTLayer name must apparently match exactly with the SrcDataSource field.

Now, using the OGR command line:
ogr2ogr test test_lon_lat_elev.vrt
Will produce a directory called test containing, among other things a .shp shapefile.

If you still have the Mercury global mosaic loaded in the QGIS workspace, you can click on the Open Shapefile button, navigate to the test directory that was created above, and then select the new .shp file. The result should be something similar to the image below.

Note that it only appears that there are two points.  There are actually four, which can be seen if you zoom in on the cluster far enough to resolve them.







Posted in Planetary Science | Tagged: , , , , , | Leave a Comment »

Planetary Mapping with QGIS Part 2 – Working With Polygons and Converting Meters to Longitude/Latitude

Posted by Jeff on September 14, 2011

After working through Part One of this series, we should now have a workspace that is properly aligned to Mercury’s surface.

A common task in planetary mapping is identification surface features. One of the most prominent features of Mercury (and of the Moon, and many other planetary bodies) is the abundance of craters on the surface; and in order to do statistical measurement of crater types, densities, sizes, etc. they often need to be converted from images to vector shapes… and this is most often done manually.

QGIS has three basic vector types that can be used… point, line, and polygon. I’m going to skip over points and lines for the moment and talk about polygons. Polygons are stored in the ESRI Shapefile format, which can be read by many different statistics and mapping software packages. As a quick introduction to the process of using polygons in QGIS, I’ll trace the perimeter of Rembrandt Basin.

Using the same image from Part One, zoom in on the area just below the equator and just to the east of the center of the image (the dark vertical band in the image.) If the CRS for the Coordinate Capture plugin is still set correctly, you can use that to find the area of 88E longitude and 33S latitude.

Now click on the New Shapefile button in the top menu bar. The dialog box that appears should be similar to the one shown below:

Change the Type to Polygon, leave the CRS ID alone (since it should have inherited the correct CRS from the project settings), and set the Name field to Rembrandt.
Leave everything else alone and click OK.
You will then be prompted to save the new shapefile. Call it Rembrandt and save it wherever you wish.

This creates a new layer, which is added to the left side layer list.
To add a polygon to this layer, click the Toggle Editing button in the menu bar (). This basically “opens” the shapefile for changes.
Now start adding the polygon by first clicking on the polygon icon in the menu bar () and then left click on the image to place the first vertex.

Each left click will add a vertex, so continue to trace the edges of the crater until you are ready to add the last point. Hover the mouse over that point and this time RIGHT click to add the last vertex and close the polygon.

At this point, a dialog box will appear prompting for an ID for this polygon. Each shapefile can have multiple polygons (or points, or lines) so each one needs to be identified uniquely. For some reason, QGIS doesn’t simply increment the ID number intelligently, so the user is forced to keep typing new (unique) numbers each time.

The area of Rembrandt that was just highlighted is probably now opaque. To solve this, go to Layer->Properties under the Style tab there is a slider which controls transparency.

After clicking OK, click on the Toggle Editing button () again. This will prompt to save the file. Click Save to save and close the new shapefile.

Remember that the new shapefile is stored in its own file, separate from the project file, so you can import that file into other programs, and if the project file is destroyed, the shapefile can still be used.

How to deal with shapefiles programatically and convert meter coordinates to longitude/latitude:

I do a lot of custom programming and scripting, either in MATLAB, Python, BASH, or C, so I often need to be able to grab the raw vertex coordinates out of the shapefile. To do this, I use the tools that QGIS was built upon: GDAL/OGR.
If OGR is installed (and it should be, since QGIS is installed), the following command dumps the shapefile data to the screen:
> ogrinfo -al rembrandt.shp

The result will be the layer name, number of features (in this case, polygons that it contains), projection information, and then the coordinates of the polygon vertices, such as the following:

id (Integer) = 1
POLYGON ((3429285.222557838074863 -1155233.898176640970632,3510078.008147859945893 -1071392.328224731609225,3721968.521299049258232 -1014989.817529810592532,3856115.033222104422748 -1013465.425348866847344,3973493.231154777575284 -1074441.112586619099602,4081725.076001788023859 -1164380.251262303907424,4174712.999039360322058 -1280234.057014033198357,4194530.097391629591584 -1412856.176756144501269,4156420.292868034448475 -1514990.45287937973626,4103066.566535001154989 -1636941.827354884473607,4007029.859135541599244 -1737551.711297175614163,3883554.092479093000293 -1781759.084544546203688,3709773.383851498831064 -1827490.849972860421985,3610687.89209015108645 -1786332.261087377555668,3545139.028309567365795 -1746698.064382838550955,3462821.850538601633161 -1708588.259859243407845,3371358.319681973196566 -1661332.102249985327944,3295138.710634782910347 -1600356.41501223295927,3301236.279358558356762 -1513466.060698435874656,3305809.455901389475912 -1371697.587870661634952,3351541.22132970392704 -1266514.527385538909584,3429285.222557838074863 -1155233.898176640970632))

Notice that the vertex coordinates are in meters rather than a more generally useful longitude/latitude pair.

To get longitude/latitude, we can use the ogr2ogr command to reproject the data. To do that, we’ll recycle the PROJ.4 string from Part One. (Recall that GDAL/OGR is built with PROJ.4 libraries.)

ogr2ogr -t_srs "+proj=longlat +a=2440000 +b=2440000 +no_defs" rembrandt_reproj.shp rembrandt.shp

A new file “rembrandt_reproj.shp” is created. If we then use the same ogrinfo command on this file, the POLYGON section now looks like this:

id (Integer) = 1
POLYGON ((80.526053278338225 -27.127060129494836,82.423219519901195 -25.158302709004971,87.398806455321036 -23.833865898857244,90.548818328104815 -23.798070309393793,93.305078716790618 -25.229893887931873,95.846565568695723 -27.341833666275548,98.0300965259663 -30.062298465497907,98.49543918899117 -33.176514748818235,97.600549452404863 -35.57481924286953,96.347703821184041 -38.438466399945703,94.09258168498657 -40.800975304533537,91.19313893844695 -41.839047398973648,87.112441739613416 -42.912915082877205,84.785728424489022 -41.946434167363996,83.246518077560594 -41.015748841314242,81.313556246534176 -40.120859104727941,79.165820878727047 -39.011195831360929,77.376041405554446 -37.579372252822843,77.519223763408263 -35.539023653406069,77.626610531798605 -32.210033833305026,78.700478215702162 -29.740138160326836,80.526053278338225 -27.127060129494836))

So… voila! We have longitude/latitude pairs! How to get this into a properly-formatted .csv file is left as an exercise for the reader. My solution would involve a few line of BASH script.

Next up: Creating shapefiles from .csv longitude/latitude lists.

Posted in Planetary Science | Tagged: , , , , , | Leave a Comment »

Getting started with QGIS and planetary mapping Part One : Finding Longitude/Latitude on a Map of Mercury

Posted by Jeff on September 13, 2011

This post marks the start to my documentation of using FOSS (Free or Open Source Software) to do planetary mapping. I’ve had to do a lot of background reading, searching, asking, and programming in order to get to a point where I can be productive doing mapping in this software environment, and I’m hoping that these posts will help others if they attempt the same.

First, most (if not all) of the people I know who are involved in planetary mapping use ArcGIS. Having essentially no background in GIS, but nearly twenty years of experience with Linux/UNIX, I opted for an open-source approach.

So, to begin, QGIS can be downloaded here. I’m currently using version 1.7.0 on both Linux and Windows. This version included some updates which I consider to be critical to making the program usable for planetary mapping, such as “on-the-fly reprojection.”

To get started, I’ll use the publicly-available Mercury mosaic which was derived from Mariner 10 and MESSENGER flybys. This can be downloaded in its 450MB glory from the USGS Astrogeology Center. I recommend grabbing the .pgw/.png file, which seems to play more nicely with QGIS than its jgw/jpg sibling. The pgw and jgw files are companion worldfiles which give the GIS software important information about the scaling and dimensions of the geo-referenced image.

To bring up the Mercury mosaic, click on “Add Raster Layer” and browse to the .png file that was loaded in the previous step. The resulting screen should look like this:

Now we need to tell QGIS how the image is projected, since neither the .pgw nor the .png contain this information.

QGIS comes with *many* choices of map projections. Unfortunately, all of them are terrestrial.
To force it to re-project to a different planet (which naturally has a different spherical radius, or alternatively, major and minor ellipsoid axes), we’ll need to add a custom Coordinate Reference System (CRS).
To do this, go to Settings->Custom CRS.
Now click on the Star (to add a new CRS definition) and type in the following (as shown in the image below):

Name: Equidistant Mercury Spherical
Parameters: proj +proj=eqc +lat_ts=0 +a=2440000 +b=2440000 +no_defs

The Name field is not critical. The Parameters field is the one doing all the work. Here, we’re telling it that the projection is in Equidistant Cylindrical (which the USGS calls “Simple Cylindrical”) that the Latitude of True Scale (that is, where the image is not stretched) is at zero degees latitude, and that the major and minor axes are both 2440 km, indicating that it is spherical.

Click the disk icon to make the CRS definition available to the workspace.

So now QGIS knows how to relate the pixels of the image to real-world distances (in meters). But for some odd reason, it doesn’t know how to use this information to convert meters into Longitude and Latitude, so we need to explicitly tell it by creating another custom CRS definition.

QGIS is built using the GDAL/OGR libraries, which handle the translation and projection of geographic file formats, including geo-referenced images. These are, in turn, built upon the PROJ.4 library, which handles the actual reprojection of coordinates.
This means that we can use the internal “+proj … ” strings to force QGIS to use projections which are not built-in, as we have done above.

Click the Star again to create a new definition, and enter the following (as shown in the image below):

Name: Geodetic Mercury Spherical
Parameters: +proj=longlat +a=2440000 +b=2440000 +no_defs

Click on the disk icon to save it, and then click OK.

Almost done for today…

All we’ve done so far is create the new CRS definitions, but we haven’t actually told QGIS to use them yet, so now click on Layer->Set CRS of Layer(s) and scroll down to the bottom. Click on the arrow next to User Defined Coordinate Systems to expand the tree and then click on Equidistant Mercury Spherical and then OK.

But that only set the CRS for that layer. Now set the workspace to that CRS by selecting Layer->Set Project CRS From Layer.

Lastly, in the Coordinate Capture window, click on the globe reticule button to select the CRS to use for coordinate display and like last time, go to the User Defined Coordinate Systems but this time select Geodetic Mercury Spherical and then OK.

You can now click the Start Capture button in the Coordinate Capture window and then click anywhere on the image to get the Longitude, Latitude.

Next Up – Working with Polygons

Posted in Planetary Science | Tagged: , , , , , | 2 Comments »

Vista grr

Posted by Jeff on November 3, 2008

So I got my spiffy (and shiny) new Dell XPS M1130 laptop. (2.1Ghz Core 2 Duo, 4GB mem, 250gig HD, LED backlit 13″ monitor… and Vista.) Each time I started the machine, or resume from hibernation, my hard drive thrashes like crazy for about 20 minutes. You’d think with 4GB of memory, the hard drive would almost never need to be accessed except when loading programs. You’d be wrong. Welcome to SuperFetch. This is a program that runs in the background and pre-loads from disk into memory to allow for faster load times for frequently-used programs. Well, it keeps trying to preload my 8GB virtual machine file (that I use to run Linux concurrently with Vista). Hence the hard drive thrashing. I can sort of see the logic of using a pre-fetch device… on a desktop. A note to Microsoft: I use a laptop, and I like to be able to suspend and re-awaken my laptop in a hurry so that I do things like… I don’t know… move it? When all of that stuff is pre-loaded into memory, before the computer can sleep or suspend, it has to either unload it, or write it out to a temporary disk file. Double the hard disk use, double the memory use, half the convenience.

Thankfully, this demon was pretty easy to exorcise from Vista with the right incantations.

Posted in Computer stuff | Tagged: , , , , | Leave a Comment »

to Yahoo…

Posted by Jeff on October 22, 2008

1) Move Jerry Yang from CEO to CTO.
2) Cut 50% of non-US staff and facilities. Retreat back to your core, original, loyal base.
3) Sell Search. Yes, I realize this was your first service and that all of Yahoo was built from it, but it’s a dog. And you won’t catch Google.
4) Buy Skype from eBay and then wait for eBay to totally tank. Resurrect Yahoo Auctions and link it to Shopping.
5) Forget developing your own peer-networking stuff. It’s been done to death, just use someone else’s API.
6) Flickr is great. Answers is pretty cool. Finance is nice. You have stuff to build on.
7) You are a portal. Make it the best, easiest, nicest place to go for news, stocks, video, etc. iGoogle is nowhere close and MSN may already have you beat. But it’s the only thing that really defines Yahoo as a single entity anymore.

Posted in Thoughts | Tagged: | Leave a Comment »

New volcano simulation

Posted by Jeff on September 22, 2008

My advisor feels strongly that all scientific publications should be peer-reviewed, so in light of Kevin’s valuable insight into my post about the finite element analysis of a point load on an elastic plate, I refined my model and re-ran the simulation. The result is clearly more realistic and can easily be verified by means of a field trip to Hawaii, or by a viewing the award-winning, critically-acclaimed documentary, Joe Versus the Volcano.

New volcano simulation

New volcano simulation

Posted in General | Tagged: | 3 Comments »

Dumb signs

Posted by Jeff on September 18, 2008

Really? As opposed to prepaying after?

Really? As opposed to prepaying after?

You know you haven't read the directions when...

You know you haven't read the directions when...

Posted in General | Tagged: , , | Leave a Comment »

Starting FEA work – flexural plate with a point load

Posted by Jeff on September 5, 2008

Deformation of a buoyant plate with a point load

Deformation of a buoyant plate with a point load

I’m starting to make some headway with the Finite Element Analysis software that we use here (MSC Marc). The image above is the result of a simple model that simulates a large load on a floating edge of a deformable plate. In a geologic sense, this is supposed to represent a volcano sitting on the edge of a crustal plate, which causes the one edge to sink into the buoyant mantle that normally supports it.

Posted in Research | Tagged: , , | 2 Comments »

Resignation Week

Posted by Jeff on August 15, 2008

I suppose there is some philosophical significance to the process of cutting ties here in Ann Arbor to start a new life in Cleveland, and I think that Kushida Sensei would have a decent amount to say about starting any new endeavor with a beginner’s mind… but it’s more than a little overwhelming to close the book on most of the things that I’ve been involved with over the last decade.

On Friday, I worked 16 hours for my last shift in dispatch. That it was a really busy and stressful shift helped to keep me from feeling too badly about walking out of the comm center for the last time. Although I’ll miss most of the people who work in there, I feel pretty relieved that I won’t have to worry about the bureaucracy and policies involved in actually working as a dispatcher.

Yesterday was both my last day in the office with UM (again, I’ll miss the people, not the work environment) and very sadly, my last day as lieutenant of the fire department. While I’ll be staying on with HVA in my capacity as a road medic, the opportunities for continuing in the fire service while in Cleveland are pretty limited. Very few departments are staffed by anything other than career firefighters and I just won’t have the time or flexibility to work a full-time shift schedule.

In any case, it is a huge change in our lives as we start an entirely new chapter and leave the familiar behind.

Posted in Day-to-day | Tagged: , | Leave a Comment »

First ride to campus

Posted by Jeff on August 5, 2008

On monday, I took the bike out for a trial run from the new house to campus. One of our neighbors (Fran) mentioned that it’s a 10-minute ride there and a 15-minute ride back. I managed to get there in about 20, with a 25-minute return trip. Here are some tips for next time to help shave that time down a little:

  • inflate the tires to something more than 20psi
  • Switch the knobby offroad tires to street slicks
  • Write down the directions or take a map
  • Check the brakes… back AND front
  • The bike goes faster when the chain stays on the chainring

For you former Bursley-ites reading this, it turns out that my ride to/from campus involves almost the same elevation changes as the ride from north campus to central campus… approximately 200 feet. However, my ride now has definite bonus of taking me right through Little Italy.

Fortunately, most of the streets I got lost on were pretty nice areas of town, but I did see some areas a few blocks away from our neighborhood that I’m really glad we didn’t buy into (Derbyshire between Lamberton and Cottage Grove). The last time I rode a bike to work or school was about 10 years ago, so it also took me about five minutes to remember how to use a U-lock to property secure a bike to a rack.

Having finally made it to campus, I had a great discussion with Dr. Hauck, my advisor. It looks like I’m going to start my research by learning the finite element analysis software MARC and using it to study the lifecycles of lunar craters. While not quite as exciting as the other option (analysis of IR, visible, and topographic data from the next Mercury Messenger pass) it gets me going on a project quickly.

Posted in Day-to-day | Tagged: , , | Leave a Comment »


Get every new post delivered to your Inbox.