Date   

Automate process of creating a mask

ts@...
 

Dear Rasterios,

beeing quite new to rasterio I'm looking for help regarding the following:

I've a large amount of 1 Channel Satellite imaginary. While the data contains no-data information around the actual motif nearly black areas do exist which I'd like to cleanup.


Following the docs masking by shapefile  It perfectly works to crop my image to the desired image area:



The red shape is quite inaccurate drawn just for the purpose of testing!

My question now is, how can I automate the creation of a shapefile to mask with based on the grey value or a defined border within my data? Having a closer look  on the edge of my desired image area my problem should get clear. The border is fuzzy, further the nearly black pixel value could exists in the desired motif.



I'm thankful if you could show me a direction to automate the creation of a shape mask (or geojson which should work as well from what I read at docs)!

– Toni


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

Sean Gillies
 

Hi Judson,

On Wed, Aug 12, 2020 at 9:42 AM Judson Buescher <judson.buescher@...> wrote:

Hi Sean,

The only way that I’m installing the package is through my requirements.txt which I include in my deploy folder. I have

 

rasterio==1.1.5

 

in the requirements.txt. For the credentials I’ve tried a couple different ways. I’ve tried setting them using

 

os.environ['GS_SECRET_ACCESS_KEY'] = secret_access_key

os.environ['GS_ACCESS_KEY_ID'] = access_key

 

I’ve also tried

 

session = boto3.Session(aws_access_key_id=access_key,

                        aws_secret_access_key=secret_access_key)

 

and then using that session object when I open the raster file. Was there any other info you’d like?

Thanks a ton for the help,
Judson


If the source data for your App Engine app is on AWS, you'll need to make sure that you're providing AWS keys to access it.  You'll want

    AWS_SECRET_ACCESS_KEY = secret_access_key
    etc.

in your environment, not GS_SECRET_ACCESS_KEY (etc.). I've never used boto3 on App Engine and wouldn't recommend that approach of configuring AWS for GDAL and rasterio.

--
Sean Gillies


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

judson.buescher@...
 

Hi Sean,

The only way that I’m installing the package is through my requirements.txt which I include in my deploy folder. I have

 

rasterio==1.1.5

 

in the requirements.txt. For the credentials I’ve tried a couple different ways. I’ve tried setting them using

 

os.environ['GS_SECRET_ACCESS_KEY'] = secret_access_key

os.environ['GS_ACCESS_KEY_ID'] = access_key

 

I’ve also tried

 

session = boto3.Session(aws_access_key_id=access_key,

                        aws_secret_access_key=secret_access_key)

 

and then using that session object when I open the raster file. Was there any other info you’d like?

Thanks a ton for the help,
Judson

 

From: <main@rasterio.groups.io> on behalf of "Sean Gillies via groups.io" <sean.gillies@...>
Reply-To: "main@rasterio.groups.io" <main@rasterio.groups.io>
Date: Saturday, August 8, 2020 at 10:43 PM
To: "main@rasterio.groups.io" <main@rasterio.groups.io>
Subject: Re: [rasterio] rasterio opens file from AWS S3 bucket on local machine, but can't find file when deployed to Google App Engine

 

Apologies for the late reply. There may be a bug, yes, but there are more likely sources of trouble. Can you explain more about how you are installing the rasterio package on GCP and how you are configuring your AWS credentials?

 

On Thu, Jul 23, 2020 at 7:06 AM <judson.buescher@...> wrote:

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:



--

Sean Gillies


Re: window_transform doesn't produce accurate coordinate transforms

pranay1117@...
 

Hi Sean, thanks for your reply - that was indeed the source of the issue on my end. I was able to fix it!

--
Best, 
Pranay


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

Sean Gillies
 

Hi,

On Tue, Jul 21, 2020 at 11:52 PM <ashnair0007@...> wrote:
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?

If your original raster's georeferencing did not involve rotation, and if the merged raster has the same origin (upper left corner) as the original, you can reuse the original raster's transform property. If the merged raster has a different origin, you must compute a new transform using https://rasterio.readthedocs.io/en/latest/api/rasterio.transform.html#rasterio.transform.from_bounds or https://rasterio.readthedocs.io/en/latest/api/rasterio.transform.html#rasterio.transform.from_origin.

Sorry about the late reply. I hope it is still useful.

--
Sean Gillies


Re: window_transform doesn't produce accurate coordinate transforms

Sean Gillies
 

Hi,

On Thu, Jul 23, 2020 at 7:06 AM pranay1117 via groups.io <pranay1117=yahoo.co.in@groups.io> wrote:
Hi, I'm new to rasterio and working with geographical data. So far, it's been a great experience working with this library - but I've run into an issue when working with the Window module. My use case is as follows: 
- I obtain 2 bounding coordinates from the user (top-left and bottom-right) and convert them to (row, col) using `source.index(lon, lat)`.
- I create a window using these bounds and read the dataset from that. Also, I get the appropriate transform for this window.
selected_window = Window.from_slices((row0, row1), (col0, col1))
data = src.read(1, window=selected_window)
affine_transform = src.window_transform(selected_window)

Afterwards, I am using this `data` numpy array to select various points in these bounds and plot them on a map. However, the (lat, long) coordinates obtained using the affine_transform above are not very accurate. When I run the following:
print("Lat/Lon bounds: {}".format(src.window_bounds(selected_window)))
print("Transform at (0,0): {}".format(affine_transform * (0,0)))
r, c = data.shape
print("Transform at {}: {}".format((r-1, c-1), affine_transform * (r-1, c-1)))

The output is as follows:
Source lat/lon bounds: (72.45 20.69, 82.39, 28.78)
Transform at (0,0): (72.42, 28.78)
Transform at (970, 1193): (80.53, 18.841)

As is evident from the output, the latitude must lie between (20.69,28.78) but the transform value 18.841 at (970, 1193) is not in that interval. 

I saw some earlier posts in this group that suggested updating to version 1.1.5. I've done that but the issue persists. 
Any help will be greatly appreciated. Thanks!

Your usage of the affine transform

    affine_transform * (r-1, c-1)

is incorrect. You must swap the order of the tuple indexes to get the correct results because mathematically speaking

    (x', y') = affine_transform * (x, y)

and row is y and column is x. Rasterio has a method that will help keep this straight: https://rasterio.readthedocs.io/en/latest/api/rasterio.transform.html#rasterio.transform.xy.

The ordering of items in (row, col) versus (x, y) is a constant source of confusion. I hope this reply hasn't come too late and that it is useful.

--
Sean Gillies


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

Sean Gillies
 

Apologies for the late reply. There may be a bug, yes, but there are more likely sources of trouble. Can you explain more about how you are installing the rasterio package on GCP and how you are configuring your AWS credentials?

On Thu, Jul 23, 2020 at 7:06 AM <judson.buescher@...> wrote:

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:




--
Sean Gillies


Re: returning subset of raster data via API

Sean Gillies
 

On Wed, Aug 5, 2020 at 5:41 AM Tom Kralidis <tom.kralidis@...> wrote:

Per https://gdal.org/user/virtual_file_systems.html#drivers-supporting-virtual-file-systems and confirmed by Even R, GDAL's NetCDF driver cannot write to /vsi*. Perhaps a workaround can be to drop to a tempfile approach when writing NetCDF.

Thanks again for the clarifications.


You're welcome. It's unfortunate that we don't have a more informative message from GDAL that we can raise to users. "No such file or directory" is too opaque.

--
Sean Gillies


Re: returning subset of raster data via API

 

Per https://gdal.org/user/virtual_file_systems.html#drivers-supporting-virtual-file-systems and confirmed by Even R, GDAL's NetCDF driver cannot write to /vsi*. Perhaps a workaround can be to drop to a tempfile approach when writing NetCDF.

Thanks again for the clarifications.


Re: returning subset of raster data via API

Sean Gillies
 

Hi Tom,

On Mon, Aug 3, 2020 at 5:32 PM Tom Kralidis <tom.kralidis@...> wrote:

(sorry, just getting back to this!)

Thanks Sean. Indeed the MemoryFile strategy works as expected:

import sys

import rasterio
from rasterio.io import MemoryFile

with rasterio.open(sys.argv[1]) as src:
    with MemoryFile() as memfile:
        profile = src.profile
        with memfile.open(**profile) as dest:
            dest.write(src.read())

Are there any known issues with reading/subsetting, and writing back a NetCDF via MemoryFile? Traceback:

(pygeoapi) (base) tomkralidis@aruba:~/Dev/pygeoapi/pygeoapi-tomkralidis$ python dd.py tests/data/CMIP5_rcp8.5_annual_abs_latlon1x1_PCP_pctl25_P1Y.nc 
Traceback (most recent call last):
  File "rasterio/_io.pyx", line 1135, in rasterio._io.DatasetWriterBase.__init__
  File "rasterio/_err.pyx", line 205, in rasterio._err.exc_wrap_pointer
rasterio._err.CPLE_OpenFailedError: Unable to create netCDF file /vsimem/fa663a74-d548-494d-a2ab-880bdff831c4/fa663a74-d548-494d-a2ab-880bdff831c4. (Error code 2): No such file or directory .

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "dd.py", line 9, in <module>
    with memfile.open(**profile) as dest:
  File "/Users/tomkralidis/opt/miniconda3/lib/python3.7/site-packages/rasterio/env.py", line 382, in wrapper
    return f(*args, **kwds)
  File "/Users/tomkralidis/opt/miniconda3/lib/python3.7/site-packages/rasterio/io.py", line 136, in open
    nodata=nodata, sharing=sharing, **kwargs)
  File "rasterio/_io.pyx", line 1139, in rasterio._io.DatasetWriterBase.__init__
rasterio.errors.RasterioIOError: Unable to create netCDF file /vsimem/fa663a74-d548-494d-a2ab-880bdff831c4/fa663a74-d548-494d-a2ab-880bdff831c4. (Error code 2): No such file or directory .

Versions:

>>> import rasterio
>>> rasterio.__version__
'1.1.5'
>>> rasterio.gdal_version()
'3.0.4'

Thanks

..Tom

I don't use netcdf files with GDAL often and have never tried to write a netcdf file to vsimem. It could be that there is a combination of creation or config options required that I don't know about. I see some mention about platform requirements in https://github.com/OSGeo/gdal/pull/786 but am not sure what to make of that.

--
Sean Gillies


Re: returning subset of raster data via API

 

(sorry, just getting back to this!)

Thanks Sean. Indeed the MemoryFile strategy works as expected:

import sys

import rasterio
from rasterio.io import MemoryFile

with rasterio.open(sys.argv[1]) as src:
    with MemoryFile() as memfile:
        profile = src.profile
        with memfile.open(**profile) as dest:
            dest.write(src.read())

Are there any known issues with reading/subsetting, and writing back a NetCDF via MemoryFile? Traceback:

(pygeoapi) (base) tomkralidis@aruba:~/Dev/pygeoapi/pygeoapi-tomkralidis$ python dd.py tests/data/CMIP5_rcp8.5_annual_abs_latlon1x1_PCP_pctl25_P1Y.nc 
Traceback (most recent call last):
  File "rasterio/_io.pyx", line 1135, in rasterio._io.DatasetWriterBase.__init__
  File "rasterio/_err.pyx", line 205, in rasterio._err.exc_wrap_pointer
rasterio._err.CPLE_OpenFailedError: Unable to create netCDF file /vsimem/fa663a74-d548-494d-a2ab-880bdff831c4/fa663a74-d548-494d-a2ab-880bdff831c4. (Error code 2): No such file or directory .

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "dd.py", line 9, in <module>
    with memfile.open(**profile) as dest:
  File "/Users/tomkralidis/opt/miniconda3/lib/python3.7/site-packages/rasterio/env.py", line 382, in wrapper
    return f(*args, **kwds)
  File "/Users/tomkralidis/opt/miniconda3/lib/python3.7/site-packages/rasterio/io.py", line 136, in open
    nodata=nodata, sharing=sharing, **kwargs)
  File "rasterio/_io.pyx", line 1139, in rasterio._io.DatasetWriterBase.__init__
rasterio.errors.RasterioIOError: Unable to create netCDF file /vsimem/fa663a74-d548-494d-a2ab-880bdff831c4/fa663a74-d548-494d-a2ab-880bdff831c4. (Error code 2): No such file or directory .

Versions:

>>> import rasterio
>>> rasterio.__version__
'1.1.5'
>>> rasterio.gdal_version()
'3.0.4'

Thanks

..Tom


Re: How to extract DATA_ENCODING creation option type from an input grib2 file for use in the output grib2 file

 

Thanks Shane/Sean. On the pygeoapi side, I also stumbled on this issue and was able to allow for user specific options in configuration and passing that on to rasterio accordingly. Definitely not dynamic, but worked for me, at least from the perspective of consistent (enough) production of GRIB2.

Cheers

..Tom


window_transform doesn't produce accurate coordinate transforms

pranay1117@...
 

Hi, I'm new to rasterio and working with geographical data. So far, it's been a great experience working with this library - but I've run into an issue when working with the Window module. My use case is as follows: 
- I obtain 2 bounding coordinates from the user (top-left and bottom-right) and convert them to (row, col) using `source.index(lon, lat)`.
- I create a window using these bounds and read the dataset from that. Also, I get the appropriate transform for this window.
selected_window = Window.from_slices((row0, row1), (col0, col1))
data = src.read(1, window=selected_window)
affine_transform = src.window_transform(selected_window)

Afterwards, I am using this `data` numpy array to select various points in these bounds and plot them on a map. However, the (lat, long) coordinates obtained using the affine_transform above are not very accurate. When I run the following:
print("Lat/Lon bounds: {}".format(src.window_bounds(selected_window)))
print("Transform at (0,0): {}".format(affine_transform * (0,0)))
r, c = data.shape
print("Transform at {}: {}".format((r-1, c-1), affine_transform * (r-1, c-1)))

The output is as follows:
Source lat/lon bounds: (72.45 20.69, 82.39, 28.78)
Transform at (0,0): (72.42, 28.78)
Transform at (970, 1193): (80.53, 18.841)

As is evident from the output, the latitude must lie between (20.69,28.78) but the transform value 18.841 at (970, 1193) is not in that interval. 

I saw some earlier posts in this group that suggested updating to version 1.1.5. I've done that but the issue persists. 
Any help will be greatly appreciated. Thanks!



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?
 

121 - 140 of 698