0

I am researching to complete a project.

The aim is by given a coordinate I can extract the image among the raster files.

How do I check if the coordinate that I have for example (51.3334198, 3.2973934) is in a raster image file - k_01.tif?

If this coordinate is indeed in k_01.tif, how do I extract a small part of if i.e MxN window?

My code:

import rasterio      
src = rasterio.open('k_01.tif')
src.bound

>BoundingBox(left=145000.0, bottom=238000.0, right=162000.0, top=247000.0)


src.crs


>CRS.from_epsg(31370)

@snowman2 gave the link Extracting data from a raster.

I followed the advice and my code as follow:

I used another .tif and another longitude and latitude.


import rasterio
from pyproj import Transformer

lon = 4.2811439
lat = 51.2010511


path = './k15.tif'


with rasterio.open(path) as rds:
    # convert coordinate to raster projection
    transformer = Transformer.from_crs("EPSG:4326", rds.crs, always_xy=True)
    xx, yy = transformer.transform(lon, lat)


# get value from grid
value = list(rds.sample([(xx, yy)]))[0]




value 
> array([11.61], dtype=float32)


xx , yy


>(143876.25575738482, 210123.58613220416)


src.bounds


> BoundingBox(left=130000.0, bottom=198000.0, right=162000.0, top=218000.0)

enter image description here

From the result above I understand that xx,yy is the conversion of the latitude and longitude and it is indeed in the boundary of the tif.

What is the value 11.6 for?

My second part of the question is how to extract the latitude and longitude with MxN window ?

As shown as picture below.

enter image description here

And crop it to just as follow?

enter image description here

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
Adam Tan
  • 9
  • 3

1 Answers1

1

I think what you are missing for is a combination of:

  1. index
  2. Window.from_slices. Also see windowed rw docs.

This should do it:

import rasterio
from rasterio.windows import Window
from pyproj import Transformer

lon = 4.2811439
lat = 51.2010511


path = './k15.tif'


with rasterio.open(path) as rds:
    # convert coordinate to raster projection
    transformer = Transformer.from_crs("EPSG:4326", rds.crs, always_xy=True)
    xx, yy = transformer.transform(lon, lat)
    row, col = rds.index(xx, yy)
    window = Window.from_slices(rows=(row-2, row+2), cols=(col-2, col+2))
    data = rds.read(window=window)
snowman2
  • 7,321
  • 12
  • 29
  • 54