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

Archive for the ‘Planetary Science’ Category

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 »