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.
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.