Search code examples
python-3.xgisrastergdalmap-projections

How can i preserve the source raster projection when using gdal_translate?


I'm currently working on a small raster refining tool. The goal is, to have a simple CLI tool, to compute tiles from a georeferenced source raster and create a corresponding index.shp. For this I'm using python 3.7 and gdal. The tool runs smoothly and generates the expected tiles and shapefile, but it gets rid of the projection, which is stored in the source raster. Qgis defaults the newly computed tiles to EPSG 4326 while informing me about an unknown projection. The original raster is in EPSG 25832.

My Setup:

Windows 10 64 bit

Python 3.7.2

Gdal I cannot access the specific version, since gdal-config is not installed and I cannot make it work, but it is 64-bit and I installed it through the binaries provided on gisinternals.com. Windows software list says GDAL 204 MSVC 2017.

While running the script, I get error messages telling me about missing files, e.g. pcs.csv, datum.csv ellipsoid.csv and so on. This indicates that having those files, would fix my problem. But oddly enough, I have used Osgeo4W to install, python 2.7 with gdal and it works like a charm, of course having adjusted the python parts. Tiles get calculated and stay in the projection of the source. Without any external files which specify a projection, in fact using the exact same data which is really confusing to me.

To my understanding, there is no flag or option which forces gdal to keep the projection. If have overlooked or missunderstood the docs im glad for advice.

Before anyone asks, i know that using the osgeo4w installer is obviously the easy and working solution here. But keeping in mind that python 2.7 will soon be discontinued and also using this as a chance to learn new things i wanted to build a 3.7 based tool with gdal installed on my machine

The corresponding code looks like this and does the following :

1.) Command string is build

2.) string is handed to os.system, which in turn executes accordingly

    for i in range(0, width, tilelenght):
        y = 0
        for j in range(0, height, tilelenght):
            gdaltranString = f'gdal_translate -of GTIFF -srcwin {i}, {j}, {tilelenght}, {tilelenght} {input_filepath} {output_filepath}{x}_{y}.tif'
            subprocess.run(gdaltranString)
            y = y+1
        x = x+1

The expected result, would be a collection of functional .tif files which have the EPSG code of the source file, in this case 25832.

But as already mentioned, the projection gets lost somewhere in the process.


Solution

  • So,i have found the solution to my problem, without really understanding how it became an issue to begin with.

    The solution was to create an user variable GDAL_DATA with the path to the projection definition files.

    The weird thing is, i now have GDAL_DATA, as system variable and user variable, both pointing to the same directory.

    If someone knows more about the mysterious ways of windows system variables, please share your wisdom, or the source of said wisdom.