Change CRS of large GeoTIFF and extract windows from large raster in Python

Question

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.

Goal/Issues:

(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.

Answer 1

Answer: With help of your answers I managed to achieve what I needed. For reference, here the relevant steps. Data is loaded, a window with given bounds in target projection is extracted. The window data is saved in an output file (notably, the CRS of the extracted data still seems to be the source's). Then I reload the file and create meshgrids for eastings and northings (suggested here: Obtain coordinates and corresponding pixel values from GeoTiff using python gdal and save them as numpy array)

raster = gdal.Open(r'C:\Users\Lenovo\Dropbox\Code\avalanche\data\DSM\copernicus_eu_dem_v11_E40N20\eu_dem_v11_E40N20\eu_dem_v11_E40N20.TIF')

bounds=[center[0]-width/2, center[1]+width/2,center[0]+width/2, center[1]-width/2]

gdal.Translate('output.tif',raster,projWin=bounds,outputBounds=bounds,projWinSRS='EPSG:21781')

with rasterio.open('output.tif') as r:
    np.warnings.filterwarnings('ignore')
    T0 = r.transform
    elev_grid = r.read()[0]

cols, rows = np.meshgrid(np.arange(elev_grid.shape[1]),np.arange(elev_grid.shape[0]))

T1 = T0*Affine.translation(0.5,0.5)
rc2en = lambda r,c: (c,r)*T1
east_grid, north_grid = np.vectorize(rc2en, otypes=[np.float, np.float]))rows,cols)