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

Posts Tagged ‘gis’

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:

5586257.05,-497123.81,1.5
5586261.30,-497072.70,1.2
5586410.36,-486191.97,0.4
5550459.23,-264084.09,-0.3

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

Longitude,Latitude,Elevation
5586257.05,-497123.81,1.5
5586261.30,-497072.70,1.2
5586410.36,-486191.97,0.4
5550459.23,-264084.09,-0.3

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:

<OGRVRTDataSource>
 <OGRVRTLayer name="test_lon_lat_elev_proj">
 <SrcDataSource>test_lon_lat_elev_proj.csv</SrcDataSource>
 <GeometryType>wkbPoint</GeometryType>
 <LayerSRS>NULL</LayerSRS>
 <GeometryField encoding="PointFromColumns" x="Longitude" y="Latitude"/>
 </OGRVRTLayer>
 </OGRVRTDataSource>

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.

 

 

 

 

 

 

Advertisements

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 »