PDS tabular data to shapefile and images (MESSENGER Magnetic Crustal Field Maps)

Quite often PDS “map” data is released in a simple Latitude/Longitude table format. Conversion to another form (say for a GIS) can be trickier than most would think. It can also be hard to visualize some of this data if not rasterized (converted to an image). Fortunately, there are actually lots of free methods to help with conversions.

When I get asked about table data, I usually ask for more context:

(1) Is the data perfectly regular? Meaning are the Lat/Lons equally spaced (e.g., every half-degree). If so, then tools like GMT’s xyz2grd or GDAL’s “XYZ” format or specialty routine gdal_grid can help. There are R routines and many methods in Python to handle this also.

(2) If the data is more random (irregular lat/lon spacing), then one usually needs to set an interpolation to fill (or not fill) holes. Then gdal_grid can also help but also GMT has many interpolation methods (like “surface”). Also R (also free) has “rasterFromXYZ” and Python users can use numpy or scipy just name a couple.

For a Python example, see this post from kidpixo for similar tabular data. How is the below method different? I am using GDAL’s PDS4 table support and the new planetary code capabilities in PROJ library (which GDAL uses) with the goal to write out a GIS shapefile and several GeoTIFF images.

Of course, other commercial apps like ArcGIS, ENVI, Global Mapper can do these conversions too.

Data: MESSENGER Magnetic Crustal Field Maps (by Hood, L.L.) can be found at the PPI node: https://pds-ppi.igpp.ucla.edu/search/view/?f=yes&id=pds://PPI/mess-mag-crustal-field-map

This data is Lat/Lon table which is regularly gridded using a 0.5 (Lat) x 1.0 (Lon) degree spacing. The specific file is Bm_output_0E360E which is 0E to 360E Longitude and 35N to 75N in Latitude. To use, download both the *.tab (table) and the *.xml (PDS4 label). You might quickly notice this data in 2D-map form (Cartesian) is rectangular (0.5x1 degree). It does make sense for this data since it is pole-ward, thus in a spherical system it is closer to square (perfectly square at Lat=60N). But storing files in a spherical system is not common for 2D mapping applications. So for this example, we will keep the original pixels as rectangular. This (non-square pixels) is not well supported in many planetary formats but is supported in PDS4 (and has been in GeoTIFF since the beginning).

So far I have only tested conversion using GDAL v3.4 routines (not GMT or R or Python routines). Fortunately, GDAL had no issue supporting this PDS4 tabular data which made the decision to use GDAL routines easy. For a GDAL env., I currently recommend using Anaconda (any OS).

Conversion to GIS point Shapefile (or many other export options available in ogr2ogr): using GDAL version 3.4 – first test read support for the PDS4 table ogrinfo Bm_output_0E360E.xml

Now convert to GIS Shapefile, where -a_srs (input spatial reference) is set to using newly added IAU_2015 code for Mercury (in degrees). The GIS shapefile, can be used directly in most GIS applications for visualization but many will still refer to visualize an image (see below).

ogr2ogr -f "ESRI Shapefile" -a_srs IAU_2015:19900 BM_output_0E360E_points.shp Bm_output_0E360E.xml

Or go straight to an image: Conversion to a GeoTIFF image (raster)

#  using GDAL version 3.4 -- convert to a PDS4 (GeoTIFF) image
#  where -a_srs is set to new IAU_2015 code for Mercury, in degrees), although for images, GDAL will
#       convert to an Equirectangular projection in meters for the PDS4 label.
#  -tr = target resolution (in degrees) for the GeoTIFF, PDS4 pixel resolution will be converted to meters.
#       note the input spacing results in a rectangular output pixel. This is supported in PDS4 but not
#       common (or supported) in other planetary formats (e.g., PDS3, VICAR, ISIS)
#  -txe and -tye are the target X/Y extents in degrees. This helps define a hard boundary during the 
#       conversion from points to pixels (using a "-a" nearest neighbor algorithm for the conversion)

gdal_grid -zfield "Field Magnitude" -a nearest:radius1=0.0:radius2=0.0:angle=0.0:nodata=0.0 -ot Float32 -of GTIFF -tr 1.0 0.5 -txe 0 360 -tye 75 35 -a_srs IAU_2015:19900 Bm_output_0E360E.xml Bm_output_0E360E_mag.tif
gdal_grid -zfield "Field Radial"    -a nearest:radius1=0.0:radius2=0.0:angle=0.0:nodata=0.0 -ot Float32 -of GTIFF -tr 1.0 0.5 -txe 0 360 -tye 75 35 -a_srs IAU_2015:19900 Bm_output_0E360E.xml Bm_output_0E360E_rad.tif
gdal_grid -zfield "Field East"      -a nearest:radius1=0.0:radius2=0.0:angle=0.0:nodata=0.0 -ot Float32 -of GTIFF -tr 1.0 0.5 -txe 0 360 -tye 75 35 -a_srs IAU_2015:19900 Bm_output_0E360E.xml Bm_output_0E360E_east.tif
gdal_grid -zfield "Field North"     -a nearest:radius1=0.0:radius2=0.0:angle=0.0:nodata=0.0 -ot Float32 -of GTIFF -tr 1.0 0.5 -txe 0 360 -tye 75 35 -a_srs IAU_2015:19900 Bm_output_0E360E.xml Bm_output_0E360E_north.tif
#  note: alt is al the same value "40", not converted
#gdal_grid -zfield "Altitude"        -a nearest:radius1=0.0:radius2=0.0:angle=0.0:nodata=0.0 -ot Float32 -of GTIFF -tr 1.0 0.5 -txe 0 360 -tye 75 35 -a_srs IAU_2015:19900 Bm_output_0E360E.xml Bm_output_0E360E_alt.tif


#  (optional) Now add PDS4 label to existing GeoTIFF using the original PDS table label as the input template
#     which helps set metadata no longer supported in the created output GeoTIFF. Note the input
#     PDS4 template was update to 1.G.0.0 (1.16.0.0) to support more recent updates in CART (optional).
#  Generally the output PDS4 label will need some minor touch-ups (but it should validate once the latest GDAL v3.4 is released )
#  -of = output format
#  -co = conversion option

gdal_translate -of PDS4 -co CREATE_LABEL_ONLY=yes -co TEMPLATE=Bm_output_0E360E.xml Bm_output_0E360E_mag.tif   Bm_output_0E360E_mag.xml
gdal_translate -of PDS4 -co CREATE_LABEL_ONLY=yes -co TEMPLATE=Bm_output_0E360E.xml Bm_output_0E360E_rad.tif   Bm_output_0E360E_rad.xml
gdal_translate -of PDS4 -co CREATE_LABEL_ONLY=yes -co TEMPLATE=Bm_output_0E360E.xml Bm_output_0E360E_east.tif  Bm_output_0E360E_east.xml
gdal_translate -of PDS4 -co CREATE_LABEL_ONLY=yes -co TEMPLATE=Bm_output_0E360E.xml Bm_output_0E360E_north.tif Bm_output_0E360E_north.xml

# (optional) convert GeoTIFFs to simple stretched 8bit jpegs (optional browse export)
gdal_translate -of JPEG -scale -ot Byte Bm_output_0E360E_mag.tif   Bm_output_0E360E_mag_browse.jpg
gdal_translate -of JPEG -scale -ot Byte Bm_output_0E360E_rad.tif   Bm_output_0E360E_rad_browse.jpg
gdal_translate -of JPEG -scale -ot Byte Bm_output_0E360E_east.tif  Bm_output_0E360E_east_browse.jpg
gdal_translate -of JPEG -scale -ot Byte Bm_output_0E360E_north.tif Bm_output_0E360E_north_browse.jpg