Date   

Re: gdal_proximity

Anderson Roberto da Silva
 

Hello,
First sorry for my bad English.
Maybe I can help you and other people with the same doubt. Currently I' ve been working on a project that I need this functionality.

def proximity(in_raster, out_raster, list_values=[]):
"""Calculate proximity
:param in_raster: File path from input raster
:param out_raster: File path from output raster
:param list_values: List of the values.
If not specified, it's considered 0 value
:return: Return the file path from output raster proximity
"""
if len(list_values) == 0:
values = 'VALUES=' + str(0)
else:
values = 'VALUES=' + ','.join(str(value) for value in list_values)

in_raster_ds = gdal.Open(in_raster, gdal.GA_ReadOnly)
x_size = in_raster_ds.RasterXSize
y_size = in_raster_ds.RasterYSize

raster_band1 = in_raster_ds.GetRasterBand(1)

driver = gdal.GetDriverByName('GTiff')

dest_ds = driver.Create(out_raster, x_size, y_size, 1, gdal.GDT_Float32,
creation_options)

dest_ds.SetProjection(in_raster_ds.GetProjection())
dest_ds.SetGeoTransform(in_raster_ds.GetGeoTransform())

dest_ds_band1 = dest_ds.GetRasterBand(1)

gdal.ComputeProximity(raster_band1, dest_ds_band1, [values,
"DISTUNITS=GEO"])

in_raster_ds, raster_band1, dest_ds = None, None, None

return out_raster

For test:
class TestProximity(unittest.TestCase):
in_raster = path.join(RASTER_FOLDER, 'rodovias_raster.tif')
out_raster = path.join(RASTER_FOLDER, 'rodovias_proximity.tif')

def test_proximity(self):
proximity(self.in_raster, self.out_raster)

def test_proximity_values(self):
proximity(self.in_raster, self.out_raster, [1])

But first, you need to rasterize your roads.
My function to this:
def rasterize_fit_by_layer(geopackage, target_layer, field, grid_size,
out_raster, extend_layer=None):
"""Rasterize vector layer to image only in the layer extents

:param geopackage: Geopackage file
:param target_layer: Name of the layer on geopackage to rasterize
:param field: Name of the field that contains the integer values
:param grid_size: Cell size of the grid raster
:param extend_layer: Name of the layer on geopackage to extents (Optional)
:param out_raster: File path to save output raster
:return: File path from output raster
"""
source = ogr.Open(geopackage)

target_layer_rasterize = source.GetLayer(target_layer)
if target_layer_rasterize is None:
return None

layer_def = target_layer_rasterize.GetLayerDefn()
if layer_def.GetFieldIndex(field) < 0:
return None

spatial_ref = target_layer_rasterize.GetSpatialRef()

if extend_layer is None:
extents = target_layer_rasterize.GetExtent()
else:
layer_extent = source.GetLayer(extend_layer)
if layer_extent is None:
return None
extents = layer_extent.GetExtent()

x_min, x_max, y_min, y_max = extents

cols = int((x_max - x_min) / grid_size)
rows = int((y_max - y_min) / grid_size)

x_res = (x_max - x_min) / float(cols)
y_res = (y_max - y_min) / float(rows)

driver = gdal.GetDriverByName('GTiff')
dest_ds = driver.Create(out_raster, cols, rows, 1, gdal.GDT_Byte,
creation_options)
dest_ds.SetGeoTransform((x_min, x_res, 0, y_max, 0, -y_res))
dest_ds.SetProjection((spatial_ref.ExportToWkt()))

options = ['ATTRIBUTE=' + field]
# options = ['ATTRIBUTE=' + field, 'ALL_TOUCHED=TRUE']
gdal.RasterizeLayer(dest_ds, [1], target_layer_rasterize, None,
options=options)

dest_ds, source = None, None
layer_extent, target_layer_rasterize = None, None
Change for your uses.

Hope this helps.

-- 
Anderson Roberto da Silva
Engenheiro Cartógrafo





On Mon, Jan 4, 2021 at 5:21 PM Spencer Gardner <spencergardner@...> wrote:
Hi all. I'm using Rasterio for a raster analysis but I need to calculate proximity as is done in gdal_proximity.

Is this implemented somewhere in Rasterio? I don't see anything like it in the docs. Or if not, is there a way to run the GDAL operation on the array? Or even better, I could run a numpy operation on it, but I'm not finding anything in numpy that is suitable.

As an example, I am rasterizing a line feature that represents a segment of road. My analysis needs to account for proximity to the road such that cells touched by the road or immediately adjacent are weighted higher than roads that are 500 ft away. My intent with the proximity analysis was to use results to generate these weights.

Any pointers are much appreciated. Thanks!

Spencer


gdal_proximity

Spencer Gardner
 

Hi all. I'm using Rasterio for a raster analysis but I need to calculate proximity as is done in gdal_proximity.

Is this implemented somewhere in Rasterio? I don't see anything like it in the docs. Or if not, is there a way to run the GDAL operation on the array? Or even better, I could run a numpy operation on it, but I'm not finding anything in numpy that is suitable.

As an example, I am rasterizing a line feature that represents a segment of road. My analysis needs to account for proximity to the road such that cells touched by the road or immediately adjacent are weighted higher than roads that are 500 ft away. My intent with the proximity analysis was to use results to generate these weights.

Any pointers are much appreciated. Thanks!

Spencer


Re: Any idea why I'd be getting a missing region when reprojecting?

hi@...
 

Another update, for posterity. I tried doing the same reprojection using the gdal.Warp and it worked perfectly with far less effort around setting up transforms and whatnot. I filed an issue with rasterio: github.com/mapbox/rasterio/issues/2074


Re: Any idea why I'd be getting a missing region when reprojecting?

hi@...
 

Just an update: I set the env variable CHECK_WITH_INVERT_PROJ to YES and was able to generate the output using the vertical slice but again, when I try to do this with the global raster I get a blank output. I feel like I may have solved the original problem I posted about but am now running into an issue with generating the global raster.


Any idea why I'd be getting a missing region when reprojecting?

hi@...
 

I'll try my best to describe the issue although it's really confusing me so I'm not sure how well I can explain it.

I'm trying to reproject some satellite imagery in the form of a raster that covers the entire world using the CRS of the imagery to EPSG:4326 using the reproject function. The raster before being reprojected uses a sinusoidal projection with meters as units and looks like this:


When I try to reproject it, I get the following:


As you can see, anything to the right of a certain horizontal value in the original raster is just missing entirely.

I really can't figure out why this is happening and I've tried changing dozens of things in my code to fix it with very little success. I've provided some of my code below.

Sorry for the code formatting, couldn't see any option for better formatting.


for label in layers.keys():
    with rasterio.open(os.path.join(OUTPUT_DIR, f'{label}_intermediate.tif')) as intermediate_file:
        # Compute transform needed for output
        final_transform, final_width, final_height = calculate_default_transform(
            intermediate_file.crs,
            FINAL_CRS,
            intermediate_file.width,
            intermediate_file.height,
            *intermediate_file.bounds,
            resolution = 0.05
        )
 
        # Extract parameters from intermediate file
        kwargs = intermediate_file.meta.copy()
 
        # Update using the final parameters
        kwargs.update({
            'crs': FINAL_CRS,
            'transform': final_transform,
            'width': final_width,
            'height': final_height
        })
 
        # Compute reprojection and write to file
        with rasterio.open(os.path.join(OUTPUT_DIR, f'{label}_final.tif'), 'w', **kwargs) as final_file:
            dst, tsf = reproject(
                source = rasterio.band(intermediate_file, 1),
                destination = rasterio.band(final_file, 1),
                src_transform = intermediate_file.transform,
                src_crs = intermediate_file.crs,
                dst_transform = final_transform,
                dst_crs = FINAL_CRS,
                dst_resolution = 0.05,
                resampling = Resampling.average,
                warp_mem_limit = 1024 * REPROJECTION_MEMORY_GB,
                kwargs = {'COMPRESSION': 'LZW', 'TILED': 'YES'}
            )

So, like I've said, I feel like I've tried everything and with very little success. Here's a summary.
 
  • I tried changing the ordering of width and height since rasterio's functions aren't consistent in the ordering but this didn't help.
  • I tried changing the resampling method... nothing.
  • I thought it might be a memory problem so increased/decreased the memory of the reproject function... nothing. Also downscaled the original raster to reduce memory but again nothing.
  • I tried processing a vertical slice of the data covering only the affected region but the output raster was still hidden. However, when I removed the resolution parameters so it would maintain the original resolution the region did manage to appear properly. But, when I then try to repeat this on the global raster I just get a blank output which seems to happen when I process files that are too large although I'm not sure why.
  • Continuing with the vertical slice... if I set the resolution parameters back to how they are in the code but then modify the output's transform I'm able to have the regions appear again. I modify the transform by rounding the upper left corner coordinates so that the bound is from -180 to 180 rather than what is originally set which has some floating point error and is more like -179.999995 to 180.000005. However, again, when I run this on the global raster I just get an empty raster. I don't think it's a memory issue since it's using the same resolution used to generate the global raster with the missing region.
  • I was able to slightly increase the amount of data that was visible but I can't recall what exactly I changed to do this. I know it involved increasing/decreasing one of the width variables or something along those lines which slightly extended the "visible" region. Not very much, but enough that there was maybe another degree that was visible. I tried further increasing/decreasing this variable but that didn't help.
I'm definitely feeling fairly lost. There's a good chance I'm simply doing something that I don't understand well as I'm new to this stuff. Let me know if there's any additional info I can provide. Any help would be greatly appreciated :)


Re: Binary raster to geojson

Sean Gillies
 

https://rasterio.readthedocs.io/en/latest/cli.html#shapes is a command line program for extracting GeoJSON shapes from a raster. It's a thin wrapper around rasterio.features.dataset_features() (see https://github.com/mapbox/rasterio/blob/master/rasterio/rio/shapes.py).


On Sat, Dec 26, 2020 at 11:32 PM <ashnair0007@...> wrote:
Given a binary numpy array (array of 0s and 1s where 0=background), what would be the best way to polygonise this raster and write it to a geojson?
_._,_._

--
Sean Gillies


Re: Get PROJ.4 representation of CRS object

Sean Gillies
 

Hi,

Rasterio's CRS.to_dict and .to_prog4 *do* return PROJ4 mappings and strings but we reduce to +init=epsg:xxxx when allowed. In the new world where there is no longer an epsg file to init from, it makes more sense to return the expanded PROJ4 form, I agree. We could do this for 1.2.0 at the risk of breaking programs with code like

epsg_code = crs.to_dict()["init"].split(":")[-1]

Nobody should need to do that as we have a .to_epsg method, but one of the laws of library development is that there will be some users critically dependent on any public property of the library, no matter if it was exposed by accident or by intent :)


On Sat, Dec 26, 2020 at 10:29 AM Denis Rykov <rykovd@...> wrote:
Thanks! This is exactly what I was looking for.

BTW GDAL has an appropriate method to export CRS to proj4:

>>> srs.ExportToProj4()
'+proj=utm +zone=5 +datum=WGS84 +units=m +no_defs '

In my opinion it would be nice to have it in rasterio.CRS as well.

On Sun, Dec 27, 2020 at 12:19 AM Alan Snow <alansnow21@...> wrote:
pyproj.CRS has a to_proj4 method that should be able to do that: https://pyproj4.github.io/pyproj/stable/examples.html,_._,_

--
Sean Gillies


Binary raster to geojson

ashnair0007@...
 

Given a binary numpy array (array of 0s and 1s where 0=background), what would be the best way to polygonise this raster and write it to a geojson?


Re: Get PROJ.4 representation of CRS object

Denis Rykov
 

Thanks! This is exactly what I was looking for.

BTW GDAL has an appropriate method to export CRS to proj4:

>>> srs.ExportToProj4()
'+proj=utm +zone=5 +datum=WGS84 +units=m +no_defs '

In my opinion it would be nice to have it in rasterio.CRS as well.

On Sun, Dec 27, 2020 at 12:19 AM Alan Snow <alansnow21@...> wrote:
pyproj.CRS has a to_proj4 method that should be able to do that: https://pyproj4.github.io/pyproj/stable/examples.html


Re: Get PROJ.4 representation of CRS object

Alan Snow
 

pyproj.CRS has a to_proj4 method that should be able to do that: https://pyproj4.github.io/pyproj/stable/examples.html


Get PROJ.4 representation of CRS object

Denis Rykov
 

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?


Re: What is the raster merge criteria algorithm on rasterio.merge?

James McBride
 

Hi Eduardo,
I'd suggest taking a look at the method argument that merge accepts. There are a few pre-defined methods, and you can also define your own function for selecting which pixel to select from the datasets to be merged. I copied the documentation below:

    method : str or callable
        pre-defined method:
            first: reverse painting
            last: paint valid new on top of existing
            min: pixel-wise min of existing and new
            max: pixel-wise max of existing and new
        or custom callable with signature:
 
        def function(old_data, new_data, old_nodata, new_nodata, index=None, roff=None, coff=None):
 
            Parameters
            ----------
            old_data : array_like
                array to update with new_data
            new_data : array_like
                data to merge
                same shape as old_data
            old_nodata, new_data : array_like
                boolean masks where old/new data is nodata
                same shape as old_data
            index: int
                index of the current dataset within the merged dataset collection
            roff: int
                row offset in base array
            coff: int
                column offset in base array
 


What is the raster merge criteria algorithm on rasterio.merge?

Eduardo <eduardosteps@...>
 

I'm creating a mosaic from many sentinel images with a bunch of overlapping areas, although I'm trying to understand what is criteria. For example:

This is a small piece totally clear of clouds

1

And this is a larger piece with many clouds

2

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?


rasterio.windows.from_bounds height and width int parameters

Paolo Corti
 

Hi all

I am trying to figure out how to pass the height and width integer
parameters to the from_bounds method. In my case they look not to be
honored and the window returned has slightly different height and
width (which are float).

You can reproduce my issue with this code:

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)
Any idea why this is happening? Thanks a lot in advance
Paolo

--
Paolo Corti
Geospatial software developer
web: http://www.paolocorti.net
twitter: @capooti
skype: capooti
#drt3jc1


rio-color 1.0.1

Sean Gillies
 

Hi all,

Version 1.0.1 of rio-color is on PyPI now: https://pypi.org/project/rio-color/1.0.1/#files. Thank you, Kyle Barron, for showing how to get started with cibuildwheel.

Please note that there are no Python 3.9 wheels for the current stable version of rasterio. If you want to install rio-color for 3.9, you'll need to get rasterio==1.2b1.

Happy December Solstice, everyone.

--
Sean Gillies


Re: Rasterio 1.2b1

Alan Snow
 

I think this is the problem: https://pyproj4.github.io/pyproj/stable/gotchas.html#internal-proj-error-sqlite-error-on-select


Re: Rasterio 1.2b1

vincent.sarago@...
 

I just installed rasterio==1.2.0b1 on a fresh virtual env (with no GDAL/PROJ env set) and I'm getting a proj error 

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

Not sure what's going on :-( 
 


Re: rio-color 1.0.1.dev0 wheels

Sean Gillies
 

Hi all,

I got reports that the 1.0.1.dev0 wheels install and work, so I've made a more stable set and a 1.0.1a1 release. Target date for 1.0.1 is 2020-12-21. See https://github.com/mapbox/rio-color/issues/71 for details and further updates.

On Wed, Dec 16, 2020 at 3:14 PM Sean Gillies <sean@...> wrote:
Hi all,

Lack of an up to date wheel building system has been blocking us from making a 1.0.1 release. See https://github.com/mapbox/rio-color/issues/65. Kyle Barron suggested using cibuildwheel and I've succeeded in making a wide variety of wheels: https://pypi.org/manage/project/rio-color/release/1.0.1.dev0/. If you're willing to `pip install rio-color==1.0.1.dev0` and report success or failure, we'd appreciate it. Be aware of a few things:

1. That command will pull in dev versions of rio-color's dependencies, such as rasterio 1.2b1, so try this in a disposable Python environment.
2. There are some rio-color wheels that have no match in rasterio 1.2b1. So, if you use pip to install rio-color 1.0.1.dev0 on Windows, pip will find no rasterio wheel for Windows, and the installation will fail.


--
Sean Gillies


Re: rasterio.features.shapes with holes in polygons

Paolo Corti
 

On Tue, Dec 15, 2020 at 6:37 PM Sean Gillies via groups.io
<sean=mapbox.com@groups.io> wrote:

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?

Hi Sean!

Thanks a lot for your hint. My problem was caused by the fact that I
was merging the polygons returned by rasterio.features.shape using
unary_union. I changed the approach and I created a multipolygon from
the returned polygons using shapely.geometry.shape and this made the
trick.
However one thing which I have noticed is that the returned geometry
is invalid in some cases. So I had to buffer to avoid errors when writing to a
shapefile. Is this a good approach or you suggest something different?

MultiPolygon([shape(geom) for geom in polygons].buffer(0)

Thanks again for your help and for this wonderful project - and for
its siblings too :)


Paolo


rio-color 1.0.1.dev0 wheels

Sean Gillies
 

Hi all,

Lack of an up to date wheel building system has been blocking us from making a 1.0.1 release. See https://github.com/mapbox/rio-color/issues/65. Kyle Barron suggested using cibuildwheel and I've succeeded in making a wide variety of wheels: https://pypi.org/manage/project/rio-color/release/1.0.1.dev0/. If you're willing to `pip install rio-color==1.0.1.dev0` and report success or failure, we'd appreciate it. Be aware of a few things:

1. That command will pull in dev versions of rio-color's dependencies, such as rasterio 1.2b1, so try this in a disposable Python environment.
2. There are some rio-color wheels that have no match in rasterio 1.2b1. So, if you use pip to install rio-color 1.0.1.dev0 on Windows, pip will find no rasterio wheel for Windows, and the installation will fail.

--
Sean Gillies

221 - 240 of 905