4

I am new to the area of Remote Sensing. I am working with the Sentinel-2 data from AWS. I've downloaded this specific tile and created an RGB GeoTiff. Now, I want to crop the sea in the bottom of the image. First, I get a specific longitude and latitude (here: 42.364456, 0.183105). Since gdalinfo tells me that the coordinate system is "WGS 84 / UTM zone 30N" I converted the longitude and longitude to UTM which results in the zone 31T. But this does not work with the coordinates of the image (it is tile 30T, I already know that the Sentinel tiles are overlapping) How can I solve this problem? Longitude 42.364456 and latitude 0.183105 result in UTM 31T 268039 4694087.

In general: How can I find the exact position (longitude and latitude given) in a tile given like the example here?

Here is the output of gdalinfo of the example tile:

Driver: GTiff/GeoTIFF
Files: 30_T_YN_2015_12_3.tif
Size is 10980, 10980
Coordinate System is:
PROJCS["WGS 84 / UTM zone 30N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-3],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","32630"]]
Origin = (699960.000000000000000,4800000.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=JPEG
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  699960.000, 4800000.000)
Lower Left  (  699960.000, 4690200.000)
Upper Right (  809760.000, 4800000.000)
Lower Right (  809760.000, 4690200.000)
Center      (  754860.000, 4745100.000)
Band 1 Block=10980x16 Type=Byte, ColorInterp=Red
Band 2 Block=10980x16 Type=Byte, ColorInterp=Green
Band 3 Block=10980x16 Type=Byte, ColorInterp=Blue

Tile: enter image description here

  • The easiest way to do this might be to create a vector layer of the lake and then use gdalwarp with the -cutline options (see this post). GDAL will handle the translation between different projections for you in that case. – Kersten Dec 5 '16 at 13:56
  • Thanks Kersten. However, I need a general solution. I will have a lot of addresses, which will be mapped to their coordinates (longitude and latitude) and then some type of cropping is performed. So I am not really interested in this lake but in a general solution of the problem. Creating a bounding box which just works for one certain are does not work for me. Sorry for the inaccurateness. – Elternhaus Dec 5 '16 at 14:01
  • Clarification: I am (more) interested in mapping the point (given by longitude and latitude) to the right position in the image. Sorry for the inaccurateness. – Elternhaus Dec 5 '16 at 14:08
  • I think that you can classify the image (water) and then extract the water polygons. After that you can mask all the lakes from the image. – Pau Dec 5 '16 at 15:50
  • Thanks Pau. However, please see my comments above. I am interested in a general solution. How can I find the position in the image which belongs to the given longitude and latitude? – Elternhaus Dec 5 '16 at 16:17
1

Here is a python script that computes that.

from functools import partial
from pprint import pprint

import pyproj
import rasterio

from shapely.geometry import Point
from shapely.ops import transform

# gdalinfo for sample image returns:
#Size is 10980, 10980
#Corner Coordinates:
#Upper Left  (  600000.000, 9900040.000) ( 56d 6' 4.50"W,  0d54'15.32"S) 
#Lower Left  (  600000.000, 9790240.000) ( 56d 6' 3.14"W,  1d53'51.08"S)
#Upper Right (  709800.000, 9900040.000) ( 55d 6'52.88"W,  0d54'13.95"S)
#Lower Right (  709800.000, 9790240.000) ( 55d 6'50.03"W,  1d53'48.20"S)
#Center      (  654900.000, 9845140.000) ( 55d36'27.64"W,  1d24' 2.33"S)
img = rasterio.open('S2B_21_M_XU_2017_8_23_0_B04.tif')

# We will test with scene center as a point in EPSG 4236 (lat/lon, WGS84)
longitude = -55.60768
latitude = -1.400647
input_coord = Point(longitude, latitude)

to_target_crs = partial(pyproj.transform, pyproj.Proj('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs '),
                    pyproj.Proj(img.crs))
output_coord = transform(to_target_crs, input_coord)

print "Coord in image projection: %f, %f" % (output_coord.x, output_coord.y)

# Affine transformation from image plain coordinates to image column, line
aff = ~img.transform

print "Column, line in image: %f, %f" % ((output_coord.x, output_coord.y) * aff)

Should give you as result:

Coord in image projection: 654899.705524, 9845140.101523
Column, line in image: 5489.970552, 5489.989848
0

The easiest way to crop a raster to certain bounds (given in an arbitrary coordinate system like {lat, lon}) is to use gdal_translate. You can specify the bounding box using -projwin and the coordinate system of that bounding box using -projwin_srs.

gdal_translate -projwin lon_min lat_max lon_max lat_min -projwin_srs EPSG:4326 \
input_data.tif output_data.tif

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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