Hi folks! Is there a way to get a PROJ.4 representation of CRS object?

I've tried to use to_proj4() method:

from rasterio.crs import CRS
wkt = 'PROJCS["WGS 84 / UTM zone 5N",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",-153],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","32605"]]'
crs = CRS.from_wkt(wkt)
crs.to_proj4()

It returns: "+init=epsg:32605" but what I want to get is "+proj=utm +zone=5 +datum=WGS84 +units=m +no_defs". Is it possible?

This is a small piece totally clear of clouds

And this is a larger piece with many clouds

This is the result of the merge using rasterio.merge

I used rasterio.merge. It brings me a result where the most recent acquired image prevails, but don't understand the algorithm method yet. How can I tell the code to keep the image with fewer clouds as possible? Which method are this algorithm using to preserve the earlier raster and not the cleanest (no clouds)?

It would be absolutely awesome if I could select the raster with fewer clouds based in pixel values, for insance, to be merged for this mosaic. Is there a way to choose them?

from affine import Affine
from rasterio.coords import BoundingBox
bbox = BoundingBox(left=30.0, bottom=-29.0, right=31.0, top=-28.00083)
transform = Affine(0.0008333333299814356, 0.0, 16.457916616, 0.0, -0.0008333333300030174, -22.127083043)
window = from_bounds(*bounds, src_pop.transform, height=1200, width=1200)
print(window) Window(col_off=16250.500126164021, row_off=7048.496376568462, width=1200.000004826732, height=1199.0040047916773)

Python 3.8.6 (v3.8.6:db455296be, Sep 23 2020, 13:31:39) 

[Clang 6.0 (clang-600.0.57)] on darwin

Type "help", "copyright", "credits" or "license" for more information.

>>> from rasterio.crs import CRS

>>> CRS.from_epsg(4326)

ERROR 1: PROJ: proj_create_from_database: cannot build geodeticCRS 4326: SQLite error on SELECT extent.description, extent.south_lat, extent.north_lat, extent.west_lon, extent.east_lon, scope.scope, (CASE WHEN scope.scope LIKE '%large scale%' THEN 0 ELSE 1 END) AS score FROM usage JOIN extent ON usage.extent_auth_name = extent.auth_name AND usage.extent_code = extent.code JOIN scope ON usage.scope_auth_name = scope.auth_name AND usage.scope_code = scope.code WHERE object_table_name = ? AND object_auth_name = ? AND object_code = ? ORDER BY score, usage.auth_name, usage.code: no such table: usage

Traceback (most recent call last):

  File "rasterio/_crs.pyx", line 278, in rasterio._crs._CRS.from_epsg

  File "rasterio/_err.pyx", line 190, in rasterio._err.exc_wrap_int

rasterio._err.CPLE_AppDefinedError: PROJ: proj_create_from_database: cannot build geodeticCRS 4326: SQLite error on SELECT extent.description, extent.south_lat, extent.north_lat, extent.west_lon, extent.east_lon, scope.scope, (CASE WHEN scope.scope LIKE '%large scale%' THEN 0 ELSE 1 END) AS score FROM usage JOIN extent ON usage.extent_auth_name = extent.auth_name AND usage.extent_code = extent.code JOIN scope ON usage.scope_auth_name = scope.auth_name AND usage.scope_code = scope.code WHERE object_table_name = ? AND object_auth_name = ? AND object_code = ? ORDER BY score, usage.auth_name, usage.code: no such table: usage

 

During handling of the above exception, another exception occurred:

 

Traceback (most recent call last):

  File "<stdin>", line 1, in <module>

  File "/Users/vincentsarago/Workspace/venv/py38/lib/python3.8/site-packages/rasterio/crs.py", line 333, in from_epsg

    obj._crs = _CRS.from_epsg(code)

  File "rasterio/_crs.pyx", line 280, in rasterio._crs._CRS.from_epsg

rasterio.errors.CRSError: The EPSG code is unknown. PROJ: proj_create_from_database: cannot build geodeticCRS 4326: SQLite error on SELECT extent.description, extent.south_lat, extent.north_lat, extent.west_lon, extent.east_lon, scope.scope, (CASE WHEN scope.scope LIKE '%large scale%' THEN 0 ELSE 1 END) AS score FROM usage JOIN extent ON usage.extent_auth_name = extent.auth_name AND usage.extent_code = extent.code JOIN scope ON usage.scope_auth_name = scope.auth_name AND usage.scope_code = scope.code WHERE object_table_name = ? AND object_auth_name = ? AND object_code = ? ORDER BY score, usage.auth_name, usage.code: no such table: usage

Hi Paolo,

Welcome! I ran the your file through rio-shapes and jq

rio shapes --bidx 1 ~/Downloads/sample_cover.tif | jq -c 'select(.properties.val == 3.0)' | fio collect

and then uploaded to a gist:

https://gist.github.com/sgillies/96ea784a00540d174e02060dde9e1a5a

It looks like the holes aren't lost there. How did you make the second image?

Is that because on the other hand GDAL 3.2.0 and/or PROJ 7.2.0 are heavier than their counterparts from release 1.1.8?

would you expect there to be problems if I were to have pyproj<3.0.0 and rasterio>=1.2 installed in my Python environment, due to their different PROJ versions, or would that still work OK?