Possibly some of the most obsfucated Python code I have ever written
Tags: dissertationMarsPythonMars Express HRSC03 Dec 2016 - MawKernewek
I have recently been reviewing the Python code from my MSc dissertation which read in the Mars Express High Resolution Stereo Camera digital terrain model tiles, and produced histograms and statistics.
The code was not very well written, as much was edited hastily towards the end of the data analysis phase of my dissertation. However I have managed to do something even more silly since then.
The command gdalinfo allows some summary statistics of a GDAL compatible file to be output to the command line:
gdalinfo h0037_0000_da4.img
Driver: PDS/NASA Planetary Data System
Files: h0037_0000_da4.img
Size is 4391, 36200
Coordinate System is:
PROJCS["SINUSOIDAL MARS",
GEOGCS["GCS_MARS",
DATUM["D_MARS",
SPHEROID["MARS",3396000,0]],
PRIMEM["Reference_Meridian",0],
UNIT["degree",0.0174532925199433]],
PROJECTION["Sinusoidal"],
PARAMETER["longitude_of_center",227],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0]]
Origin = (-55150.000000000000000,2645890.000000000000000)
Pixel Size = (50.000000000000000,-50.000000000000000)
Metadata:
BANDWIDTH=-1e+32
CENTER_FILTER_WAVELENGTH=-1e+32
DATA_SET_ID="MEX-M-HRSC-5-REFDR-DTM-V1.0"
FILTER_NAME=
INSTRUMENT_ID=HRSC
INSTRUMENT_NAME="HIGH RESOLUTION STEREO CAMERA"
MISSION_NAME="MARS EXPRESS"
NOTE=
PRODUCER_INSTITUTION_NAME=
PRODUCT_CREATION_TIME=2008-06-11T06:53:45.000Z
PRODUCT_ID="H0037_0000_DA4.IMG"
PRODUCT_TYPE=
SPACECRAFT_NAME=
TARGET_NAME=MARS
Corner Coordinates:
Upper Left ( -55150.000, 2645890.000) (134d18'27.69"W, 44d38'24.94"N)
Lower Left ( -55150.000, 835890.000) (133d57'33.77"W, 14d 6' 9.93"N)
Upper Right ( 164400.000, 2645890.000) (129d 6' 6.55"W, 44d38'24.94"N)
Lower Right ( 164400.000, 835890.000) (130d 8'24.44"W, 14d 6' 9.93"N)
Center ( 54625.000, 1740890.000) (131d56'32.83"W, 29d22'17.44"N)
Band 1 Block=4391x1 Type=Int16, ColorInterp=Undefined
Min=-2721.000 Max=21234.000
Minimum=-2721.000, Maximum=21234.000, Mean=2131.470, StdDev=6496.150
NoData Value=-32768
Metadata:
STATISTICS_MAXIMUM=21234
STATISTICS_MEAN=2131.47
STATISTICS_MINIMUM=-2721
STATISTICS_STDDEV=6496.15
The original purpose of using this in my Python was to find out which hemisphere we were in, for the processing of the aspect data, which was to separate the northern and southern hemisphere tiles.
So therefore I wrote this:
findHemisph_cmd = "gdalinfo *da4*img | grep 'Center' > CentreLoc.txt"
os.system(findHemisph_cmd)
HemiSph = CentreLoc.split('"')[-1][0] # N or S hopefully
Recently, I though it useful to actually get the locations of the centres of each tile, so that output histograms for individual tiles can be labelled.
So I took the "Centre" line from the output of gdalinfo, and converted the degrees and minutes to decimal degrees accurate to one decimal place:
centrelocation = CentreLoc.split("(")[-1].strip()
centrelocation = centrelocation[:-1]
centrelocation = DMSstr2dec(centrelocation)
The next thing I wanted to do, was to make the tiles appear in a specific order, namely northern hemisphere first, going from -180° to 180° longitude, then the same for the southern hemisphere:
The filename of the histogram above is prefixed by frame0048000606h0037_0000 which since it is at 132°W, this is converted to 48°, which produces the first part of the string as "00480". The latitude is 29.4°N, which is 60.6° counting from the north pole down which converts to "00606".
I recently discovered that gdalinfo has a -json switch so that it produces output in JSON format, which produces this:
{
"description":"h0037_0000_da4.img",
"driverShortName":"PDS",
"driverLongName":"NASA Planetary Data System",
"files":[
"h0037_0000_da4.img"
],
"size":[
4391,
36200
],
"coordinateSystem":{
"wkt":"PROJCS[\"SINUSOIDAL MARS\",\n GEOGCS[\"GCS_MARS\",\n DATUM[\"D_MARS\",\n SPHEROID[\"MARS\",3396000,0]],\n PRIMEM[\"Reference_Meridian\",0],\n UNIT[\"degree\",0.0174532925199433]],\n PROJECTION[\"Sinusoidal\"],\n PARAMETER[\"longitude_of_center\",227],\n PARAMETER[\"false_easting\",0],\n PARAMETER[\"false_northing\",0]]"
},
"geoTransform":[
-55150.0,
50.0,
0.0,
2645890.0,
0.0,
-50.0
],
"metadata":{
"":{
"BANDWIDTH":"-1e+32",
"CENTER_FILTER_WAVELENGTH":"-1e+32",
"DATA_SET_ID":"\"MEX-M-HRSC-5-REFDR-DTM-V1.0\"",
"FILTER_NAME":"",
"INSTRUMENT_ID":"HRSC",
"INSTRUMENT_NAME":"\"HIGH RESOLUTION STEREO CAMERA\"",
"MISSION_NAME":"\"MARS EXPRESS\"",
"NOTE":"",
"PRODUCER_INSTITUTION_NAME":"",
"PRODUCT_CREATION_TIME":"2008-06-11T06:53:45.000Z",
"PRODUCT_ID":"\"H0037_0000_DA4.IMG\"",
"PRODUCT_TYPE":"",
"SPACECRAFT_NAME":"",
"TARGET_NAME":"MARS"
}
},
"cornerCoordinates":{
"upperLeft":[
-55150.0,
2645890.0
],
"lowerLeft":[
-55150.0,
835890.0
],
"upperRight":[
164400.0,
2645890.0
],
"lowerRight":[
164400.0,
835890.0
],
"center":[
54625.0,
1740890.0
]
},
"wgs84Extent":{
"type":"Polygon",
"coordinates":[
[
[
-134.3076928,
44.6402621
],
[
-133.9593814,
14.1027589
],
[
-129.1018186,
44.6402621
],
[
-130.1401215,
14.1027589
],
[
-134.3076928,
44.6402621
]
]
]
},
"bands":[
{
"band":1,
"block":[
4391,
1
],
"type":"Int16",
"colorInterpretation":"Undefined",
"min":-2721.0,
"max":21234.0,
"minimum":-2721.0,
"maximum":21234.0,
"mean":2131.47,
"stdDev":6496.15,
"noDataValue":-32768.0,
"metadata":{
"":{
"STATISTICS_MAXIMUM":"21234",
"STATISTICS_MEAN":"2131.47",
"STATISTICS_MINIMUM":"-2721",
"STATISTICS_STDDEV":"6496.15"
}
}
}
]
}
This would of course have been much easier to process, with the standard Python JSON library.
I will write another post on the histograms themselves later.
The code was not very well written, as much was edited hastily towards the end of the data analysis phase of my dissertation. However I have managed to do something even more silly since then.
The command gdalinfo allows some summary statistics of a GDAL compatible file to be output to the command line:
gdalinfo h0037_0000_da4.img
Driver: PDS/NASA Planetary Data System
Files: h0037_0000_da4.img
Size is 4391, 36200
Coordinate System is:
PROJCS["SINUSOIDAL MARS",
GEOGCS["GCS_MARS",
DATUM["D_MARS",
SPHEROID["MARS",3396000,0]],
PRIMEM["Reference_Meridian",0],
UNIT["degree",0.0174532925199433]],
PROJECTION["Sinusoidal"],
PARAMETER["longitude_of_center",227],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0]]
Origin = (-55150.000000000000000,2645890.000000000000000)
Pixel Size = (50.000000000000000,-50.000000000000000)
Metadata:
BANDWIDTH=-1e+32
CENTER_FILTER_WAVELENGTH=-1e+32
DATA_SET_ID="MEX-M-HRSC-5-REFDR-DTM-V1.0"
FILTER_NAME=
INSTRUMENT_ID=HRSC
INSTRUMENT_NAME="HIGH RESOLUTION STEREO CAMERA"
MISSION_NAME="MARS EXPRESS"
NOTE=
PRODUCER_INSTITUTION_NAME=
PRODUCT_CREATION_TIME=2008-06-11T06:53:45.000Z
PRODUCT_ID="H0037_0000_DA4.IMG"
PRODUCT_TYPE=
SPACECRAFT_NAME=
TARGET_NAME=MARS
Corner Coordinates:
Upper Left ( -55150.000, 2645890.000) (134d18'27.69"W, 44d38'24.94"N)
Lower Left ( -55150.000, 835890.000) (133d57'33.77"W, 14d 6' 9.93"N)
Upper Right ( 164400.000, 2645890.000) (129d 6' 6.55"W, 44d38'24.94"N)
Lower Right ( 164400.000, 835890.000) (130d 8'24.44"W, 14d 6' 9.93"N)
Center ( 54625.000, 1740890.000) (131d56'32.83"W, 29d22'17.44"N)
Band 1 Block=4391x1 Type=Int16, ColorInterp=Undefined
Min=-2721.000 Max=21234.000
Minimum=-2721.000, Maximum=21234.000, Mean=2131.470, StdDev=6496.150
NoData Value=-32768
Metadata:
STATISTICS_MAXIMUM=21234
STATISTICS_MEAN=2131.47
STATISTICS_MINIMUM=-2721
STATISTICS_STDDEV=6496.15
The original purpose of using this in my Python was to find out which hemisphere we were in, for the processing of the aspect data, which was to separate the northern and southern hemisphere tiles.
So therefore I wrote this:
findHemisph_cmd = "gdalinfo *da4*img | grep 'Center' > CentreLoc.txt"
os.system(findHemisph_cmd)
HemiSph = CentreLoc.split('"')[-1][0] # N or S hopefully
Recently, I though it useful to actually get the locations of the centres of each tile, so that output histograms for individual tiles can be labelled.
So I took the "Centre" line from the output of gdalinfo, and converted the degrees and minutes to decimal degrees accurate to one decimal place:
centrelocation = CentreLoc.split("(")[-1].strip()
centrelocation = centrelocation[:-1]
centrelocation = DMSstr2dec(centrelocation)
 |
| An example of the histograms, showing total area within elevation bins of width 200m, and amount of that area that has slope exceeding 5° and 10°. |
The next thing I wanted to do, was to make the tiles appear in a specific order, namely northern hemisphere first, going from -180° to 180° longitude, then the same for the southern hemisphere:
The filename of the histogram above is prefixed by frame0048000606h0037_0000 which since it is at 132°W, this is converted to 48°, which produces the first part of the string as "00480". The latitude is 29.4°N, which is 60.6° counting from the north pole down which converts to "00606".
I recently discovered that gdalinfo has a -json switch so that it produces output in JSON format, which produces this:
{
"description":"h0037_0000_da4.img",
"driverShortName":"PDS",
"driverLongName":"NASA Planetary Data System",
"files":[
"h0037_0000_da4.img"
],
"size":[
4391,
36200
],
"coordinateSystem":{
"wkt":"PROJCS[\"SINUSOIDAL MARS\",\n GEOGCS[\"GCS_MARS\",\n DATUM[\"D_MARS\",\n SPHEROID[\"MARS\",3396000,0]],\n PRIMEM[\"Reference_Meridian\",0],\n UNIT[\"degree\",0.0174532925199433]],\n PROJECTION[\"Sinusoidal\"],\n PARAMETER[\"longitude_of_center\",227],\n PARAMETER[\"false_easting\",0],\n PARAMETER[\"false_northing\",0]]"
},
"geoTransform":[
-55150.0,
50.0,
0.0,
2645890.0,
0.0,
-50.0
],
"metadata":{
"":{
"BANDWIDTH":"-1e+32",
"CENTER_FILTER_WAVELENGTH":"-1e+32",
"DATA_SET_ID":"\"MEX-M-HRSC-5-REFDR-DTM-V1.0\"",
"FILTER_NAME":"",
"INSTRUMENT_ID":"HRSC",
"INSTRUMENT_NAME":"\"HIGH RESOLUTION STEREO CAMERA\"",
"MISSION_NAME":"\"MARS EXPRESS\"",
"NOTE":"",
"PRODUCER_INSTITUTION_NAME":"",
"PRODUCT_CREATION_TIME":"2008-06-11T06:53:45.000Z",
"PRODUCT_ID":"\"H0037_0000_DA4.IMG\"",
"PRODUCT_TYPE":"",
"SPACECRAFT_NAME":"",
"TARGET_NAME":"MARS"
}
},
"cornerCoordinates":{
"upperLeft":[
-55150.0,
2645890.0
],
"lowerLeft":[
-55150.0,
835890.0
],
"upperRight":[
164400.0,
2645890.0
],
"lowerRight":[
164400.0,
835890.0
],
"center":[
54625.0,
1740890.0
]
},
"wgs84Extent":{
"type":"Polygon",
"coordinates":[
[
[
-134.3076928,
44.6402621
],
[
-133.9593814,
14.1027589
],
[
-129.1018186,
44.6402621
],
[
-130.1401215,
14.1027589
],
[
-134.3076928,
44.6402621
]
]
]
},
"bands":[
{
"band":1,
"block":[
4391,
1
],
"type":"Int16",
"colorInterpretation":"Undefined",
"min":-2721.0,
"max":21234.0,
"minimum":-2721.0,
"maximum":21234.0,
"mean":2131.47,
"stdDev":6496.15,
"noDataValue":-32768.0,
"metadata":{
"":{
"STATISTICS_MAXIMUM":"21234",
"STATISTICS_MEAN":"2131.47",
"STATISTICS_MINIMUM":"-2721",
"STATISTICS_STDDEV":"6496.15"
}
}
}
]
}
This would of course have been much easier to process, with the standard Python JSON library.
I will write another post on the histograms themselves later.