1

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.

4
  • Welcome to GIS SE. As a new user, please take the Tour, which emphasizes the importance of asking One question per Question. Please place the relevant code block in the question body (formatted with the {} button) -- external links change and disappear, damaging the question; besides, if you place it in the question it's much more likely to be read.
    – Vince
    Apr 9, 2019 at 13:28
  • 1
    why not use gdal ? with gdal_translate, the -projwin option can define the coordinates of you geotif in the same ccordinate system as the output (which is straightforward , and it will reproject your data at the same time (using -t_srs EPSG:21781
    – radouxju
    Apr 9, 2019 at 15:05
  • Please post the answer to your question as an answer, instead of an edit, and accept it. This way the question also appears as answered for people searching for this problem in the future.
    – Kersten
    Apr 29, 2019 at 8:44
  • Exuse me, I wasn't aware of the self-answering/accepting mechanism. Done.
    – bebissig
    Apr 29, 2019 at 8:56

1 Answer 1

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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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