1. main@rasterio.groups.io
   Date   
Reprojecting GEOTIFF to epsg:3857 with scrambled CRS data (possibly GCJ-02)

toku2288@...
 

Hello, I am trying to reproject a population distribution raster file of China into 3857 space. I tried but the images are not aligning. Looking at the CRS it has a lot of different standards, I am not sure what to make of it.

The meta of the source file looks something like this:
{'driver': 'GTiff', 
'dtype': 'int32',
 'nodata': 2147483647.0,
 'width': 5049,
 'height': 5582,
 'count': 1,
 'crs': CRS.from_wkt('PROJCS["Albers_Conic_Equal_Area",GEOGCS["GCS_Unknown_datum_based_upon_the_Krassowsky_1940_ellipsoid",DATUM["Not_specified_based_on_Krassowsky_1940_ellipsoid",SPHEROID["Krasovsky_1940",6378245,298.3,AUTHORITY["EPSG","7024"]],AUTHORITY["EPSG","6024"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["standard_parallel_1",25],PARAMETER["standard_parallel_2",47],PARAMETER["latitude_of_center",0],PARAMETER["longitude_of_center",105],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]'),
 'transform': Affine(1000.0, 0.0, -2638055.0614613006,
        0.0, -1000.0, 5922165.484718928)}
PopTIF.crs.data returns:
{'proj': 'aea',
 'lat_1': 25,
 'lat_2': 47,
 'lat_0': 0,
 'lon_0': 105,
 'x_0': 0,
 'y_0': 0,
 'ellps': 'krass',
 'units': 'm',
 'no_defs': True}
after running the reprojection algorithm:
{'driver': 'GTiff',
 'dtype': 'int32',
 'nodata': 2147483647.0,
 'width': 6708,
 'height': 6102,
 'count': 1,
 'crs': CRS.from_dict(init='epsg:3857'),
 'transform': Affine(1184.7466003556895, 0.0, 7551104.353792087,
        0.0, -1184.7466003556895, 7299154.069454485)}




I am then using PopTIF.dataset_mask()  to mask the area of a Heightfield GEOTIFF (I requested the raster area of this file to be the same as the lat/lon of the PopTIF data:  'lat_1': 25,'lat_2': 47, 'lat_0': 0, 'lon_0': 105
 {'driver': 'GTiff', 
'dtype': 'float32',
 'nodata': -32767.0,
 'width': 2816,
 'height': 2288,
 'count': 1,
 'crs': CRS.from_dict(init='epsg:4326'),
 'transform': Affine(0.02197265625, 0.0, 73.4765625,
        0.0, -0.02197265625, 53.7890625)}

 EPSG:4326 to EPSG:3857

{'driver': 'GTiff',
 'dtype': 'float32',
 'nodata': -32767.0,
 'width': 2594,
 'height': 2537,
 'count': 1,
 'crs': CRS.from_dict(init='epsg:3857'),
 'transform': Affine(2655.775407243914, 0.0, 8179373.5227401415,
        0.0, -2655.775407243914, 7130308.041557059)}


So, After moving both GEOTIFFs into 3857 space. I then needed to upscale the Heightfield so the array matched the the size of the Population data:

import skimage.transform as st

pathToRaster  = r'./DEMTIF/cnl10_3857.tif'
newsize = reprojectedmask.shape


with rio.open(pathToRaster) as src:
    arr = src.read()
    arr[arr ==-9999] = np.nan
    arr  = arr[0,:,:]
    new_shape = newsize
    newarr= st.resize(arr, new_shape, mode='constant')

Resulting array can be masked, and looks something like this, Notice how there is no height data in areas such as Taiwan, There has been a mismatch with the projections, However I cant tell where/how. This may be due to the fact I specified the lat/lon in AEA(plus other CRS)? Would appreciate any guidance on how to find the correct area to sample and how to transform it 
Re: Different values when I use a window

Sean Gillies
 

Hi,

We fixed a window shift bug in 1.1.5. See


I wonder if this is another bug or if somehow your installation didn't get the fix. Can you say what your source of the rasterio package is?

On Sun, Jul 5, 2020 at 5:11 PM <adrianocorbelinoii@...> wrote:
I am using rasterio 1.1.5.



--
Sean Gillies
Re: Different values when I use a window

adrianocorbelinoii@...
 

I am using rasterio 1.1.5.
Re: Merge resulting in one pixel shift for half the output file

melda.salhab.14@...
 

Thanks for pointing that out. But it's really odd. I've been trying to update rasterio using conda, but the version remains at v1.1.0. I just tried doing a new install using conda install -c conda-forge rasterio but the version isn't changing. Any tips on how to get the latest version? I'm using a virtual environment 
Re: Merge resulting in one pixel shift for half the output file

Denis Rykov
 

The latest version is 1.1.5 (Jun 3, 2020). 1.1.0 has been released in 2019.
Re: Merge resulting in one pixel shift for half the output file

melda.salhab.14@...
 

Yep - I'm using 1.1.0
Re: Merge resulting in one pixel shift for half the output file

Denis Rykov
 

Do you use the latest version of rasterio?
Merge resulting in one pixel shift for half the output file

melda.salhab.14@...
 

I'm using merge with method = 'max' to merge two rasters with identical metadata. When I open the output file and the input files in QGIS, the top part of the output raster is perfectly merged, but the bottom part of the output raster seems to have the merged values of the pixel to it's left. Very bizarre. Has anyone experienced anything like this? Code below.


    combinedFloodList = [FluvialCropped_outputfile] + [PluvialCropped_outputfile] 
 
    # Iterate over raster files and add them to list in 'read mode'            
    allFiles = []
    for fp in combinedFloodList:
        src = rs.open(fp)
        allFiles.append(src)
 
    # Merge flood files. Using max pixel by pixel method
    floodArray, out_trans = merge(allFiles, method='max')
 
    # Save results as raster
    ## Update the metadata
    out_meta = allFiles[0].meta.copy()
    out_meta.update({"height": floodArray.shape[1], 
                      "width": floodArray.shape[2],
                      "transform": out_trans})
 
    # Output file name
    floodName = f"Flood_{countryCode}_{adm1Index}.tif"
 
    ## export as a new geotiff 
    FloodMerged_outputfile = Flood_outputFolderPath / floodName
    with rs.open(FloodMerged_outputfile, 'w', **out_meta, compress = 'LZW') as dest:
        dest.write(floodArray)
Different values when I use a window

adrianocorbelinoii@...
 

Hello,  I've been having a problem using windows, probably a misunderstanding on my part.
My goal is to dump a specific data from a band into a csv file, first I tried to use a window but the results differ from the results that I see on Qgis(the results appear to be shifted 1 position).
So I decided to try to read the whole band, instead to use a window, and the result match with Qgis.
I think that I am doing something wrong when I create/use the window, what I tried to do was:
  1. create a window from the bounds of the desired area: 
    bandWin = rasterio.windows.from_bounds(*boundingBox,band.transform)
  2. get the transform matrix of that window:
    win_transform = band.window_transform(bandWin)
  3. Obtain the window data:
    win_data = band.read(1, window = bandWin)
  4. Get the desired window indexes:
    indexes = rasterio.transform.rowcol(win_transform, [Xs], [Ys])
  5. Get the values that I want:
    csvData.append(win_data[indexes[0],indexes[1]][0])
The only differences between when I don't use a window is that instead use win_transform I use band.transform and I read the whole band.

All help is appreciated.
List of numpy arrays to raster

ashnair0007@...
 

I originally had a geotiff (size= 4490 x 7341) which I tiled into tiles of size 256 x 256 to perform an operation. The result is a list of 522 numpy array. From the original geotiff I have the geo transform as well.  What would be the recommended way to merge them to get the raster? 
Re: Window from bounds is shifted by 1 Pixel

thomaswanderer@...
 

Thanks for letting me know! Thumbs up for your work!
Re: Window from bounds is shifted by 1 Pixel

Sean Gillies
 

Hi Thomas,

I think Denis found the root of the problem and we have a fix at https://github.com/mapbox/rasterio/pull/1938 that will be a part of the 1.1.5 release. Thank you for the test cases and for your patience!


On Tue, May 26, 2020 at 4:39 PM <thomaswanderer@...> wrote:

I am not 100 sure, but I now simply rounded all windows (origin + offsets) and this seems to have fixed it.



--
Sean Gillies
Re: Create raster aligned with one created by gdal_translate

Denis Rykov
 

Re: can resampling not populate cell if source has nodata?

Nick Forfinski-Sarkozi
 

FYI...I was able to achieve the desired behavior by "lying" about the nodata value.  The actual nodata was nan, but when I specified something else (e.g., nodata=-9999) in the profile, the resampled raster didn't contain values for cells corresponding to source nans.


On Wed, May 27, 2020 at 4:20 PM Sean Gillies via groups.io <sean=mapbox.com@groups.io> wrote:
I see what you mean. Rasterio doesn't have a built in option for this, unless there is a GDAL resampling trick that I am not aware of.

On Wed, May 27, 2020 at 9:55 AM <nick.forfinski@...> wrote:
Hi Sean.  Thanks for the message.  I understand your example, but I may have been confusing in my email.  I think I'm asking something different.  Is there a way to force the resampling process to not populate a resampled value where there are corresponding nodata values in the original raster?  For example, in the image below, the lighter 5-m purple raster is the resampled version of the darker 1-m blue raster.  Many of the resampled raster's pixels correspond to gaps in the original raster.  How could I resample so that all 25 source cells are required to produce a resampled value?  
 
Thanks again for your help and a great module,
Nick
 



--
Sean Gillies
Re: can resampling not populate cell if source has nodata?

Alan Snow
 

Here is an example using rioxarray: https://corteva.github.io/rioxarray/stable/examples/interpolate_na.html

It uses scipy griddata.
Re: can resampling not populate cell if source has nodata?

Sean Gillies
 

I see what you mean. Rasterio doesn't have a built in option for this, unless there is a GDAL resampling trick that I am not aware of.

On Wed, May 27, 2020 at 9:55 AM <nick.forfinski@...> wrote:
Hi Sean.  Thanks for the message.  I understand your example, but I may have been confusing in my email.  I think I'm asking something different.  Is there a way to force the resampling process to not populate a resampled value where there are corresponding nodata values in the original raster?  For example, in the image below, the lighter 5-m purple raster is the resampled version of the darker 1-m blue raster.  Many of the resampled raster's pixels correspond to gaps in the original raster.  How could I resample so that all 25 source cells are required to produce a resampled value?  
 
Thanks again for your help and a great module,
Nick
 



--
Sean Gillies
Re: can resampling not populate cell if source has nodata?

Nick Forfinski-Sarkozi
 

Hi Sean.  Thanks for the message.  I understand your example, but I may have been confusing in my email.  I think I'm asking something different.  Is there a way to force the resampling process to not populate a resampled value where there are corresponding nodata values in the original raster?  For example, in the image below, the lighter 5-m purple raster is the resampled version of the darker 1-m blue raster.  Many of the resampled raster's pixels correspond to gaps in the original raster.  How could I resample so that all 25 source cells are required to produce a resampled value?  
 
Thanks again for your help and a great module,
Nick
 
Re: can resampling not populate cell if source has nodata?

Sean Gillies
 

Hi Nick,

Ignoring nodata values is the default behavior for rasterio, as it is for the GDAL library that rasterio uses. Here's an example using one of rasterio's test files.

$ rio insp tests/data/RGB.byte.tif
Rasterio 1.1.5dev Interactive Inspector (Python 3.6.6)
Type "src.meta", "src.read(1)", or "help(src)" for more information.
>>> show(src.read())

That command uses matplotlib to display the raw data and I've zoomed in to the bottom to show the hard edge along the image collar. Valid data on the inside, nodata on the outside.

Figure_1.png
Now, I'll dilate the image by a factor of 3 on reading and use bilinear resampling to fill in the extra pixels.

>>> from rasterio.enums import Resampling
>>> resampled = src.read(out_shape=(3, src.height * 3, src.width * 3), resampling=Resampling.bilinear)
>>> show(resampled)

Zooming in again to the collar, you can see that the valid data has been smoothed, but the hard nodata edge remains. Nodata values were ignored by the resampling system.
Figure_2.png

On Tue, May 26, 2020 at 11:30 AM <nick.forfinski@...> wrote:
Hi rasterio.  New user here.  Thanks for a fantastic module.  During resampling, is there a way to force the cells of the resampled raster to nodata if any of the corresponding original cells are nodata?   (If by chance anyone is familiar with it, I'm looking for functionality analogous to arcpy's Aggregate function with ignore_data=True https://pro.arcgis.com/en/pro-app/tool-reference/spatial-analyst/aggregate.htm )

Thanks,
Nick



--
Sean Gillies
Re: Get bounds after rasterio merge

Luke
 


On Wed, 27 May 2020 at 04:18, <ashnair0007@...> wrote:
I merged two tiffs via the rasterio.merge function which returned the merged numpy array and the geo transform. Is there a way to get the bounds (in latitude and longitude) of this new merged image?

Re: Window from bounds is shifted by 1 Pixel

thomaswanderer@...
 

I am not 100 sure, but I now simply rounded all windows (origin + offsets) and this seems to have fixed it.
141 - 160 of 698  
More
Home Messages Hashtags Subgroups About Features Pricing Updates Terms Help Twitter