I have a really annoying problem working with QGIS for visualization and python (fiona / OGR) for the actual operations. I am, by example, doing a simple dissolving operation with fiona. The input is a EPSG : 26931 defined by :
+proj=omerc +lat_0=57 +lonc=-133.6666666666667 +alpha=323.1301023611111 +k=0.9999 +x_0=5000000 +y_0=-5000000 +gamma=323.1301023611111 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
But then my output systematically switch to "Generated CRS" 10000X defined by :
+proj=omerc +lat_0=57 +lonc=-133.666666667 +alpha=323.130102361 +k=0.9999 +x_0=5000000 +y_0=-5000000 +no_uoff +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
So obviously fiona ignore the gamma value. My data being really huge I would like to avoid to correct it with gdalwarp every time as it is really long.
Here is an example of code that modify the CRS : def dissolve(base_watershed):
from shapely.geometry import shape, mapping
from shapely.ops import unary_union
import fiona ,sys, os
import itertools
from osgeo import ogr
watersheds = "watersheds.shp"
#from http://gis.stackexchange.com/questions/149959/dissolve-polygons-based-on-attributes-with-python-shapely-fiona
with fiona.open(base_watershed) as input:
# preserve the schema of the original shapefile, including the crs
meta = input.meta
with fiona.open(watersheds, 'w', **meta) as output:
# groupby clusters consecutive elements of an iterable which have the same key so you must first sort the features by the 'HUC_8' field
e = sorted(input, key=lambda k: k['properties']['HUC_8'])
# group by the 'HUC_8' field
for key, group in itertools.groupby(e, key=lambda x:x['properties']['HUC_8']):
properties, geom = zip(*[(feature['properties'],shape(feature['geometry'])) for feature in group])
# write the feature, computing the unary_union of the elements in the group with the properties of the first element in the group
output.write({'geometry': mapping(unary_union(geom)), 'properties': properties[0]})
Some other part of my code using OGR does the same thing. Is there an easy solution or do I really have to do the full process and reproject everything at the end? I am wondering if overwriting the 10000X CRS definition by the correct 26931 would be an issue since it looks like there is no transformation between the two.
