1. main@rasterio.groups.io
   Date   
rasterio 1.2b4

Sean Gillies
 

Hi all,

I think this may be it, the last pre-1.2.0 release. Would love to know if you find any issues with it.

What a day!

--
Sean Gillies
rio-color 1.0.2

Sean Gillies
 

Hi all,

Version 1.0.2 of rio-color is on PyPI today. It has one bug fix, a change to the rasterio requirement's version specifier which filters out pre-releases. Installing rio-color 1.0.1 could entail one of the rasterio pre-releases. Thank you, Vincent Sarago and Kyle Barron, for pointing out the problem.

--
Sean Gillies
Re: reading with a window doesn't honour window shape

Loïc Dutrieux
 

Note that setting out_shape to (1200, 1200) will likely read a (1200, 1199) shaped chunk of data and resample it to (1200, 1200); which is not exactly the same as reading a 1200*1200 window.

Kind regards,
Loïc
________________________________________
From: main@rasterio.groups.io <main@rasterio.groups.io> on behalf of Paolo Corti <pcorti@gmail.com>
Sent: 15 January 2021 20:43:21
To: main@rasterio.groups.io
Subject: Re: [rasterio] reading with a window doesn't honour window shape

Thanks Guillaume, that did the trick :)
Would also love to know the reason why that is needed.
best
Paolo

On Fri, Jan 15, 2021 at 1:53 PM Guillaume Lostis <guillaume@lostis.com<mailto:guillaume@lostis.com>> wrote:
Hi,

If you add `out_shape=(1200, 1200)` to your `read()` call it returns an array with the right shape. I don't know exactly why you need to specify this though.

Best,

Guillaume

Le ven. 15 janv. 2021 à 17:48, Paolo Corti <pcorti@gmail.com<mailto:pcorti@gmail.com>> a écrit :
Hi all

I have the following problem: I want to read a dataset using a window with a specific shape, but the window shape isn't always honored from the returned array.

You can replicate this problem using this dataset: https://www.worldpop.org/geodata/summary?id=6334<;https://urldefense.com/v3/__https://www.worldpop.org/geodata/summary?id=6334__;!!DOxrgLBm!R8DE-M0LoCdGJkxHEE7d81JPP8UIGkAXjLZ_GJL5tBBF3QRGWwWHOBaSpab9fGiBA4Bpvcs$>

with rasterio.open(raster_pop_path) as src_pop:
window = Window(col_off=1470.4999917338748, row_off=2902.4961942222235, width=1200, height=1200)
pop_arr = src_pop.read(1, window=window)
print(pop_arr.shape)

print returns (1200, 1199) while it should be (1200, 1200)

Any idea why this could happen?
Thanks in advance

Paolo

--
Paolo Corti
Geospatial software developer
web: http://www.paolocorti.net<;https://urldefense.com/v3/__http://www.paolocorti.net__;!!DOxrgLBm!R8DE-M0LoCdGJkxHEE7d81JPP8UIGkAXjLZ_GJL5tBBF3QRGWwWHOBaSpab9fGiBh8FrJuw$>
twitter: @capooti
skype: capooti
#drt3jc1


--
Paolo Corti
Geospatial software developer
web: http://www.paolocorti.net<;https://urldefense.com/v3/__http://www.paolocorti.net__;!!DOxrgLBm!R8DE-M0LoCdGJkxHEE7d81JPP8UIGkAXjLZ_GJL5tBBF3QRGWwWHOBaSpab9fGiBh8FrJuw$>
twitter: @capooti
skype: capooti
#drt3jc1
Re: reading with a window doesn't honour window shape

Paolo Corti
 

Thanks Guillaume, that did the trick :)
Would also love to know the reason why that is needed.
best
Paolo

On Fri, Jan 15, 2021 at 1:53 PM Guillaume Lostis <guillaume@...> wrote:
Hi,

If you add `out_shape=(1200, 1200)` to your `read()` call it returns an array with the right shape. I don't know exactly why you need to specify this though.

Best,

Guillaume

Le ven. 15 janv. 2021 à 17:48, Paolo Corti <pcorti@...> a écrit :
Hi all

I have the following problem: I want to read a dataset using a window with a specific shape, but the window shape isn't always honored from the returned array.

You can replicate this problem using this dataset: https://www.worldpop.org/geodata/summary?id=6334

with rasterio.open(raster_pop_path) as src_pop:
window = Window(col_off=1470.4999917338748, row_off=2902.4961942222235, width=1200, height=1200)
pop_arr = src_pop.read(1, window=window)
print(pop_arr.shape)

print returns (1200, 1199) while it should be (1200, 1200)

Any idea why this could happen?
Thanks in advance

Paolo

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



--
Paolo Corti
Geospatial software developer
web: http://www.paolocorti.net
twitter: @capooti
skype: capooti
#drt3jc1
Re: reading with a window doesn't honour window shape

Guillaume Lostis
 

Hi,

If you add `out_shape=(1200, 1200)` to your `read()` call it returns an array with the right shape. I don't know exactly why you need to specify this though.

Best,

Guillaume

Le ven. 15 janv. 2021 à 17:48, Paolo Corti <pcorti@...> a écrit :
Hi all

I have the following problem: I want to read a dataset using a window with a specific shape, but the window shape isn't always honored from the returned array.

You can replicate this problem using this dataset: https://www.worldpop.org/geodata/summary?id=6334

with rasterio.open(raster_pop_path) as src_pop:
window = Window(col_off=1470.4999917338748, row_off=2902.4961942222235, width=1200, height=1200)
pop_arr = src_pop.read(1, window=window)
print(pop_arr.shape)

print returns (1200, 1199) while it should be (1200, 1200)

Any idea why this could happen?
Thanks in advance

Paolo

--
Paolo Corti
Geospatial software developer
web: http://www.paolocorti.net
twitter: @capooti
skype: capooti
#drt3jc1
reading with a window doesn't honour window shape

Paolo Corti
 

Hi all

I have the following problem: I want to read a dataset using a window with a specific shape, but the window shape isn't always honored from the returned array.

You can replicate this problem using this dataset: https://www.worldpop.org/geodata/summary?id=6334

with rasterio.open(raster_pop_path) as src_pop:
window = Window(col_off=1470.4999917338748, row_off=2902.4961942222235, width=1200, height=1200)
pop_arr = src_pop.read(1, window=window)
print(pop_arr.shape)

print returns (1200, 1199) while it should be (1200, 1200)

Any idea why this could happen?
Thanks in advance

Paolo

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

Spencer Gardner
 

Thanks. I'll give that a go.


On Wed, Jan 13, 2021 at 7:45 AM Sean Gillies <sean.gillies@...> wrote:
Hi Spencer,

On Mon, Jan 4, 2021 at 1: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

You're not overlooking anything. The rasterio project doesn't surface GDALComputeProximity. There are pointers to equivalent functionality in scipy at https://github.com/mapbox/rasterio/issues/531.

Hope this helps,

--
Sean Gillies
Re: Rasterio installation on Amazon/AWS/EC2/notebook instance

Sean Gillies
 

Hi,

On Wed, Jan 13, 2021 at 7:33 AM <robmarkcole@...> wrote:
I am attempting to install rasterio on an AWS notebook instance (essentially managed EC2 with Jupyter installed).
This should be straightforward but I have already sunk several hours without success. 

The instance has a bunch of conda environments available so I wish to install into one named python3.
The following results:

sh-4.2$ conda install -c conda-forge rasterio --name python3
Collecting package metadata (current_repodata.json): done
Solving environment: failed with initial frozen solve. Retrying with flexible solve.
Solving environment: failed with repodata from current_repodata.json, will retry with next repodata source.
Collecting package metadata (repodata.json): done
Solving environment: - 
The environment is inconsistent, please check the package plan carefully
The following packages are causing the inconsistency:
 
  - defaults/noarch::jupyterlab==1.2.6=pyhf63ae98_0
  - defaults/linux-64::python-language-server==0.31.7=py36_0
  - defaults/linux-64::nb_conda==2.2.1=py36_0
  - defaults/noarch::numpydoc==0.9.2=py_0
  - defaults/noarch::anaconda-project==0.8.4=py_0
  - defaults/noarch::boto3==1.9.162=py_0
  - defaults/linux-64::s3transfer==0.2.1=py36_0
  - defaults/linux-64::nbconvert==5.6.1=py36_0
  - defaults/linux-64::bokeh==1.4.0=py36_0
  - defaults/noarch::jupyterlab_server==1.0.6=py_0
  - defaults/noarch::botocore==1.12.189=py_0
  - defaults/linux-64::jupyter==1.0.0=py36_7
  - defaults/linux-64::scikit-image==0.16.2=py36h0573a6f_0
  - defaults/linux-64::imageio==2.6.1=py36_0
  - defaults/linux-64::nb_conda_kernels==2.2.4=py36_0
  - defaults/linux-64::spyder==4.0.1=py36_0
  - defaults/linux-64::requests==2.22.0=py36_1
  - defaults/noarch::dask==2.11.0=py_0
  - defaults/noarch::ipywidgets==7.5.1=py_0
  - defaults/linux-64::widgetsnbextension==3.5.1=py36_0
  - defaults/noarch::s3fs==0.4.2=py_0
  - defaults/linux-64::notebook==6.0.3=py36_0
  - defaults/linux-64::anaconda-client==1.7.2=py36_0

At this point conda hangs. I also tried running:

$ conda update all

And repeating the process, but with the same result.
Any guidance greatly appreciated

I don't know anything about AWS notebook instances and am not a regular conda user, so my help might be a little sketchy. In my experience, one cannot incrementally upgrade or install new packages into a conda environment without breaking it. Instead, create a completely new environment, and expect that to be resolved properly.

--
Sean Gillies
Re: gdal_proximity

Sean Gillies
 

Hi Spencer,

On Mon, Jan 4, 2021 at 1: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

You're not overlooking anything. The rasterio project doesn't surface GDALComputeProximity. There are pointers to equivalent functionality in scipy at https://github.com/mapbox/rasterio/issues/531.

Hope this helps,

--
Sean Gillies
Rasterio installation on Amazon/AWS/EC2/notebook instance

robmarkcole@...
 

I am attempting to install rasterio on an AWS notebook instance (essentially managed EC2 with Jupyter installed).
This should be straightforward but I have already sunk several hours without success. 

The instance has a bunch of conda environments available so I wish to install into one named python3.
The following results:

sh-4.2$ conda install -c conda-forge rasterio --name python3
Collecting package metadata (current_repodata.json): done
Solving environment: failed with initial frozen solve. Retrying with flexible solve.
Solving environment: failed with repodata from current_repodata.json, will retry with next repodata source.
Collecting package metadata (repodata.json): done
Solving environment: - 
The environment is inconsistent, please check the package plan carefully
The following packages are causing the inconsistency:
 
  - defaults/noarch::jupyterlab==1.2.6=pyhf63ae98_0
  - defaults/linux-64::python-language-server==0.31.7=py36_0
  - defaults/linux-64::nb_conda==2.2.1=py36_0
  - defaults/noarch::numpydoc==0.9.2=py_0
  - defaults/noarch::anaconda-project==0.8.4=py_0
  - defaults/noarch::boto3==1.9.162=py_0
  - defaults/linux-64::s3transfer==0.2.1=py36_0
  - defaults/linux-64::nbconvert==5.6.1=py36_0
  - defaults/linux-64::bokeh==1.4.0=py36_0
  - defaults/noarch::jupyterlab_server==1.0.6=py_0
  - defaults/noarch::botocore==1.12.189=py_0
  - defaults/linux-64::jupyter==1.0.0=py36_7
  - defaults/linux-64::scikit-image==0.16.2=py36h0573a6f_0
  - defaults/linux-64::imageio==2.6.1=py36_0
  - defaults/linux-64::nb_conda_kernels==2.2.4=py36_0
  - defaults/linux-64::spyder==4.0.1=py36_0
  - defaults/linux-64::requests==2.22.0=py36_1
  - defaults/noarch::dask==2.11.0=py_0
  - defaults/noarch::ipywidgets==7.5.1=py_0
  - defaults/linux-64::widgetsnbextension==3.5.1=py36_0
  - defaults/noarch::s3fs==0.4.2=py_0
  - defaults/linux-64::notebook==6.0.3=py36_0
  - defaults/linux-64::anaconda-client==1.7.2=py36_0

At this point conda hangs. I also tried running:

$ conda update all

And repeating the process, but with the same result.
Any guidance greatly appreciated
Rasterio 1.2b3

Sean Gillies
 

Hi all,

Rasterio 1.2b3 is on PyPI and has a fix for https://github.com/mapbox/rasterio/issues/2079.

--
Sean Gillies
Georeference and warp a drone image based on sensor orientation and location using transform.from_gcps()

pgrayobx@...
 

I have been flying some drone surveys over the ocean and need to properly project and georeference and warp my images. I have all the information I think I need: lat, lon, altitude, yaw, pitch, and roll along with sensor optics params. I have to imagine this can be done with some existing package and I had thought some rasterio functions were suitable, but I can't find any. So I've been using the cameratransform package to get the GPS positions of the corners of my image to feed into rasterio:
 
import cameratransform as ct
 
# camera parameters
 
cam = ct.Camera(ct.RectilinearProjection(focallength_mm=f,
                                         sensor=sensor_size,
                                         image=image_size),
               ct.SpatialOrientation(elevation_m=img.altitude,
                                     tilt_deg=pitch+sensor_offset,
                                     roll_deg=roll,
                                    heading_deg=yaw))
 
# gps pts are lat lon
cam.setGPSpos(img.latitude, img.longitude, img.altitude)
 
# these are the coordinates of the image corners
coords = np.array([cam.gpsFromImage([0 , 0]), \
    cam.gpsFromImage([image_size[0]-1 , 0]), \
    cam.gpsFromImage([image_size[0]-1, image_size[1]-1]), \
    cam.gpsFromImage([0 , image_size[1]-1])])

From there I'm not certain what to do. I've tried basically saying these corners are GCPs and trying to warp the image in rasterio:
 
import rasterio
 
gcp1 = rasterio.control.GroundControlPoint(row=0, col=0, x=coords[0,1], y=coords[0,0], z=coords[0,2], id=None, info=None)
gcp2 = rasterio.control.GroundControlPoint(row=image_size[0]-1, col=0, x=coords[1,1], y=coords[1,0], z=coords[1,2], id=None, info=None)
gcp3 = rasterio.control.GroundControlPoint(row=image_size[0]-1, col=image_size[1]-1, x=coords[2,1], y=coords[2,0], z=coords[2,2], id=None, info=None)
gcp4 = rasterio.control.GroundControlPoint(row=0, col=image_size[1]-1, x=coords[3,1], y=coords[3,0], z=coords[3,2], id=None, info=None)
 
# Register GDAL format drivers and configuration options with a
# context manager.
with rasterio.Env():
    
    # open the original image to get some of the basic metadata
    with rasterio.open(path_name, 'r') as src:
        profile = src.profile
                    
    # create rasterio transform
    tsfm = rasterio.transform.from_gcps([gcp1,gcp2,gcp3,gcp4])
    
    # I also tried this function but to no avail
    #tsfm = rasterio.warp.calculate_default_transform(rasterio.crs.CRS({"init": "epsg:4326"}), rasterio.crs.CRS({"init": "epsg:4326"}), img.size()[0], img.size()[1], gcps=[gcp1,gcp2,gcp3,gcp4])
 
    crs = rasterio.crs.CRS({"init": "epsg:4326"})
 
    profile.update(
        dtype=rasterio.uint8,
        transform = tsfm,
        crs=crs)
        
    with rasterio.open('example.tif', 'w', **profile) as dst:
        dst.write(src.read().astype(rasterio.uint8), 1)

This does produce an image but it is not properly warped. It is just in the correct location. Since my sensor is a 40 deg off nadir the warping should be reasonably significant, when I warp it with cameratransform's cam.getTopViewOfImage() function I get this:

example image

So I'm curious if I'm missing some rasterio function or ability, isn't what the GCPs should be doing in the transform.from_gcps() function?
 
Rasterio 1.2b2

Sean Gillies
 

Hi all,

Rasterio 1.2b2 is on PyPI now. We've made metadata reading more forgiving of unexpected values and changed the Cython language level for .pyx files to "3". See https://github.com/mapbox/rasterio/blob/master/CHANGES.txt.

Unless there are any problems discovered, we'll aim for a 1.2.0 in the middle of next week.

--
Sean Gillies
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?

nart@...
 

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?

nart@...
 

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?

nart@...
 

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
1 - 20 of 698  
More
Home Messages Hashtags Subgroups About Features Pricing Updates Terms Help Twitter