I am resampling a raster to become a different resolution, but it seems like the resolutions are not entirely accurate.

Here is the code I have (the code is very standard, taken from the documentation).

        orig_res = ds.res[0]

        # Calculate scale factor for this specific input layer
        upscale_factor = orig_res / 15.0
        print("upscale factor: ", upscale_factor)

        # Scale while reading in array
        scaled_arr = ds.read(
            out_shape=(
                ds.count,
                int(ds.height * upscale_factor),
                int(ds.width * upscale_factor)
            ),
            resampling=Resampling.bilinear
        )

        # Also have to update the transform to be scaled
        scaled_transform = ds.transform * ds.transform.scale(
            (ds.height / scaled_arr.shape[1]),
            (ds.width / scaled_arr.shape[2])
        )

        scaled_profile = ds.profile
        scaled_profile.update(transform=scaled_transform, height=scaled_arr.shape[1], width=scaled_arr.shape[2])



        with rio.open(f_out_path, 'w', **scaled_profile) as scaled_ds:
            scaled_ds.write(scaled_arr)

            print(f"Orig arr shape: ({ds.count}, {ds.shape[0]}, {ds.shape[1]})")

            print(f"Scaled ds shape: ({scaled_ds.count}, {scaled_ds.shape[0]}, {scaled_ds.shape[1]})")


            print("Orig ds res: ", ds.res)
            print("Scaled ds res: ", scaled_ds.res)

Code output:

upscale factor:  0.3333333333333333
Orig arr shape: (1, 9240, 2671)
Scaled ds shape: (1, 3080, 890)
Orig ds res:  (5.0, 5.0)
Scaled ds res:  (15.0, 15.00561797752809)

But in the final scaled ds resolution, why is the vertical scale direction slightly different? (15.0, 15.00561797752809)?

I'm not sure if it's a bad thing or not. I just want to make sure it's okay, or if there's a way to correct it, that would be good too.

I am attempting to use rasterio to merge 3 bands from Sentinel-2 data, sourced as .JP2 files, into a georeferenced GeoTiff.  However, I am running into 2 problems doing this.

First of all, here is my code for the merge:

    #Setup the bands
    b4 = rio.open(path + Bands.B04.value[0])
    b3 = rio.open(path + Bands.B03.value[0])
    b2 = rio.open(path + Bands.B02.value[0])
   
    #Setup the path for the constructed RGB image
    tif = path + "RGB.tif"
   
    #Create the RGB image
    with rio.open(tif, 'w', driver='Gtiff', width=b4.width, height=b4.height, count=3, crs=b4.crs, tranform=b4.transform, dtype=b4.dtypes[0], photometric="RGB", tile=True, compress='lzw', bigtiff=True) as rgb:
        rgb.write(b2.read(1), 1)
        rgb.write(b3.read(1), 2)
        rgb.write(b4.read(1), 3)
        rgb.close()

Now, how do I copy the source projection information and bounds in the new GeoTiff I am creating?  Second, why is it showing up as all black when viewing the image in most common image browsers, but just fine in mapping software such as ArcMap or ArcGIS Pro?

Any ideas?

No. Exactly as I show in my code example. No need to write it out again, as you are just updating the georeferencing with the dataset opened in update mode ("r+").  The dataset is saved automatically when you exit the `with` statement context.

Luke

What is the best way to do this warping?

Should I not actually warp and just set the `.crs` property of the jp2 raster and then re-save it in the same file?
Or should I warp it the crs, and then re-save as a .jp2 or should I save as a .tif?

Should I use virtual warping like here https://rasterio.readthedocs.io/en/latest/topics/virtual-warping.html
or should I use reprojection like here https://rasterio.readthedocs.io/en/latest/topics/reproject.html 

I just want to make sure I don't mess anything up, so thank you!

@ Luke Pinner,

How should I then save the file with the updated CRS?
Like this? Or is there a better way?

with rasterio.open(in_tif, 'r+') as src:
    src.crs = CRS.from_epsg(32629)
    with raster.open("out.tif", 'w', **src.profile) as dst:
        dst.write(src.read(1))

To override the projection with rasterio:

On lundi 9 mars 2020 14:42:32 CET himat15 via Groups.Io wrote:

Hi, I have a geotiff with a "Transverse_Mercator_Complex" CRS.

Here is the output of the `.crs` property after reading in the tif with

rio.read():

PROJCS["WGS_1984_Complex_UTM_Zone_29N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_
1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0
],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator_Compl
ex"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],P
ARAMETER["Central_Meridian",-9.0],PARAMETER["Scale_Factor",0.9996],PARAMET
ER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]

The rest of my files have the crs EPSG:32629, so I tried to reproject this
tif into EPSG:32629 following the standard given example in the docs here (
https://rasterio.readthedocs.io/en/latest/topics/reproject.html ).

But, I got this error:
No translation for Transverse_Mercator_Complex to PROJ.4 format is known

I found that this projection is the same (?) as the regular Transverse
Mercator from
https://gis.stackexchange.com/questions/226679/complex-utm-projection. (
https://gis.stackexchange.com/questions/226679/complex-utm-projection ) It
seems to just be an ESRI specific one. It's also listed in this repo
https://github.com/Esri/projection-engine-db-doc/blob/master/text/pe_list_p
rojection.txt. (
https://github.com/Esri/projection-engine-db-doc/blob/master/text/pe_list_p
rojection.txt )

What should I do in this case to transform my file in this case
to EPSG:32629?

My understanding is that Transverse_Mercator_Complex corresponds to PROJ
"etmerc" method, which is the default since PROJ 4.9.3 for UTM. PROJ just
lacks a mapping for the Transverse_Mercator_Complex term.

So you could just override the projection information with for example:
gdal_translate in.tif out.tif -a_srs EPSG:32629
(or gdal_edit.py in.tif -a_srs EPSG:32629)
There's no need to do a reprojection.

Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com

Hi, I have a geotiff with a "Transverse_Mercator_Complex" CRS.

Here is the output of the `.crs` property after reading in the tif with rio.read():

PROJCS["WGS_1984_Complex_UTM_Zone_29N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator_Complex"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-9.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]

The rest of my files have the crs EPSG:32629, so I tried to reproject this tif into EPSG:32629 following the standard given example in the docs here.

But, I got this error:
No translation for Transverse_Mercator_Complex to PROJ.4 format is known

I found that this projection is the same (?) as the regular Transverse Mercator from https://gis.stackexchange.com/questions/226679/complex-utm-projection. It seems to just be an ESRI specific one. It's also listed in this repo https://github.com/Esri/projection-engine-db-doc/blob/master/text/pe_list_projection.txt.

What should I do in this case to transform my file in this case to EPSG:32629?

I have some .jp2 files that don't have a `.crs` property after I used `rio.open` on those files. 
I read https://gis.stackexchange.com/a/233417/143989 and after using the `.gcps` property, I was able to get the crs.

But, the second answer also suggested that I could use `gdalwarp` to embed the projection info in the raster directly.
My question is, do I need to do the `gdalwarp`? I don't know if not having the crs in the normal property will affect rasterio's ability to do other operations on the jp2 file raster such as transformations to change the cell size. Or will all rasterio operations work as normal with the jp2 file as it is?

--

--

Add some history for reference: https://rasterio.groups.io/g/dev/topic/switch_from_proj_4_format_to/28991934

Yeah, this of course works. I'm more thinking from a philosophical standpoint of "which object works better to pass around?"


It seems like the pyproj objects work a bit better. And that's what's natively used in geopandas and a bunch of other places. I guess the question becomes "what's the benefit of rasterio having its own CRS object type"...I guess the main argument against is that adding pyproj is adding an additional dependency?

See this gist for a bit of exploration of this.

FYI: https://github.com/corteva/rioxarray/issues/92 and https://github.com/corteva/rioxarray/issues/92