Date   

rasterio opens file from AWS S3 bucket on local machine, but can't find file when deployed to Google App Engine

Judson Buescher
 

Hello,

 

I'm trying to run a python flask app, that opens a file hosted on an S3 bucket and analyzes the data. When I deploy the flask app on my local host, it runs like expected. Once I deploy the app to the cloud, it suddenly cannot find the file. I feel like this might be a bug, but i'm too afraid to post in the github. I'll attach the relevant snippet of code, and the error I'm receiving in the GCP Logs. I apologize for terrible formatting as I'm unaware of how to do it properly in this forum.
Here's the code: 

session = boto3.Session(aws_access_key_id=access_key, aws_secret_access_key=secret_access_key)    

def start_process(map, lon, lat):    
    with rasterio.Env(AWSSession(session)):        
        return getCoordinatePixel(map,lon,lat)  url = 's3://metcon-test-bucket/high_plus.tif'  

app = Flask(__name__) @app.route('/api')
def api_id():    
    lat = int(request.args['lat'])    
    lon = int(request.args['lon'])          
    return jsonify(start_process(url, lon, lat))  

if __name__ == '__main__':    
    app.run()

 and here's the error:



Re: Reprojecting GEOTIFF to epsg:3857 with scrambled CRS data (possibly GCJ-02)

Even Rouault
 

On mardi 21 juillet 2020 18:52:10 CEST toku2288 via groups.io wrote:

> Hi, Just an update. Turns GeoTIFF does use GCJ encoding, However, There seem

> to be some open-source projects online which can reverse the transformation

> to high accuracy, Is there some way to use this algorithm to reproject the

> GeoTIFF file?

 

This had been discussed for PROJ in

https://github.com/OSGeo/PROJ/issues/982

 

Summary: apparently technically doable, but for legal reasons for folks operating in China, cannot be included in PROJ.

 

--

Spatialys - Geospatial professional services

http://www.spatialys.com


Reprojecting output image to input image crs where the image sizes differ

ashnair0007@...
 

I have a geotiff (size= 4490 x 7341) that I cut into 256 x 256 crops to do an operation. Once the operation is done, I merge them back together. Since the original dimensions aren't divisible by 256 the resulting image will be of size 4608 x 7424. However the output tif and the input tif are of the same location, so they need to be geo-referenced the same. Is it possible to do so?


Re: Reprojecting GEOTIFF to epsg:3857 with scrambled CRS data (possibly GCJ-02)

toku2288@...
 

Hi, Just an update. Turns GeoTIFF does use GCJ encoding, However, There seem to be some open-source projects online which can reverse the transformation to high accuracy, Is there some way to use this algorithm to reproject the GeoTIFF file?


Re: Convert BoundingBox data into lat/lon

toku2288@...
 

Hi Guillaume,

Thanks, It worked. I previously tried using the CRS EPSG:3857 using that method, and it returned lat/lon values out of range. Is there some way to get the EPSG:3857 values?


Re: Convert BoundingBox data into lat/lon

Guillaume Lostis
 

Hi,

The rasterio.warp.transform_bounds function should do what you want: https://rasterio.readthedocs.io/en/latest/api/rasterio.warp.html#rasterio.warp.transform_bounds

Pass it the CRS of your dataset as src_crs, CRS({'init': 'EPSG:4326'}) as dst_crs, and your bounds, and you will get your bounds in lon/lat as output.

Regards, 

Guillaume Lostis 

On Tue, 21 Jul 2020, 09:36 toku2288 via groups.io, <toku2288=protonmail.com@groups.io> wrote:
Hi, I need to convert the bounding box data of a TIF file to lat/lon co-ordinates so I can request data from another API

GridTIF.bounds returns a BoundingBox object like the following:
BoundingBox(left=-2638055.0614613006, bottom=340165.4847189281, right=2410944.9385386994, top=5922165.484718928)
Is there some way to convert this to lat/lon of the two corners in it's current CRS?
 


Convert BoundingBox data into lat/lon

toku2288@...
 

Hi, I need to convert the bounding box data of a TIF file to lat/lon co-ordinates so I can request data from another API

GridTIF.bounds returns a BoundingBox object like the following:
BoundingBox(left=-2638055.0614613006, bottom=340165.4847189281, right=2410944.9385386994, top=5922165.484718928)
Is there some way to convert this to lat/lon of the two corners in it's current CRS?
 


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