Date   

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


Re: Rasterio 1.2b1

Alan Snow
 

Is that because on the other hand GDAL 3.2.0 and/or PROJ 7.2.0 are heavier than their counterparts from release 1.1.8?
The new wheel now includes the proj.db, which is close to the size of grids that were removed. Also, PROJ requires more dependencies, so that is also a likely factor.

would you expect there to be problems if I were to have pyproj<3.0.0 and rasterio>=1.2 installed in my Python environment, due to their different PROJ versions, or would that still work OK?
It might work, but I wouldn't recommend it. At the very least, there are different versions proj.db in each one. If you do go that route, I would recommend you proceed with caution.


Re: Rasterio 1.2b1

Guillaume Lostis
 

Hi Sean,

Thanks a lot for this upcoming release, I am really looking forward to using some of its new features!

Out of interest, the 1.2b1 wheels do not seem to be lighter than the 1.1.8 wheels even though they do not contain PROJ data. Is that because on the other hand GDAL 3.2.0 and/or PROJ 7.2.0 are heavier than their counterparts from release 1.1.8?

Also Sean and Alan, would you expect there to be problems if I were to have pyproj<3.0.0 and rasterio>=1.2 installed in my Python environment, due to their different PROJ versions, or would that still work OK?

Best,

Guillaume

On Tue, Dec 15, 2020 at 6:59 PM Alan Snow <alansnow21@...> wrote:
Thanks for working on this Sean!

I verified that PROJ_NETWORK=ON works with a basic test:

rio_geom.py:
from rasterio.warp import transform_geom

geometry = [
{
"type": "Polygon",
"coordinates": [
[
[-94.07955380199459, 41.69085871273774],
[-94.06082436942204, 41.69103313774798],
[-94.06063203899649, 41.67932439500822],
[-94.07935807746362, 41.679150041277325],
[-94.07955380199459, 41.69085871273774],
]
],
}
]


print(transform_geom("+proj=latlon", "+proj=utm +zone=10 +datum=NAD27", geometry))

$ python rio_geom.py
[{'type': 'Polygon', 'coordinates': [[(2914912.501340564, 5039833.18814224), (2916474.3810722474, 5040429.817378172), (2916971.9920050735, 5039127.068439782), (2915410.0220987815, 5038530.542053506), (2914912.501340564, 5039833.18814224)]]}]
$ PROJ_NETWORK=ON python rio_geom.py
[{'type': 'Polygon', 'coordinates': [[(2914912.2992929644, 5039834.208342302), (2916474.1765534407, 5040430.837641376), (2916971.7831669645, 5039128.087552843), (2915409.815547156, 5038531.561302773), (2914912.2992929644, 5039834.208342302)]]}]


Re: rasterio.features.shapes with holes in polygons

Sean Gillies
 

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:


It looks like the holes aren't lost there. How did you make the second image?

On Tue, Dec 15, 2020 at 3:47 PM Paolo Corti <pcorti@...> wrote:
Hello users and devs

First: thanks a lot for this wonderful project, I am really enjoying using it for my geospatial needs.

For a specific process I am using rasterio.features.shapes and it is working greatly for my purpose. However I have a few cases where I am experiencing strange results.
For example, see this 1 band raster: 

Screen Shot 2020-12-15 at 5.33.13 PM.png
I would like to extract the polygons for the yellow area (value: 3).

When running the following code:

with rasterio.open(cover_path) as src:
cover_arr = src.read(1)
mask = (cover_arr == 3)
polygons = shapes(cover_arr, mask=mask, transform=src.transform)

I get these results (in blue):

Screen Shot 2020-12-15 at 5.35.37 PM.png

Apparently holes are not correctly parsed.

It would be great to know if I am doing something wrong here or if it is something which is not working properly.

For your reference, I am attaching a copy of the raster I am using here so you can reproduce this (sample_cover.tif)

Thanks in advance!
Paolo

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



--
Sean Gillies


rasterio.features.shapes with holes in polygons

Paolo Corti
 

Hello users and devs

First: thanks a lot for this wonderful project, I am really enjoying using it for my geospatial needs.

For a specific process I am using rasterio.features.shapes and it is working greatly for my purpose. However I have a few cases where I am experiencing strange results.
For example, see this 1 band raster: 

Screen Shot 2020-12-15 at 5.33.13 PM.png
I would like to extract the polygons for the yellow area (value: 3).

When running the following code:

with rasterio.open(cover_path) as src:
cover_arr = src.read(1)
mask = (cover_arr == 3)
polygons = shapes(cover_arr, mask=mask, transform=src.transform)

I get these results (in blue):

Screen Shot 2020-12-15 at 5.35.37 PM.png

Apparently holes are not correctly parsed.

It would be great to know if I am doing something wrong here or if it is something which is not working properly.

For your reference, I am attaching a copy of the raster I am using here so you can reproduce this (sample_cover.tif)

Thanks in advance!
Paolo

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


Re: Rasterio 1.2b1

Alan Snow
 

Thanks for working on this Sean!

I verified that PROJ_NETWORK=ON works with a basic test:

rio_geom.py:
from rasterio.warp import transform_geom

geometry = [
{
"type": "Polygon",
"coordinates": [
[
[-94.07955380199459, 41.69085871273774],
[-94.06082436942204, 41.69103313774798],
[-94.06063203899649, 41.67932439500822],
[-94.07935807746362, 41.679150041277325],
[-94.07955380199459, 41.69085871273774],
]
],
}
]


print(transform_geom("+proj=latlon", "+proj=utm +zone=10 +datum=NAD27", geometry))

$ python rio_geom.py
[{'type': 'Polygon', 'coordinates': [[(2914912.501340564, 5039833.18814224), (2916474.3810722474, 5040429.817378172), (2916971.9920050735, 5039127.068439782), (2915410.0220987815, 5038530.542053506), (2914912.501340564, 5039833.18814224)]]}]
$ PROJ_NETWORK=ON python rio_geom.py
[{'type': 'Polygon', 'coordinates': [[(2914912.2992929644, 5039834.208342302), (2916474.1765534407, 5040430.837641376), (2916971.7831669645, 5039128.087552843), (2915409.815547156, 5038531.561302773), (2914912.2992929644, 5039834.208342302)]]}]


Rasterio 1.2b1

Sean Gillies
 

Hi all,

Wheels for Python versions 3.6-3.8 and manylinux1 and macos (built with xcode 9.3) are on PyPI this morning. I'd appreciate it if you could give it a try: python -m pip install rasterio==1.2b1.

At last, these wheels include GDAL 3.2.0 and PROJ 7.2.0. Please note that these wheels include no PROJ data. You must provide your own data or set PROJ_NETWORK=ON in your environment to use the PROJ CDN. See https://proj.org/resource_files.html for details.

I'm particularly interested in feedback about including leaving out the PROJ data. If it's going to impact the way you deploy rasterio, this is the time to let us know.

Yours,

--
Sean Gillies


resampling raster after clipping/masking

Eyal Saiet
 

Hello,
I have found references in rasterio.readthedocs.io both on how to crop/clip a raster and how to resample. I generally understand the steps and sought to combine the two steps into one:

#Clipping Raster
with rasterio.open(dem_raster_pd['rasters'][0]) as raster:
    #create a mask
    out_image,out_transform=rio.mask.mask(raster,shapes,crop=True)
    #adding meta of source
    out_meta=raster.meta


    #updading meta
    out_meta.update({"driver":"GTiff","height":out_image.shape[1],
                                     "width":out_image.shape[2],
                                    "transform":out_transform})

#Because I am resampling next I skipped writing crp_raster
   ##Rescaling to 1 m
    ras_xy=1 #(m)

   # figuring out the factor
    rs_factor=((ras_xy/out_image.shape[1],ras_xy/out_image.shape[2])
    rs_in_factor=(1/rs_factor[0],1/rs_factor[1])

    #But because I don't have a crp_raster raster object I cannot use read and apply yhr out_shape
     data_resample=crp_raster.read(out_shape=(crp_raster.count,int(crp_raster.height*rs_in_factor[0]),int(crp_raster.width*rs_in_factor[1])),resampling=Resampling.bilinear,masked=True)#,masked=True


    #Scaling image trnasform
    transform_resmp=out_transform * out_transform.scale(
       (out_image.shape[2]/data_resample.shape[-1]),
        (out_image.shape[1]/data_resample.shape[-2])
        )

    #updating meta
    out_meta_res=crp_raster.meta
    out_meta_res.update({"height":data_resample.shape[-2],"width":data_resample.shape[-1],"transform":transform_resmp})


    #New filename
    new_ras_nm=products_dir+os.sep+raster.name.split('.')[0]+"_crp_rsmp_1m.tif"

    #Saving
with rasterio.open(new_ras_nm,"w",**out_meta_res) as dest_crp_resamp:
    print('Resampled')
    dest_crp_resamp.write(data_resample)     

I am sure I am stuck over something obvious that currently I cannot see. If I could create a raster object,  I tried crp_raster=rasterio.open("w"..), then I could use the raster_objec.read(out_shape=).



21 - 40 of 698