search

Toblerity / Fiona

Fiona reads and writes geographic data files
https://fiona.readthedocs.io/
BSD 3-Clause "New" or "Revised" License
1.16k stars 203 forks source link

Importing fiona causes subprocess calls to gdalwarp to fail to do datum grid shifting properly #439

Closed JesseCrocker closed 7 years ago

JesseCrocker commented 7 years ago

Expected behavior and actual behavior.

Calling subprocess.call("gdalwarp -t_srs") to warp a file from nad27 based projection to a wgs84 based projection doesn't properly do the datum grid shifting after calling import fiona.

It seems like importing fiona is doing something to env that is causing gdalwarp to fail to load the proj grid shift files.

Steps to reproduce the problem.

https://gist.github.com/JesseCrocker/d43016ca0b66ba8824da9254758ab048

source_image = "source.tiff"
subprocess.call("gdalwarp -t_srs epsg:3857 -r cubic " + source_image + " pre-import.tiff", shell=True)
import fiona
subprocess.call("gdalwarp -t_srs epsg:3857 -r cubic " + source_image + " post-import.tiff", shell=True)

Test file that goes with the script https://www.dropbox.com/s/bc1klyb73ep074r/source.tiff?dl=0

post-import.tiff file overlayed on pre-import.tiff file with 50% opacity showing the shift.

screen shot 2017-04-14 at 11 41 19 am

Projection info for source

$ gdalsrsinfo source.tiff

PROJ.4 : '+proj=longlat +datum=NAD27 +no_defs '

OGC WKT :
GEOGCS["NAD27",
    DATUM["North_American_Datum_1927",
        SPHEROID["Clarke 1866",6378206.4,294.9786982138982,
            AUTHORITY["EPSG","7008"]],
        AUTHORITY["EPSG","6267"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4267"]]

projection info in the same for both result files

$ gdalsrsinfo pre-import.tiff

PROJ.4 : '+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs'

OGC WKT :
PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],
    AUTHORITY["EPSG","3857"]]

Output from test case with --debug on added to gdal_warp call

GDAL: GDALOpen(source.tiff, this=0x7f88064373f0) succeeds as GTiff.
OGRCT: PROJ >= 4.8.0 features enabled
OGRCT: Using locale-safe proj version
OGRCT: Wrap source at -109.469.
OGRCT: Source: +proj=longlat +datum=NAD27 +no_defs
OGRCT: Target: +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs
OGRCT: Wrap target at -109.469.
OGRCT: Source: +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs
OGRCT: Target: +proj=longlat +datum=NAD27 +no_defs
OGRCT: Wrap source at -109.469.
OGRCT: Source: +proj=longlat +datum=NAD27 +no_defs
OGRCT: Target: +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs
OGRCT: Wrap target at -109.469.
OGRCT: Source: +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs
OGRCT: Target: +proj=longlat +datum=NAD27 +no_defs
Creating output file that is 574P x 540L.
GDAL: GDALDriver::Create(GTiff,pre-import.tiff,574,540,4,Byte,0x0)
Processing input file source.tiff.
WARP: Copying metadata from first source to destination dataset
Using band 4 of source image as alpha.
GDAL: GDALDefaultOverviews::OverviewScan()
Using band 4 of destination image as alpha.
GDAL: GDAL_CACHEMAX = 1638 MB
GDAL: GDALWarpKernel()::GWKGeneralCase()
Src=0,0,658x434 Dst=0,0,574x540
0...10...20...30...40...50...60...70...80...90...100 - done.
GDAL: GDALClose(source.tiff, this=0x7f88064373f0)
GDAL: GDALClose(pre-import.tiff, this=0x7f8806500290)
GDAL: GDALOpen(source.tiff, this=0x7fad8fc37410) succeeds as GTiff.
OGRCT: PROJ >= 4.8.0 features enabled
OGRCT: Using locale-safe proj version
OGRCT: Wrap source at -109.469.
OGRCT: Source: +proj=longlat +datum=NAD27 +no_defs
OGRCT: Target: +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs
OGRCT: Wrap target at -109.469.
OGRCT: Source: +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs
OGRCT: Target: +proj=longlat +datum=NAD27 +no_defs
OGRCT: Wrap source at -109.469.
OGRCT: Source: +proj=longlat +datum=NAD27 +no_defs
OGRCT: Target: +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs
OGRCT: Wrap target at -109.469.
OGRCT: Source: +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs
OGRCT: Target: +proj=longlat +datum=NAD27 +no_defs
Creating output file that is 574P x 540L.
GDAL: GDALDriver::Create(GTiff,post-import.tiff,574,540,4,Byte,0x0)
Processing input file source.tiff.
WARP: Copying metadata from first source to destination dataset
Using band 4 of source image as alpha.
GDAL: GDALDefaultOverviews::OverviewScan()
Using band 4 of destination image as alpha.
GDAL: GDAL_CACHEMAX = 1638 MB
GDAL: GDALWarpKernel()::GWKGeneralCase()
Src=0,0,658x434 Dst=0,0,574x540
0...10...20...30...40...50...60...70...80...90...100 - done.
GDAL: GDALClose(source.tiff, this=0x7fad8fc37410)
GDAL: GDALClose(post-import.tiff, this=0x7fad8fc3c3a0)

Operating system

Tried on macOS 10.12.4 and ubuntu 16.04

Fiona version and provenance

Testing on macOS was with fiona version '1.7.1.post1' and 1.7.5 installed from pip, Fiona-1.7.5-cp27-cp27m-macosx_10_6_intel.whl

geowurster commented 7 years ago

Confirmed against master. Interestingly, importing Rasterio instead of Fiona produces the same result.

JesseCrocker commented 7 years ago

After some additional debugging, I've found the cause of this, setting os.environ['PROJ_LIB'] to a non existing directory at https://github.com/Toblerity/Fiona/blob/master/fiona/_drivers.pyx#L120

And the workaround is to set PROJ_LIB to a valid value(usually /usr/local/share/proj/) before importing fiona.

JesseCrocker commented 7 years ago

Maybe a better title for this issue would be "importing fiona cause PROJ_LIB environment variable to be set to a non-existing directory", but maybe im not understanding the cause correctly, so not renaming yet