This is the first time I am working with GeoTIFF formats and I want to do this in Python. I am trying to use rasterio and pyproj for the following tasks.
I have already made some first steps based on numerous Stackexchange questions. Now I am stuck at a point, where I am lost - too many suggestions how to approach this - which appear to be close to my problem, but not close enough. I could possibly copy-paste more together, but I am afraid that I am not approaching a clean solution.
Therefore I share a jupyter notebook via GitHub where you can see what I have already done and where I formulate the questions right where they come up. https://github.com/bebissig/geotiff_basic_manipulations/blob/master/copernicus_preprocess.ipynb
Overview: I have a GeoTIFF (~5GB) from the Copernicus database (https://land.copernicus.eu/imagery-in-situ/eu-dem/eu-dem-v1.1?tab=download) that contains elevation data of central Europe and need to achieve the following goals.
(These are the main goals, specific questions are in the notebook.)
First; The data comes in EPSG:3035 projection. I need to extract data from Switzerland in projection EPSG:21781. So in my thinking I would be "nice" to translate the whole raster to EPSG:21781. How to do that, do I only have to change the .crs accordingly? Or also manually change the .transform? Or would you suggest to approach this totally differently?
Second; Say I have coordinates of a point in the destination projection (EPSG:21781). I want to extract a window around such a given point say coords+-2000m (in vertical and horizontal direction). Ideally, this would be still a GeoTiff, such that it still contains x,y coordinates in that window (in destination projeciton EPSG:21781). (This is also in view that possibly, later I need to add other data-bands from another database so, I guess I want to stick to a GeoFormat as long as possible).
That's actually it. If I had these windows I could read the data to np-arrays and go on with processing.
--> Se my answer below for a snippet how I made it work.