Simple JPEG2000 GML data parser / extractor
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 

640 lines
21 KiB

#!/usr/bin/env python3
# Simple JPEG2000 GML data parser
# Copyright (C) 2019 Pekka Helenius
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.
################################################################
import sys
from warnings import warn as Warn
import os.path
import argparse
import re
import xmltodict
import json
import urllib.request as URL
import math
# TODO import csv
# TODO retrieve place name by using metadata coordinates & ref.system info
# requires internet connection and connection to a valid server
#
# TODO rename JPEG2000 file based on metadata entries, syntax given by user
#
# TODO fix tfw export for JPEG2000 files
################################################################
#
# INPUT ARGUMENTS
argparser = argparse.ArgumentParser()
argparser.add_argument('-i', '--input', help = 'Input JPEG2000 image file', nargs = '?', dest = 'inputfile')
argparser.add_argument('-f', '--dataformat', help = 'Output format (Default: xml; Available: xml | json | [tfw|worldfile] | info)', nargs = '?', dest = 'outputformat')
argparser.add_argument('-o', '--output', help = 'Output file name', nargs = '?', dest = 'outputfile')
argparser.add_argument('-l', '--formatting', help = 'Data formatting (Default: raw; Available: raw, pretty)', nargs = '?', dest = 'formatting')
args = argparser.parse_args()
# Formatting defaults to pretty format
#
if args.formatting is None:
args.formatting = 'pretty'
################################################################
if not len(sys.argv) > 1:
argparser.print_help()
exit(0)
if not args.inputfile.endswith('.jp2'):
Warn("Warning: Not a valid JPEG2000 file suffix")
if args.outputformat is None:
raise ValueError("Error: No output format or file specified")
elif args.outputformat not in ('json', 'xml', 'tfw', 'worldfile', 'info'):
raise ValueError("Error: Not a valid output format")
################################################################
#
# JPEG2000 CHECK STRINGS
# Look for these jp2 strings in file header
#
jp2_header_str = ['ftypjp2', 'jp2 jpx', 'jp2', 'jp2h', 'jp2 J2P1']
# First string to look for in the selected jp2 file
# Basically, we are looking for the start of footer of the file
# which is indicated by this string
#
# We need to convert the string into bytes with encode method
# for the following while loop.
#
mdata_start_str = [ str.encode('gml.data'), str.encode('gxml'), str.encode('fxml') ]
mdata_end_str = str.encode('uuid')
################################################################
#
# OPEN JPEG2000 file
# Open the image file in read-only binary mode
with open(args.inputfile, 'rb') as f:
################################################################
#
# JPEG2000 header check
# Check for the first 4 file lines
header_lines = f.readlines()[0:3]
# Declare a variable to store header text string
header_str = ''
# For each header line 1-4...
for a in header_lines:
# Decode binary formatted string (bytes -> string conversion) and ignore any errors we encounter
header_decode = a.decode('utf-8', errors='ignore')
# Store decoded string into header_str variable
header_str += header_decode
# Check existence of each jp2 specific string we have defined in 'jp2_header_str' list above
for jp2_str in jp2_header_str:
# If a jp2 specific string is found, set up a new variable t and break the for loop
if jp2_str in header_str:
t = ''
break
# Variable t is not defined if any valid jp2 string is not found. Thus, this test gives
# us and exception (NameError) if no any jp2 string is found.
try:
t
except NameError:
raise ValueError("Error: Not a valid JPEG2000 file")
################################################################
#
# PARSE METADATA LINES
# Enumerate all lines, look for metadata start line, using
# string 'mdata_start_str' as a reference
# Break the loop when found. If not found, abort.
#
# Return to the first line again in order to parse footer lines
f.seek(0)
for mstart_num, mstart_line in enumerate(f):
# TODO better formatting for this if statement:
if mdata_start_str[0] in mstart_line or mdata_start_str[1] in mstart_line or mdata_start_str[2] in mstart_line:
break
# else
# TODO echo cannot found metadata start. Abort
# TODO should return value 2
#print(mstart_num)
#sys.exit()
# Enumerate all lines, look for metadata end line, using
# string 'mdata_end_str' as a reference
# Break the loop when found. If not found, abort.
#
f.seek(0)
for mend_num, mend_line in enumerate(f):
if mdata_end_str in mend_line:
break
# else
# TODO echo cannot found metadata end. Abort
# Reset readlines
#
# Convert metadata start line from 'str' type to 'list' with split method
# and merge it with the rest of the metadata lines, defined by readlines method.
# Type of this line list is 'list', thus we use + operator to combine these
# lists.
#
f.seek(0)
metadata_lines = mstart_line.split() + f.readlines()[mstart_num:mend_num]
#mdata_lines = mstart_line.split() + f.readlines()
# Create a new metadata_str variable where we will store our extracted footer strings.
metadata_str = ''
for byteline in metadata_lines:
# Try decode each metadata line to UTF-8 format.
# As these lines are binary code, the conversion will fail for some
# of them. In a case of failure, we let the for loop pass to the next
# line
#
# Add each decoded line into 'footer_str' variable
#
try:
byteline_decoded = byteline.decode('utf-8', errors='strict')
metadata_str += byteline_decoded
except Exception:
pass
f.close()
metadata_xml_all = re.sub(r'(^[^<]*)|([^>]$)', '', metadata_str)
# Create a list element from extracted metadata strings
metadata_xml_all_list = metadata_xml_all.split()
# Find the last element containing <> symbols in metadata_xml_list,
# get the string between them and store it to new variable
for i in reversed(metadata_xml_all_list):
if re.match('</*.*>', i):
last_tag = re.sub('</|>', '', i)
break
# In the original metadata list, find the first occurence of the 'last_tag'
for firstxml_index, value in enumerate(metadata_xml_all_list):
if re.match('<' + last_tag + '>?', value):
break
# For joined metadata list, delete all list entries presented before our 'last_tag'
# Convert list to string format
metadata_parsed_list = metadata_xml_all_list[firstxml_index:]
metadata_joined_list = ' '.join(metadata_parsed_list)
################################################################
#
class GMLDataParser(object):
def __init__(self, datalist):
self.datalist = datalist
def xmlraw(self):
return xmltodict.parse(self.datalist)
def xmlpretty(self):
return xmltodict.unparse(xmltodict.parse(self.datalist),
pretty=True,indent=" ",newl="\n")
def jsonraw(self):
return json.dumps(xmltodict.parse(self.datalist),
separators=(',', ':'))
def jsonpretty(self):
return json.dumps(xmltodict.parse(self.datalist),
indent=2, sort_keys=True)
# Convert GML metadata to JSON tree object
def jsontree(self):
return json.loads(self.jsonpretty())
# Function to get nested key values from JSON data
# by arainchi
# https://stackoverflow.com/a/19871956
def findkey(self, tree, keyvalue):
if isinstance(tree, list):
for i in tree:
for x in self.findkey(i, keyvalue):
yield x
elif isinstance(tree, dict):
if keyvalue in tree:
yield tree[keyvalue]
for j in tree.values():
for x in self.findkey(j, keyvalue):
yield x
gmlparser = GMLDataParser(metadata_joined_list)
gml_json = gmlparser.jsontree()
def findgmlkey(data, gmlkey, num):
try:
return list(gmlparser.findkey(data, gmlkey))[num]
except:
# In a case we can't parse GML data for this element, return string 'Unknown'
return str("Unknown")
################################################################
#
# Extract relevant values for TFW file/Worldfile
class GML_Pos_offsetVectors():
# Sample metadata structure of JPEG2000 files (may differ!):
# offsetVector_1 and offsetVector_2:
#
# gml:FeatureCollection
# gml:featureMember
# gml:FeatureCollection
# gml:featureMember
# gml:RectifiedGridCoverage
# gml:rectifiedGridDomain
# gml:RectifiedGrid
# gml:offsetVector[0]
# #text
# gml:offsetVector[1]
# #text
#
# gml_pos:
#
# gml:FeatureCollection
# gml:featureMember
# gml:FeatureCollection
# gml:featureMember
# gml:RectifiedGridCoverage
# gml:rectifiedGridDomain
# gml:RectifiedGrid
# gml:origin
# gml:Point
# gml:pos
# Find offsetVector elements in the file metadata
# These elements include field #text which we are searching for
gml_offsetVector_1 = findgmlkey(gml_json, '#text', 0)
gml_offsetVector_2 = findgmlkey(gml_json, '#text', 1)
# Check whether we have gml:pos or gml:coordinates element in the file metadata
# gml:coordinates is a deprecated type according to opengis.net
try:
gml_pos = findgmlkey(gml_json, 'gml:pos', 0)
except:
gml_pos = findgmlkey(gml_json, 'gml:coordinates', 0)
# Convert gml_pos to list type in a case it is string type
if type(gml_pos) is str:
# Split values, use any other symbol as a separator except for dot, minus prefix and numbers.
gml_pos = re.split('[^\-^\d^\.]+', gml_pos)
# Get semi-major axis of the Earth from ESPG metadata
# TODO get this actually from metadata!
#try:
# Try to get the value
# Fallback value
#except:
earth_axis_semimajor = 6378137
# Estimated meters for one degree on Earth surface for used ellipsoid model
dec_mult = float((2 * math.pi * earth_axis_semimajor) / 360)
# Declare a new list 'l'
l = []
for d in (gml_offsetVector_1, gml_offsetVector_2):
if type(d) is str:
d = re.split('[^\-^\d^\.]+', d)
# Add extracted value to list 'l'
l += d
# Assumed length of list gml_pos is either 4 (gml:pos) or 6 (gml:coordinates).
# We must treat these list types differently.
# In a case length is either of those, return error.
#
# Map correct gml_pos values into new array 'g'
#
g = [0] * 4
if len(l) == 4:
g[0] = l[0]
g[1] = l[1]
g[2] = l[2]
g[3] = l[3]
elif len(l) == 6:
g[0] = l[3]
g[1] = l[4] # TODO is this correct index?
g[2] = l[2] # TODO is this correct index?
g[3] = l[1]
else:
raise ValueError("Error: Incorrect worldfile metadata definition for rotational and pixel size values")
# World file definition
# https://en.wikipedia.org/wiki/World_file
#
# g[0] = pixel size of X-axis in map units
# g[1] = Y-axis rotation
# g[2] = X-axis rotation
# g[3] = pixel size of Y-axis in map units
# gml_pos[0] = X-coordinate of the center of the upper left pixel
# gml_pos[1] = Y-coordinate of the center of the upper left pixel
# TODO should gml_pos[1] value be decreased by -1?
gml_posinfo = GML_Pos_offsetVectors()
################################################################
#
# ESPG INFORMATION RETRIEVAL
class ESPGRetrieval():
#def __init__(self):
#try:
espg_number = int(findgmlkey(gml_json, '@srsName', 0).split(':')[-1])
#except:
#Warn("Warning: Not a valid ESPG number found")
#return
espg_file = str(espg_number) + '.xml'
def ESPG_retrieve():
if not os.path.isfile('./' + espg_file):
# ESPG XML data URL
espg_url = 'http://epsg.io/' + espg_file
urlreq = URL.Request(
espg_url,
data = None,
headers={
'User-Agent': 'Mozilla/5.0 (Macintosh; Intel Mac OS X 10_9_3) AppleWebKit/537.36 #(KHTML, like Gecko) Chrome/35.0.1916.47 Safari/537.36'
}
)
# Try to download the XML file and save it
with open(espg_file, 'w') as espg_of:
try:
espg_of.write(str(URL.urlopen(urlreq).read().decode('utf-8')))
except:
Warn("Warning: Could not download ESPG metadata")
espg_of.close()
@staticmethod
def ESPG_read():
with open(espg_file, 'r') as espg_rf:
espg_metadata_list = espg_rf.read()
espgparser = GMLDataParser(espg_metadata_list)
espg_json = espgparser.jsonpretty()
gml_datum = findgmlkey(espg_json, 'gml:datumName', 0)
gml_ellipsoid = findgmlkey(espg_json, 'gml:ellipsoidName', 0)
gml_coordsys = findgmlkey(espg_json, 'gml:srsName', 0)
gml_axis_1_abbrev = findgmlkey(espg_json, 'gml:axisAbbrev', 0)
gml_axis_1_dir = findgmlkey(espg_json, 'gml:axisDirection', 0).capitalize()
gml_axis_2_abbrev = findgmlkey(espg_json, 'gml:axisAbbrev', 1)
gml_axis_2_dir = findgmlkey(espg_json, 'gml:axisDirection', 1).capitalize()
# TODO. Have child element #text which contains the actual value
gml_semimajor_axis = findgmlkey(espg_json, 'gml:semiMajorAxis', 0)
gml_inverse_flat = findgmlkey(espg_json, 'gml:inverseFlattening', 0)
#espg_data = ESPGRetrieval()
#espg_data.ESPG_read()
#sys.exit()
################################################################
#
# PHYSICAL AREA SIZE CALCULATOR
def axisCalculator():
# Axis-based data
try:
x_high = float(findgmlkey(gml_json, 'gml:upperCorner', 0).split()[0])
x_low = float(findgmlkey(gml_json, 'gml:lowerCorner', 0).split()[0])
y_high = float(findgmlkey(gml_json, 'gml:upperCorner', 0).split()[1])
y_low = float(findgmlkey(gml_json, 'gml:lowerCorner', 0).split()[1])
except:
# Pixel-based data
try:
x_high = float(findgmlkey(gml_json, 'gml:high', 0).split()[0])
x_low = float(findgmlkey(gml_json, 'gml:low', 0).split()[0])
y_high = float(findgmlkey(gml_json, 'gml:high', 0).split()[1])
y_low = float(findgmlkey(gml_json, 'gml:low', 0).split()[1])
except:
x_high = "Unknown"
x_low = "Unknown"
y_high = "Unknown"
y_low = "Unknown"
for t in (x_high, x_low, y_high, y_low):
if type(t) is not float:
return list([
'Unknown',
'Unknown',
'Unknown',
'Unknown',
'Unknown'
])
def RadtoGrad(num):
rad_to_deg = 180 * num / math.pi
deg_to_grad = 10 * rad_to_deg / 9
return deg_to_grad
###############################
# X and Y lengths
# Area size in km^2
x_length = x_high - x_low
y_length = y_high - y_low
xy_area = (x_length * y_length) / 1000000
###############################
# Inverse geodetic calculation
xy_hypotenuse = math.sqrt(x_length ** 2 + y_length ** 2)
inverse_geod_angle = RadtoGrad(math.atan2(y_length, x_length))
###############################
return list([
format(x_length, '.2f'),
format(y_length, '.2f'),
format(xy_area, '.2f'),
format(xy_hypotenuse, '.2f'),
format(inverse_geod_angle, '.2f')
])
gml_calc = axisCalculator()
################################################################
#
# TFW FORMAT PARSE
def tfwparse():
worldfile_values = gml_posinfo.g + gml_posinfo.gml_pos
worldfile_out = ''
for value in worldfile_values:
worldfile_out += format(float(value)) + '\n'
# Return gml_out, remove last empty line
return worldfile_out[:-1]
################################################################
#
# INFORMATION PARSE
# Extract all important metadata elements
def infoparse():
def getkeys():
# TODO these might or might not be defined in JSON data!
infolist = [
['Image Name', args.inputfile.split('.')[0] ],
['Source Name', findgmlkey(gml_json, '@srsName', 0) ],
['GML File Name', findgmlkey(gml_json, 'gml:fileName', 0) ],
['File Structure', findgmlkey(gml_json, 'gml:fileStructure', 0) ],
['Rectified Grid Coverage ID', findgmlkey(gml_json, '@dimension', 0) ],
#['Axis Names', ' '.join(findgmlkey('gml:axisName', 0)) ],
['Map Scale', ],
['Upper Corner Coordinates', findgmlkey(gml_json, 'gml:upperCorner', 0) ],
['Lower Corner Coordinates', findgmlkey(gml_json, 'gml:lowerCorner', 0) ],
['X-axis Length in Meters', gml_calc[0] ],
['Y-axis Length in Meters', gml_calc[1] ],
['Area Size in Square Kilometers', gml_calc[2] ],
['Distance of Corners Points in Meters', gml_calc[3] ],
['Azimuth Angle of Corner Points in Gradians', gml_calc[4] ],
['Grid Envelope High', findgmlkey(gml_json, 'gml:high', 0) ],
['Grid Envelope Low', findgmlkey(gml_json, 'gml:low', 0) ],
['X-axis Pixel Size in Map Units', gml_posinfo.g[0] ],
['Y-axis pixel size in Map Units', gml_posinfo.g[3] ],
['X-axis Rotation', gml_posinfo.g[1] ],
['Y-axis Rotation', gml_posinfo.g[2] ],
['Upper Left Pixel X-coordinate Center in Map Units', gml_posinfo.gml_pos[0] ],
['Upper Left Pixel Y-coordinate Center in Map Units', gml_posinfo.gml_pos[1] ]
#['EPSG Projection Code',
#['Projection Name',
#['Projection Area',
#['Image Area',
]
#row_format ="{:>15}" * (len(teams_list) + 1)
#print(row_format.format("", *teams_list))
#for team, row in zip(teams_list, data):
# print row_format.format(team, *row)
#for i in infolist:
# print(i)
#sys.exit()
for i in range(len(infolist)):
for j in range(len(infolist[i])):
print(infolist[i][j], end=' ')
print('')
getkeys()
#print(gml_source)
################################
#
# OUTPUT WRITING
try:
args.outputfile
with open(args.outputfile, 'w') as o:
if args.outputformat in 'xml':
if args.formatting in 'pretty':
o.write(gmlparser.xmlpretty())
elif args.formatting in 'raw':
o.write(gmlparser.jsonraw())
else:
raise ValueError("Error: Undefined formatting")
elif args.outputformat in 'json':
if args.formatting in 'pretty':
o.write(gmlparser.jsonpretty())
elif args.formatting in 'raw':
o.write(gmlparser.jsonraw())
else:
raise ValueError("Error: invalid data formatting")
elif args.outputformat in 'tfw' or args.outputformat in 'worldfile':
o.write(tfwparse())
elif args.outputformat in 'info':
o.write(infoparse())
else:
raise ValueError("Error: invalid data format")
except:
if args.outputformat in 'xml':
if args.formatting in 'pretty':
print(gmlparser.xmlpretty())
elif args.formatting in 'raw':
print(gmlparser.xmlraw())
else:
raise ValueError("Error: Undefined formatting")
elif args.outputformat in 'json':
if args.formatting in 'pretty':
print(gmlparser.jsonpretty())
elif args.formatting in 'raw':
print(gmlparser.jsonraw())
else:
raise ValueError("Error: Undefined formatting")
elif args.outputformat in 'tfw' or args.outputformat in 'worldfile':
print(tfwparse())
elif args.outputformat in 'info':
print(infoparse())
else:
raise ValueError("Error: invalid data format")