2

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.

3

There is no projection problem here

  1. Shapely uses a cartesian plane system for computing geometries and the result is an unary_union of the original geometry (without any reprojection: "Shapely does not support coordinate system transformations. All operations on two or more features presume that the features exist in the same Cartesian plane", Shapely User manual)
  2. Fiona and the script you use for dissolving your shapefile (Dissolve polygons based on attributes with Python (shapely - fiona)) gives the same crs as the original shapefile

The official EPSG: definition is

+proj=omerc +lat_0=57 +lonc=-133.6666666666667 +alpha=323.1301023611111 +k=0.9999 +x_0=5000000 +y_0=-5000000 +no_uoff +gamma=323.1301023611111 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

But the content of the ESRI .prj file (EPSG:26931) used by Fiona or ogr has no alpha value:

PROJCS["NAD83_Alaska_zone_1",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Hotine_Oblique_Mercator_Azimuth_Natural_Origin"],PARAMETER["latitude_of_center",57],PARAMETER["longitude_of_center",-133.6666666666667],PARAMETER["azimuth",323.1301023611111],PARAMETER["scale_factor",0.9999],PARAMETER["false_easting",5000000],PARAMETER["false_northing",-5000000],UNIT["Meter",1]]

If you have any doubts about that, simply replace the content of the resulting watersheds.prj by the content of the original.prj file.

If you want to change the projection of a geometry, you need to use pyproj or gdal.osr and there are many examples in GISse.

  • Thanks, My issues is that my full script is supposed to provide easily usable file so I can't afford to have a 10000x and 26931 so I guess I just have to replace the CRS definition by the one recognized by QGIS. Thanks for the explanation – User18981898198119 Jun 9 '15 at 17:35
  • The .prj file is recognized by QGIS – gene Jun 9 '15 at 17:53
2

There has been a lot of trouble with +omerc, +gamma and +no_uoff parameters in PROJ.4 and GDAL.

Latest tickets are

https://trac.osgeo.org/gdal/ticket/4910

https://trac.osgeo.org/gdal/ticket/5511

https://trac.osgeo.org/proj/ticket/114

I'm not sure if your GDAL and fiona versions are working with the latest releases of PROJ and GDAL. It seems that at least latest GDAL builds do the projection correctly by now.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.