3

I'm trying to merge some rasters with rasterio.merge.merge. They are in the same CRS, EPSG:3857. I want the merged raster to cover the whole earth. Exactly the range (in coordinates): xmin, ymin, xmax, ymax = -180.0225, -90.0225, 179.9775, 90.0225.

The problem is that when using the bounds parameter of this function I have to give these limits in the same CRS of the rasters. How can I convert the tuple xmin, ymin, xmax, ymax from regular coordinates to projected coordinates in EPSG:3857?

1
  • For EPSG:4326 or OGC::CRS84 the whole Earth extents are limited to +/-180 for longitude and +/-90 for Latitude
    – nmtoken
    Commented May 19, 2023 at 9:39

3 Answers 3

4

For reasons discussed here and with the introduction of rasterio.warp.transform_bounds I think the way to do this would be:

from rasterio.warp import transform_bounds

bounds = (-180.0225, -90.0225, 179.9775, 90.0225)

xmin, ymin, xmax, ymax = transform_bounds(sar_img_crs.to_epsg(), 4326, *bounds)

Also, the answer by Pierrick above needs the always_xy=True flag to be set in the transformer.transform call for it to return values in the correct order.

1
  • 1
    The correct order for EPSG:4326 would be -90, -180, 90, 180
    – nmtoken
    Commented May 19, 2023 at 9:27
4

This is how to achieve CRS transformation using rasterio (and optionally shapely):

import rasterio.warp
from rasterio.crs import CRS

# In GeoJSON format
x_min, y_min, x_max, y_max = bounds
feature = {
    "type": "Polygon",
    "coordinates": [
       [[x_max, y_min], 
        [x_max, y_max], 
        [x_min, y_max], 
        [x_min, y_min], 
        [x_max, y_min]]
   ]
}

# Project the feature to the desired CRS
feature_proj = rasterio.warp.transform_geom(
    CRS.from_epsg(4326),
    CRS.from_epsg(3857),
    feature
)

With the latest version of rasterio (tested with 1.2+) you can also use shapely, which simplify the code:

import rasterio.warp
from rasterio.crs import CRS
from shapely.geometry import box

# Project the feature to the desired CRS
feature_proj = rasterio.warp.transform_geom(
    CRS.from_epsg(4326),
    CRS.from_epsg(3857),
    box(*bounds)
)
2
1

simply transform your bb as 2 corners using pyproj:

from pyproj import CRS
from pyproj import Transformer

transformer = Transformer.from_crs(4326, 3857)

xmin, ymin = transformer.transform(-180.0225, -90.0225, always_xy=True)
xmax, ymax = transformer.transform(179.9775, 90.0225, always_xy=True)
2

Your Answer

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

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