texts: set CRS for numpy array that has no geographical data using rasterio and bbox data

texts:
I have 1 dimension numpy array that has no geographical information. The numpy array was originally satellite image that I have manipulated in order to get this 1 dimension arrray.
I have the crs and transform of this array.
However, until now , in order to have geographic data on this array I had to write it as tiff and only then I could get back the geographical data and do things such as crop the raster and more.
I know there is option to reproject an array with no geographical data ,as mentioned here . I have tried to do this with my data as following,based on the doc:
src_shape=(886, 1138)
rows, cols = (886, 1138)
d = 1.0/240 # decimal degrees per pixel
# The following is equivalent to
# A(d, 0, -cols*d/2, 0, -d, rows*d/2).
src_transform = A.translation(-cols*d/2, rows*d/2) * A.scale(d, -d)
src_crs = {'init': 'EPSG:4326'}
source = img_pred

dst_shape = (886, 1138)

dst_transform = rasterio.transform.from_bounds(*bbox, width= bbox_size[1], height= bbox_size[0])
# this gets the transform which is in this case : Affine(0.00011618510158013974, 0.0, -44.75374, 0.0, -7.050087873462236e-05, -6.80845)

dst_crs = rasterio.crs.CRS.from_dict(init='epsg:4326')
destination = np.zeros(dst_shape, np.uint8)

reproject(
    img_pred,
    destination,
    src_transform=src_transform,
    src_crs=src_crs,
    dst_transform=dst_transform,
    dst_crs=dst_crs,
    resampling=Resampling.nearest)

Then I get new array with values 0:

so my question is : how do I get my array with the values reprojected correctly?
I assume that the problem might be in the definition of the original CRS, (src_crs), as my numpy array has no crs at all, but not sure if this is the reason.
In the end I want to reproject the array so I can clip it with shapefile and do  use other overlay tools.


texts: Enabling or disabling WMS capabilities with python

texts:
I am trying to enable Feature Access for all of my services on an ArcGISServer with python.  Below is my code.  The variablejson is the Feature Access dictionary.  How do I make the post targeting that value in the JSON.  The status code of my post returns 200 but when I look at the json of the service it remains false and when I log into manager it is not enabled.

url = "https://xxx.com/arcgis/admin/services/folder/serviceXXX.MapServer/edit"
json_open = requests.get(url)
json_load = json_open.json()

json =  json_load["extensions"][5]
key = json_load["extensions"][5]['enabled']
update = json_load["extensions"][5]['enabled'] = "true"
print json
print "__________"

post = requests.request("POST", edit_url, data={'token': token, key: update, "f":"json"} )
print "________"

print post.status_code

json structure returned from variable json:

{u'allowedUploadFileTypes': u'', u'enabled': 'false', u'maxUploadFileSize': 0, u'capabilities': u'', u'typeName': u'FeatureServer', u'properties': {u'allowGeometryUpdates': u'false', u'realm': u'', u'allowOthersToUpdate': u'false', u'dataInGdb': u'true', u'datasetInspected': u'false', u'allowOthersToQuery': u'true', u'enableOwnershipBasedAccessControl': u'false', u'editorTrackingTimeInUTC': u'true', u'xssPreventionEnabled': u'true', u'allowTrueCurvesUpdates': u'false', u'zDefaultValue': u'0', u'editorTrackingRespectsDayLightSavingTime': u'false', u'allowOthersToDelete': u'false', u'enableZDefaults': u'false', u'creatorPresent': u'false', u'syncVersionCreationRule': u'versionPerDownloadedMap', u'editorTrackingTimeZoneID': u'UTC', u'syncEnabled': u'false'}}


texts:
I can't find the option to delete a virtual field on a layer, not even in the table manager plugin.

Where is that option?


texts: PostGIS shapefile and DBF loader 2.0 problem in Windows

texts:
I am trying to import a shapefile with polygons into PostGIS (Windows 2008 server). If I import 100 polygons it works fine. If I try to import the full map of 8,100 polygons, about halfway I get an error window "Shape to PostGIS Converter has stopped working".
Problem details (as shown in that window) are:

  Problem signature:   Problem Event Name:  APPCRASH   Application
  Name: shp2pgsql-gui.exe   Application Version:    1.0.0.0   Application
  Timestamp:    505c6e19   Fault Module Name:   libglib-2.0-0.dll   Fault
  Module Version:   2.28.8.0   Fault Module Timestamp:  4e253544
  Exception Code:   40000015   Exception Offset:    0004c2d8   OS
  Version:  6.1.7600.2.0.0.272.7   Locale ID:   1033   Additional
  Information 1:    4a99   Additional Information
  2:    4a99e4ac53bd5e60be05464c3d1a217e   Additional Information 3:    794a
  Additional Information 4: 794abc993e1c90d5202dd7b66fde9bb9

Any suggestions as to what might be going wrong?


texts: Can i code a surface difference between two XYZ's in Python?

texts:
I would tend to utilize ArcGIS to perform surface difference between two TIN's or Terrains but my recent attempts have been too costly computationally. Instead, what I would like to accomplish is to perform a surface difference between two raw XYZ files and produce an output XYZ file. Is there a mathematical equation (via Python code) that I can use to generate this desired output? Obviously it would be focused on the difference in elevation (Z) of one file versus the other.

Once i have the new XYZ file, i would just proceed to create a TIN or Terrain of this XYZ file.


texts:
In a Python script, I'm cycling through thousands of Excel files and using arcpy.MakeTableView to create in memory table views from the first spreadhseet in each file.  How do you delete the reference to the in-memory view?  I've tried deleting the variable that stores the table name, but to no avail.

tblView = sheetName
arcpy.MakeTableView_management(excelSheet, tblView, .....)

## Process the table here

del tblView

The reason for deleting the reference is that in some cases the sheet names are duplicated in the Excel files.  If a table view has been created with the name "Sheet1", you can't create another table view with the same name.  In addition, due to the number and size of the Excel files being processed, I'm concerned about an impact on performance if all of those table views are preserved through the life of the script.


texts: Bad Layer in QGIS switching from PostGIS ST_DWithin from ST_Within

texts:
I can create a view in PostGIS using ST_Within and render it in QGIS without errors.

If I change to ST_DWithin (note the D) it creates a valid view in the database, but bad layers in QGIS.

This line works fine, no bad layers:
...cacities2014 d WHERE st_dwithin(c.geom, d.geom, 0.100) ...

This line causes the bad layers:
...cacities2014 d WHERE st_within(c.geom, d.geom) ...

Why would ST_DWithin have different characteristics in QGIS than ST_Within?


texts:
In geoserver i have a layer in EPSG:4326.
In OpenLayers i add it on map as WFS:

            markerLayer = new OpenLayers.Layer.Vector('Markers',{
        strategies: [new OpenLayers.Strategy.BBOX()],
        protocol: new OpenLayers.Protocol.WFS({
             version: "1.1.0",
            url:  "/geoserver/wfs",
            srsName:"EPSG:4326",
            featureType: "markers",
            featureNS: "http://www.opengeospatial.net/cite",
        }),
        projection:new OpenLayers.Projection("EPSG:900913"),
        styleMap: new OpenLayers.StyleMap({
            externalGraphic: url_servlet+'img/marker.png',
            graphicWidth: 20, graphicHeight: 24, graphicYOffset: -24,
            title: '${tooltip}'
        })/*,
        eventListeners: {
                beforefeatureadded: function(e) {
                    e.feature.attributes.annotation =
                        window.prompt("annotation", "no comment");
                }
            }*/
    });
app.mapPanel.map.addLayers([markerLayer]);

My map in OL in EPSG:900913. So when i says to WFS projection=900913 i exepcted that features reproject and will show on my map. But its show somewhere outside OSM extent. And in WFS layer i see features with coords like this:

    POINT(62.01759636987174 56.767130970984326)

So how reproject wfs correctly?

UPDATE

I change protocol to HTTP and get gml:

            protocol:new OpenLayers.Protocol.HTTP({
            url: "/geoserver/ows?service=WFS&request=GetFeature&typeName=cite:markers&srsName=EPSG:900913&outputFormat=GML2",
            format: new OpenLayers.Format.GML()
        }),

And now i get something interesting:
My layers and Google Base layer:

Same place with OSM base layer:

You see that there is no more red markers. And i have only question What a hell going on?


texts:
In the Feature Tutorial Example

GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory();

generates error:

  method getGeometryFactory in class
  org.geotools.geometry.jts.JTSFactoryFinder cannot be applied to given
  types;
        required: org.geotools.factory.Hints
        found: no arguments
        reason: actual and formal argument lists differ in length

Is this the right way?

GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory(null);


texts:
I am trying to get a raster file of nightlights data for Brazil and subdivide/regroup it into several rasters for each state, in a list or data.frame in R

I managed to download the nightlights file, import into R using the raster package, and crop it for a bounding box of Brazil. Reproducible code:

library(tidyverse)
library(R.utils)
library(raster)

#Source: 'https://ngdc.noaa.gov/eog/data/web_data/v4composites/F182013.v4.tar'
#Download this file manually ( untar() command does not work if I use download.file() to download)
fo <- 'C:/.../NightLights/'
f <- list.files(fo,full.names = T)
untar(f,exdir = fo)
tif <- paste0(fo,'/F182013.v4c_web.stable_lights.avg_vis.tif')
gunzip(paste0(tif,'.gz'),tif )

bra = getData('GADM', country="BRA", level=0) #country division
bra1 = getData('GADM', country="BRA", level=1) #state disions

F18_bra <- raster(tif) %>% crop(bra) #imports raster for a Bunding Box of Brazil

As the above code indicates I also imported the polygons for Brazilian states (object bra1)

I want to do a spatial join of the raster and the state polygons.

In the end, I want a data.frame, or list, with: state_id, state_raster. I.e., one line/element per state, with a raster object with the pixels inside the polygon of each state on each cell of state_raster (like a bucket or geom in Postgis)

Can I do this in R, preferably using the sp, sf and/or raster packages?


texts: Change the Default geodatabase in ArcGIS 10?

texts: In a multi-user environment I have a need to be able to configure ArcGIS so that it uses an SDE connection as the default geodatabase instead of the Default.gdb in each users' profile.  What is the best way to do this?  An add-in?


texts:
I've downloaded the service pack but have not yet installed it. Today I stumbled over a post on the ESRI forums where problems with reference scale and labeling are discussed. Apparently, this is now a known bug (NIM067942).

Normally, I would install a Service Pack to correct existing issues and to possibly improve performance but it seems that new bug(s) were introduced as well. This particular bug could potentially be a nuisance.

How buggy is it? Will I benefit from installing it? What are your experiences with this particular service pack?


texts:
I did a google search and a search in the ArcGIS forum. Does anyone know why when ArcGIS Desktop standardize an address such as "285 LAKE FRONT RD" and turn it into "Frontage Road"?

It's a fairly small problem but it really shouldn't happen. Right?

Edit: I'm using ArcGIS Desktop but I'm not sure why that would make a difference. I now realize that I should say that when I was standardizing the address field, I was using the US Single House locator.

The problem IS that "Lake Front Rd" is converted to "Frontage Road" after using the Standardize Addresses tool.

Frontage is not recognized by USPS as a street suffix.


texts:
I have a shapefile that has the outlines of the United States and I want to create hex grids that over that layer.

I tried using this function that will create hexes over a layer. I created an empty table and then added a geom column and then ran the custom function as shown.

CREATE TABLE IF NOT EXISTS hex_grid (hid serial NOT NULL PRIMARY KEY)
SELECT AddGeometryColumn('hex_grid', 'geom', 0, 'POLYGON', 2);

-- creating the hex grid
    SELECT genhexagons(0.05,ST_xmin(ST_Extent(SELECT geom FROM usa)),
    ST_ymin(ST_Extent(SELECT geom FROM usa)),
     ST_xmax(ST_Extent(SELECT geom FROM usa)),
     ST_ymax(ST_Extent(SELECT geom FROM usa)));

However I ran into this error that I'm not sure how to deal, I'm sure it deals with the creation of the hexes when it exceeds the US outline boundary. See below:


texts: ogr2ogr (or other) simplification of polygons with shared boundaries

texts:
I am trying to simplify a collection of 100,000 or so polygons containing nearly 4 million vertices.  The data is stored in a 2GB shapefile so there is no explicit representation of the topology between the polygons, but they consist almost exclusively of shared arc/vertices, and the coordinates of the shared vertices line up extremely well - to the point that I cannot see any difference no matter how far in I zoom in QGIS, at least in the many areas where I've looked.

I have tried ogr2ogr with the -simplify option, for doing the simplification, and although this preserves the topology of each polygon, it's not aware of shared arcs between polygons and so the result is that an arc that is shared between two polygons gets simplified in different ways for each polygon, resulting in overlaps and gaps between the resulting polygons.

Is there any such thing as an OGR-compatible data format that explicitly represents shared topology between polygons, so that I can convert this data into that into that format, so that when I simplify it, the resulting set of polygons has the same topology (and therefore has no gaps or overlaps)?

The reason I'm doing this is that I ultimately want to represent the (simplifed) polygons in a topoJSON file for web display in leaflet.js.
I've tried doing the conversion and simplification using the JavaScript topojson command, which does almost exactly what I want, except that it doesn't seem to be able to handle this large an input file. I've tried increasing the memory limit with the --max_old_space_size flag as described in https://github.com/mbostock/topojson/issues/71, to no avail.  I even tried this on a big AWS EC2 instance with 30GB of memory and it ran for a day with no output.  If I break the data up into smaller geographic regions and process each region separately with topojson, the results look wonderful, but then there are problems along the boundaries between the regions.  So I'm looking for solutions that can (preferably rapidly) handle the original huge data set.

I've seen references to the GRASS v.generalize command in related posts and am considering it as an option, but I've never used GRASS and am pretty familiar with the gdal / ogr tools and am wondering if there's a way to get this done with those tools, rather than having to climb up the steep GRASS learning curve.


texts: Creating Road Edge of Pavements from Road Centerlines Data using QGIS

texts:
Creating Road Edge of Pavements from Road Centerlines using QGIS
I have a Road Centerline data set and would like to use it to create a Road Edge Of Pavements data set that I could then export to a CAD format. Is there an easy way to do this in QGIS?


texts: How can I calculate decibel unit (dB) values for PALSAR in GoogleEarthEngine

texts:
I use the yearly Palsar mosaic of the GoogleEarthEngine.
In a first step I was able to calculate the HH-HV value and add it as a new Band to my image.
But now I would like to transform in a first step the 16-bit digital numbers (DN) to decibel unit (dB). Therefore, the following equation is provided by the GoogleEarthEngine:

γ₀ = 10log₁₀(DN²) - 83.0 dB

Does anybody know how I can calculate this transformation in the GoogleEarthEngine?

var dataset = ee.Image('JAXA/ALOS/PALSAR/YEARLY/SAR/2017');

// Compute the HH-HV value and add as new Band
var hhhv = dataset.expression(
  '(HH - HV)', {
    'HH': dataset.select('HH'),
    'HV': dataset.select('HV')
}).rename('HHHV');
var dataset = dataset.addBands(hhhv);


texts: I performed a spatial join based on point features that fell within polygons of interest. I now want to combine these features by polygon (field ID) and create a new field that combines information from one field. In the image below, I want to combine by TARGET_FID (ex: TARGET_FID = 7) and preserve the information from WS_Type to get a new field that would have, "Storage, Storage, Spring, Spring". I tried dissolve but the tool won't run because the "file is not supported".


texts: QGIS 2.8 Composer - extent / overview ... change style

texts:
I have a working overview setup with a general map and a detail map. Since scales are pretty much different - general 1:100.000 and detail 1:750 - I would like to change from the standard style (red, semi transparent rectangle) to a symbol ... like a crosshairs. has someone figured out how to change (i.e. fully customize) overview-map styles ?

To be more precise ... I would like to use symbols instead of any kind of fills!


texts: I have layer of points of which I need to build lines. Standart Points to path function doesn't suit me.  I have data on time, but it does not reflect the order of connection that I need. I need to connect them by distance. That is, if the distance between two points is less than 10 meters - draw a line. Is it possible?


texts: How to determine storage type (ST_Geom/Oracle) from arcpy/python

texts: We're migrating some feature classes from SDELOB (or some older binary format) to ST_Geometry.  As we write our arcpy script to use the ESRI Migrate Storage Tool we'd like to test whether or not something already is ST_Geom...Other than writing hooks into the database, using the comtypes package etc., does anyone know of a simple way to determine this?


texts: QGIS 2.8 Change Colour of Selected features

texts: How can I change the colour of the selected features in QGIS.  The light yellow is not easy to see.


texts: Why PgRouting NodeNetwork function doesn't create cost and reverse_cost columns?

texts:
I'm following the instrcutions here step by step;

Since my data is a custom road data that I couldn't find in Open Street Map I've followed these steps;

Road data in SHP format transferred to PostgreSQL using Postgis Shapefile Loader.
MULTISTRINGLINE data transformed to LINESTRING
Topology created for the data using pgr_createtopology function
Topology analyzed using pgr_analyzeGraph function
Errors found in analyze are fixed using pgr_nodeNetwork function
Topology recreated with new data table and analysis have been done similarly to first and second steps.

The problem starts from here; in documentation it says a column named old_id should be added to data table to keep track of original id's. Only segmented edges are added to table using the SQL statement;

insert into edge_table (old_id,dir,cost,reverse_cost,the_geom)
    (with
    segmented as (select old_id,count(*) as i from edge_table_noded group by old_id)
    select  segments.old_id,dir,cost,reverse_cost,segments.the_geom
            from edge_table as edges join edge_table_noded as segments on (edges.id = segments.old_id)
            where edges.id in (select old_id from segmented where i>1) );

When I run the statement above I get this error;

  ERROR:  column "cost" of relation "path" does not exist

In fact none of the functions I've runned while creating the nodes adds a cost column to data table. Should I add cost manually? If so how can this be done ? Documentation doesn't mention about cost. It sounds like it is automatically generated like the other columns such as source and target.


texts: Create ERDAS IMG Statistics using GDAL C++ API?

texts: I want to create very many ERDAS .img files using GDAL under C++. I don't want to have to process them all in ArcMap to create statistics for them (mainly because ArcMap fumbles when calculating statistics for very many input images in batch mode).  I want to create the .img files complete with statistics calculated for each.  I see python utilities for this, but nothing under the C++ API.  Does anyone have a solution?


texts: Counting number of points in particular polygon for list of individuals using R?

texts:
I have distribution points for a set of individuals and a file of polygons-with matching IDs.

What I would like to do is to go through and for each individual find the total number of points and the number which fall within the appropriate polygon, then have the output in spreadsheet form with a column for ID; Total number of points; points in polygon.


texts: Add PostGIS layer to QGIS via Python Console

texts:
How can I add a PostGIS layer to QGIS via the Python Console?

With some server setups, especially with a large number of tables, it can be very slow to load the list of layers via the Browser or Data Source Manager.


texts: How to transform feature geometry to 4326

texts:
I am working on WFS-T update.below is the code I am using to update modified data

   map.addInteraction(interactionSelect);
        interaction = new ol.interaction.Modify({
            features: interactionSelect.getFeatures()
        });
        map.addInteraction(interaction);
        map.addInteraction(interactionSnap);
        dirty = {};
        interactionSelect.getFeatures().on('add', function (e) {
            e.element.on('change', function (e) {
                dirty[e.target.getId()] = true;
            });
        });

        interactionSelect.getFeatures().on('remove', function (e) {
            var f = e.element;
            if (dirty[f.getId()]) {
                delete dirty[f.getId()];
                var featureProperties = f.getProperties();
                delete featureProperties.boundedBy;
                var clone = new ol.Feature(featureProperties);
                clone.setId(f.getId());
                transactWFS_polygon('update', clone);
            }
        });

before calling transactWFS_polygon function i need to transform projection to 4326.
The movement I modified feature its not visible.When I reload page also data is not visible.
My post request parameter:

<Transaction xmlns="http://www.opengis.net/wfs" service="WFS" version="1.1.0" xsi:schemaLocation="http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.1.0/wfs.xsd" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"><Update typeName="feature:polygon" xmlns:feature="http://spatial_editing"><Property><Name>gid</Name><Value>7</Value></Property><Property><Name>geometry</Name><Value><MultiPolygon xmlns="http://www.opengis.net/gml"><polygonMember><Polygon><exterior><LinearRing><posList>8355484.43590919 1768447.08640584 8267428.979324663 1719527.388303325 8370160.34533994 1684060.607179 8393397.20193863 1697513.52415719 8397066.17929632 1718304.39585076 8375052.31515019 1720750.38075589 8392174.20948607 1757440.15433277 8355484.43590919 1768447.08640584</posList></LinearRing></exterior></Polygon></polygonMember></MultiPolygon></Value></Property><Filter xmlns="http://www.opengis.net/ogc"><FeatureId fid="polygon.7"/></Filter></Update></Transaction>

Update1

Please check with jsFiddler same code i am using for WFS-T edit,delete and update.Every thing working except update.can any one tell me whats wrong here.layers are in local host so not displaying there


texts: Making Leaflet tileLayer grid starts at specified lat lon

texts:
My app needs to show a map within a box defined by

    const corner1 = L.latLng(5.21244812011719, -96.5979537963867);
    const corner2 = L.latLng(34.086555480957, -55.4220504760742);

    const map = L.map(el,
      {
        crs: L.CRS.EPSG4326,
        center: [19.6875, -76.201171875],
        zoom: 3.5,
        maxBounds: new L.latLngBounds(corner1, corner2),
      });
    L.tileLayer.wms("http://localhost:8080/wms?", {
        layers: 'xxx',
        tileSize: 700,
        transparent: true,
        maxBounds: new L.latLngBounds(corner1, corner2),
        bounds: new L.latLngBounds(corner1, corner2)
      }).addTo(map)

The thing is that the layer is requesting the tiles starting outside the borders defined (for example, this request was made http://localhost:8080/wms?&service=WMS&request=GetMap&layers=Mano%3AT22&styles=&format=image%2Fjpeg&transparent=true&version=1.1.1&maxBounds=%5Bobject%20Object%5D&time=2018-05-02T00%3A00%3A00.000Z&width=700&height=700&srs=EPSG%3A4326&bbox=-87.71484375,28.4765625,-56.953125,59.23828125 the bbox was -87.71484375, 28.4765625, -56.953125, 59.23828125, this contains the corners defined by me but is not optimal as it is requesting a map area that is not going to be seen on my app)

I need that leaflet request tiles that do not exceed the corners defined.

How can I archieve that?


texts: A library/software to create DEM from a topographic map

texts:
In one of my projects I need to digitize video clips showing some natural landscape. Specifically I need to know distances between different areas and elevation levels.

I have a topographic map of the area so distances are not a problem.

Is there any tool / library that allows to create DEM from a common topographic map?
A language of the library doesn't matter, but I prefer open source solutions if possible.

I know about QGIS and GRASS, but I cannot found how to use them for this specific purpose.

EDIT:
My map looks like this:


texts: Overcoming zoom level problems for layer offsets in QGIS

texts:
In this question I was able to align some labels with polylines representing hectometer points along a railroad. I did this using custom expressions in the "data defined" portion of the labelling tab in the layer properties:

The problem now is that the label spacing with respect to the line end points is zoom dependent. Here's a zoomed in portion of the map:

Here's the same portion of the map, zoomed out:

This is kind of a problem when it comes to making layouts, as usually you want to have separate layouts coming from the same canvas to have slightly different zoom labels but to be consistent when it comes to symbology.

Perhaps one way to overcome this issue is by specifying canvas units for the offset rather than map units in the X and Y portions of the data defined section in the Layer Properties Dialog:

X: x(end_point($geometry))+5*canvas_units
Y: y(end_point($geometry))+3*canvas_units

But I'm not sure whether this is possible. If anyone has any idea how to do this please let me know. Thanks.


texts: Connecting Python script external to Grass GIS 7 program in Windows 10?

texts:
I've tried to find resources online to link and run my Python script externally from Grass GIS 7, however, I have not been successful. I found resources here that provide instructions (e.g., https://grasswiki.osgeo.org/wiki/GRASS_and_Python#MS-Windows), but this is my first time using Python, and I haven't been able to figure it out. I am using Grass GIS 7.0.4, iPython 3.5, and Windows 10.

I tried to set environment variables in iPython and IDLE (don't know if it matters which I use to link it with Grass) with the following code:

import os
import sys

gisbase = os.environ['GISBASE'] = 'C:\OSGeo4W\apps\grass\grass-7.0.4'  #GISBASE needs to point the root of the GRASS installation directory
gisrc = 'GISRC= C:\Documents and Settings\Suzaku1\.grassrc6'
gisdbase = 'C:\Users\Suzaku1\Documents\grassdata'
location = 'demolocation'
mapset = 'PERMANENT'
LD_LIBRARY_PATH = 'C:\OSGeo4W\apps\grass\grass-7.0.4\lib'
PATH = 'C:\OSGeo4W\apps\grass\grass-7.0.4\etc';'C:\OSGeo4W\apps\grass\grass-7.0.4\etc\python';'C:\OSGeo4W\apps\grass\grass-7.0.4\lib';'C:\OSGeo4W\apps\grass\grass-7.0.4\bin';'C:\Python27';'C:\OSGeo4W\apps\Python27';'C:\OSGeo4W\apps\grass\grass-7.0.4\scripts';'C:\OSGeo4W\bin';'C:\OSGeo4W\apps\Python27\Lib\site-packages'
PYTHONLIB = 'C:\Python27'
PYTHONPATH = 'C:\OSGeo4W\apps\grass\grass-7.0.4\etc\python'
GRASS_SH = 'C:\OSGeo4W\apps\msys\bin\sh.exe'

sys.path.append(os.path.join(os.environ['GISBASE'], 'etc', 'python'))

import grass.script as grass

The gisrc variable links to a text document I created following the recommendation here (see Update#2). I saved the text file as 'grass.pth' and placed it in (C:\OSGeo4W\apps\python27\lib\site-packages). Not sure if my gisrc variable is correct here since it doesn't link to the location referenced in the code above: Problem with python script to control GRASS GIS from outside - How to import grass.script under Win 8.1?

When I run the code, I receive an error saying that, 'import grass.script as grass' importerror: no module named grass.script'.

What am I doing wrong?

I have iPython (v3.5) installed via Anaconda and a separate installation of Python that was provided with the Grass GIS 7 download (v2.7). Not sure if this matters for linking.


texts:
Overpass has some field names to store certain meta attributes of an object to a CSV file, including its coordinates.
However, bounding box coordinates are not among them.
When running a query in out bb mode, how can I retrieve the coordinates of the bounding box for each object returned?


texts: Translating tif rasters to asc taking a long time

texts:
I need to use the 19 WorldClim variables in MAXENT, so I have to take the .tif files to .asc files using the translate tool. It's taking forever for me (2% after 2 hours). I've given up. I know its possible, and people use WorldClim for MAXENT, is there a trick I'm missing?

Edit to add: My rasters vary in size from 1 GB to 512 MB. Obviously these are large files, and I expected them to take some time, but maybe not upwards of an entire day. Curiously, when I clipped them to a smaller area of interest, the files came out larger than the untouched version.


texts: Is there an R function for creating neighbors list from polygons based on both proximity and adjacency (contiguity)?

texts:
I am trying to perform hotspot analysis on non-contiguous US shapefile (48 states plus Alaska and Hawaii) using Getid Ord G stat (at the county level). In order to use getid Ord G, I need to create 1. neighbors list and 2. weight matrix for those neighbors.
neighbors list can be created by two methods: 1. proximity, using a distance band ("dnearneigh" function and then "nb2listw" function) 2. Adjacency, using borders shared between units ("poly2nb" function and then "nb2listw" function)
Adding Alaska and Hawaii causes two issues. First, there are several very large unit areas in Alaska (it makes proximity an inappropriate method). Second, there are unit areas in Hawaii without any neighbors (it makes adjacency (contiguity) an inappropriate method).
I need to define neighbors list based on both proximity and adjacency. I am looking for a function that defines neighbors based on a distance band, and if there are any cases that do not fall in the distance band (like the ones in Alaska), defines neighbors based on adjacency.
Thank you.


texts: Doing spatial join which is also based on time using ArcGIS for Desktop?

texts:
I have a time series plot of house sales for 4 years, but in those 4 years some districts change. I was wondering if I need to split the data into 4 layers for each year in order to join them to their respective district at the time of the sale.

For example, house 1 sells 2 times, once in 2007 and again in 2010. During that time the district changed from district 229 to 271. Is there any way to join them so that the right observation takes the right district identifying when doing a spatial join?

Districts are polygons and addresses are points. The address points are geolocated based on address listed in the assessor information based on house number, road name, road direction, and zip code.

I am using ArcMap 10.1.


texts: Every time I pan the map content in the QGIS composer window the scale changes. Is there a way to make it stick at specific scale?


texts:
I have downloaded and installed GeoServer 2.1-RC1 which shipped with RESTconfig. How can i check whether REST service is installed with GeoServer or Not ?

I put this "http://localhost:8080/geoserver/rest" in my browser, but it gives error! ... It means REST service is not installed as i understand. Or is there any other method to check whether the service is installed or not ?


texts: Create feature datasets from group layer name

texts:
Is it possible to create feature dataset, based on the name from group layers. I have two separate codes, one for finding out the names of group layers, and other to create feature datasets.

>>> import arcpy
    mxd = arcpy.mapping.MapDocument("CURRENT")
    for lyr in arcpy.mapping.ListLayers(mxd):
    if lyr.isGroupLayer == True:
    if lyr.longName.find('\\') == -1:
    print lyr.longName

This is the code for searching the name of datasets, and this is the code to create feature datasets.

>>> arcpy.env.overwriteOutput = True
    fdList  = ["Dataset_A", "Dataset_B", "Dataset_C"]
    folList = ["D:\\GIS_Temp\Folder A", "D:\\GIS_Temp\\Folder B", "D:\\GIS_Temp\\Folder C"]
    workRange = range(len(fdList))
    for thisIndex in workRange:
    fd = fdList[thisIndex]
    arcpy.env.workspace = folList[thisIndex]
    arcpy.CreateFeatureDataset_management("D:\\GIS_Temp\\TEST.gdb", fd, "D:\\GIS_Temp\\Projection.prj")
    for impFC in arcpy.ListFeatureClasses():
    fcName,fcExt = os.path.splitext(impFC)
    fcName.replace(" ","_")
    arcpy.FeatureClassToFeatureClass_conversion(os.path.join(folList[thisIndex],impFC),os.path.join("D:\\GIS_Temp\\TEST.gdb",fd),fcName)

Is it possible to merge the two codes, into one, so the name fdList  = ["Dataset_A", "Dataset_B", "Dataset_C"], would actually be names returned from the group layers based on the first code? I am trying to merge these two codes, so the searched names would actually be names of datasets.


texts: Is Dynamic Traffic Format supported by QGIS?

texts: I have to merge (Dynamic Traffic FORMAT) DTF file into QGIS. This data is only available from NAVTEQ or TOM TOM. When I have searched ArcGIS network analysis tool, there is a clear option for adding a DTF file. I don't see the same feature for QGIS. I would like to know if there is a similar option in QGIS.


texts: Show previously loaded tiles while loading (zoom in or zoom out) map

texts:
How to remove the white tile part(map loading area).

While zoom in and zoom out, the white tile part is appearing on the map(map loading area). How can we fill that area with the history image?

Is there any better way to load maps faster than the GeoWebCache.


texts: 'native' vs. 'qgis' geoprocessing algorithm provider in (Py)QGIS 3

texts:
In QGIS 3 processing toolbox among all the algorithms with the QGIS 3 icon there are some with Algorithm-ID starting with qgis: and some starting with native:.

What is the difference between processing algorithm providers qgis and native?


texts: How do I rotate a single Landsat imagery, so it fits the frame of image?

texts:
I want to generate a JPEG from a GeoTIFF file which has a stack of RGB bands. Considering the JPEG won't necessarily have any geo-referencing information about any pixels, how do I get rid of the black borders? I want the data-pixels(the colored ones here) to fit the frame as a rectangle.
I have experimented that the angle of inclination increases as we go to the polar regions from equator.


texts: Using Create Random points ignores number of points required in ArcPy

texts:
I am trying to generate around 117 random points over patchy forested landscape. As my forested landscape contaisn many polygons (28.000), and I would like to generate specific number of random points, I first dissolve my polygon, to create multipart polygon, as speciefied in ArcGIS help.

Example of my very patchy landscape, polygons are forests, white parts are not forested:

So, to generate random points, I am using as constraints a multipart polygon and minimal distance points.

However, no matter what required number of points in number_of_points_or_field I specify (should by 117 if I quess correctly), I always get back number around 75.

I understand that if my points cannot fit the polygon under minimal distance restriction, there is less of them.

But, I got back the same number of 75 points, even if I set lower number, i.e. 10: number_of_points_or_field = 10.

Why does this happen?

Here is my code:

# Read input forest data
inForest = os.path.join(inWD, "patchyForest")

# Process: Dissolve feature to create single multipart polygon
outForestDiss = os.path.join(outDB, "forest_diss")
arcpy.Dissolve_management(in_features = inForest ,
                          out_feature_class = outForestDiss,
                          dissolve_field = "",
                          statistics_fields = "",
                          multi_part = "MULTI_PART",
                          unsplit_lines = "")

# Process: generate random points within the forested area
outPtsName = "randomPoints"
numberPoints = 117  # tried with 10, 75, 117 and I always get back number 75
arcpy.CreateRandomPoints_management(out_path = outDB,
                                    out_name = outPtsName,
                                    constraining_feature_class = outForestDiss,
                                    # constraining_extent = "", # skip this one out
                                    number_of_points_or_field = numberPoints,
                                    minimum_allowed_distance = distanceBuff * 2,
                                    create_multipoint_output = "POINT")


texts: How to close and fill Gulfs and inlets in QGIS or Grass?

texts:
I have a raster layer from which I extracted only zones with a certain color, then polygonized it with QGis. I just generalized the new vector layer in GRASS to smooth things a bit.
There are many holes, "gulfs", "inlets" and cuts because the color zone in the raster overlays a black and white cartographic map. I would like to fill this spots, but I can't find a easy way. I already closed valid holes, but with "gulfs" and "inlets" I cannot use the same command.

In that image, only the hole on the right is suposed to exist. Is there a easy way to solve this?


texts: Web map graphics and performance optimization

texts:
With a web map (in this case, Bing Ajax API v7), what's more important for performance - the number of individual map shapes/graphics drawn on the map, or the total number of vertices drawn on the map?  Put another way, does a single map graphic with 100 vertices use the same amount of resources as 50 map graphics with two vertices each?

Thanks.


texts: Add default classification to Esri ASCII Raster?

texts:
I have some .ASC data files from Wales Govt. LiDAR data.  I converted it to LAZ using LasTools, but when I try to use the data in a third party app, I get an error saying that the data is unclassified.
I am good with data files and processing, but I do not know how to change the ASC data to make it classified.
Could anyone direct me on how to add some default classification to the data?


texts:
First time dealing with cursor.
I have lines intersecting polygons, I want to get a table with all attributes of the lines, add a field and calculate the field as the name of the polygon.
I was able to get a table for each polygon but I can't get to calculate the field and add the polygon name.
How do I do?
here is my code:

import arcpy
arcpy.env.overwriteOutput = True

arcpy.env.workspace = 'C:\\Users\\a391665\\Documents\\ArcGIS\\Default.gdb'
ws = 'C:\\Users\\a391665\\Documents\\ArcGIS\\Default.gdb'
Comuni = 'C:\\Users\\a391665\\Documents\\ArcGIS\\Sardegna3.gdb\\Comuni'
Linee = 'C:\\Users\\a391665\\Documents\\ArcGIS\\Sardegna3.gdb\\Linee'

rows = arcpy.SearchCursor(Comuni)

for row in rows:
    print row.NOME
    feat = row.Shape
    arcpy.SelectLayerByLocation_management("Linee", "INTERSECT", feat, selection_type="NEW_SELECTION");
    out_feature_class = arcpy.ValidateTableName("Linee_Clip_" + row.NOME);
    table = arcpy.TableToTable_conversion("Linee", ws, out_feature_class);
    arcpy.AddField_management(out_feature_class, 'Clip', 'TEXT', '', '', '60');
    expression = arcpy.ValidateTableName(row.NOME);
    arcpy.CalculateField_management(table, "Clip", expression);

here is the error I get:

  Empty value for ObjectID = 1
  The calculated value is invalid for the row with ObjectID = 1. For example, the calculated value may be too large for the field or you may be trying to add a string to a number field. This row will not be updated.


texts: Creating raster with defined value range in R

texts:
I want to create a raster with min value -1 and max value of 1. The values can be distributed in any form and the extent and resolution can also be random. I just want my raster to have values between -1 and 1.

I want this for purpose of symbology.


texts:
https://maps.six.nsw.gov.au/arcgis/rest/services/public/NSW_5M_Elevation/ImageServer/WMTS/1.0.0/WMTSCapabilities.xml is some elevation data of interest to me.

In QGIS, I've set that up as a WMTS layer, but what do I need to do to visualize it, for example as grey scale?


texts: Why does QueryFeatureIDs return an empty array when querying a Route Event Source?

texts:
I'm struggling with using the ESRI SOAP api to perform a spatial filter against a layer with a datatype of 'Route Event Source'

Because the number of features returned by a single spatial query could easily exceed the configured ArcGIS server max, I first grab all the features Ids of the features within the selection envelope, and then retrieve the feature attributes in chunks.

No matter what I do, the array of Feature Ids returned by QueryFeatureIds is always empty.  With the same SpatialFilter, though, I receive a correct count from QueryFeatureCount and I receive the correct records from QueryFeatureData.

Here's the test case that illustrates the problem I'm having.  Am I missing something, or is that the intended behavior? (I've changed the selection envelope in the sample)

var proxy = new MapServerProxy("http://mymapserver/ArcGIS/services/myservice/MapServer");

var envelope = new EnvelopeN();
envelope.XMax = 100010.0;
envelope.XMin = 100000.0;
envelope.YMax = 200010.0;
envelope.YMin = 200000.0;

var spatialFilter = new SpatialFilter();
spatialFilter.SpatialRel = esriSpatialRelEnum.esriSpatialRelIntersects;
spatialFilter.GeometryFieldName = "Shape";
spatialFilter.FilterGeometry = envelope;

int routeEventLayerIndex = 0;
var count = proxy.QueryFeatureCount("Layers", routeEventLayerIndex, spatialFilter);
Assert.AreEqual(2, count); // PASSES, 2 features within envelope

var records = proxy.QueryFeatureData("Layers", routeEventLayerIndex, spatialFilter);
Assert.AreEqual(2, records.Records.Length); // PASSES, data for 2 features is returned

var fidSet = proxy.QueryFeatureIDs("Layers", routeEventLayerIndex, spatialFilter);
Assert.AreEqual(2, fidSet.FIDArray.Length); // FAILS!!! Length of FIDArray is 0

Edit
I've tested this against ArcGIS Server 10.0 and 10.0 sp 4; both return the same results.


texts: OSM all relevant keys and tags to extract Features for routable maps

texts:
My aim is to create a routable network layer with the openstreetmap editor from Esri.
First of all, I need to know all relevant tags and keys for highways and cycleways to extract the data. Even pedestrian ways.
At the wiki of openstreetmap, there is a list of all tags, but what about the keys. I cannot find a systematic overview.

Does anyone know (almost) all tags and features who should be considered?


texts:
I have a html page with 2 WMS layers.
The first is one Base Layer , & the second is my personal WMS feature. I load it from Geoserver. The problem when I add the 2 layers I do not have the superposition

The problem is maybe on the projection. please can you help me this is my code

<script src='http://maps.google.com/maps?file=api&amp;v=2&amp;key=ABQIAAAAjpkAC9ePGem0lIq5XcMiuhR_wWLPFku8Ix9i2SXYRVK3e45q1BQUd_beF8dtzKET_EteAjPdGDwqpQ'></script>
  var map;
        OpenLayers.ProxyHost = "/cgi-bin/proxy.cgi?url=";
        var host = "http://'my_host':8080/geoserver/wms";

        function init() {

            map = new OpenLayers.Map('map', {
                    controls: [
                        new OpenLayers.Control.Navigation(),
                        new OpenLayers.Control.PanZoomBar(),
                        new OpenLayers.Control.LayerSwitcher({'ascending':false}),
                        new OpenLayers.Control.ScaleLine(),
                        new OpenLayers.Control.MousePosition()]});

            var reg = new OpenLayers.Layer.WMS(
                "Région",
                host,
                {layers: 'cite:matregion_view', format: 'image/gif',transparent: true}
            );

               var gmap = new OpenLayers.Layer.Google(
                "Google Streets", // the default
                {numZoomLevels: 20}
            );

 map.addLayers([gmap , reg]);
//map.setCenter(new OpenLayers.LonLat(146.7, -41.8), 6);
        map.zoomToExtent(new OpenLayers.Bounds(-16.884,20.748,-1.026,35.924));

        }


texts:
I am working in ArcMap wanting to rotate a parcel in parcel fabric.

NOT working in parcel fabric you can click "S" for another rotating point to come up.

However this does not work in parcel fabric. I can rotate the parcel, however its not letting me create another rotate point


texts: ArcGIS Pro VectorTiles export style json placeholders/variables error

texts:
When we export our vector-tiles from ArcGIS Pro we have what seems to be placeholders or variables in certain layout properties:
"text-rotate": {
     "property": "_label_angle2",
     "default": 0
 },
 "text-field": "{_name2}"

When uploading and checking the validity of the style.json in Maputnik these get marked with missing required property "stops". I can manually fix this error, but I wonder why these exist in the first place and how to avoid them.


texts: Compilation from source for QGIS 2.2.0 in OpenSuse

texts:
Someone can help me for installing from source QGIS valmiera 2.2 on linux OpenSuse 12.2.

I needed of understand the library and relative version of a tested enviroment.

The version of OpenSuse is not influent but I needed to use valmiera 2.2 version.


texts: Why does FME Reprojector create coordinate values 0?

texts:
I'm using the NeighborFinder transformer in FME Desktop 2013, and I need it to use a maximum distance in meters. My input data, however, are points and lines in WGS84, so before sending them into the NeighborFinder, I'm using Reprojector to convert the geometry into FME's dynamic coordinate systems AZMED or AZMEA.

The lines reproject okay. The points, however, all reproject to 0,0. So when the points go to NeighborFinder, they create meaningless output.

Here's a part of the setup before NeighborFinder:

Using Inspection, I can see a point before entering the Reprojector:

And after leaving the Reprojector:

Note the loss of geometry.

I've tried leaving the Source Coordinate System as "Read from Feature", but still no gold.

The Transformer Description for Reprojector gives an example of reprojecting a point to AZMED or AZMEA. And I can't find any indication online that points aren't supported.

I can reproject the points by changing the Destination Coordinate System to the relevant UTM zone (EPSG:32617, UTM Zone 17). But I'd prefer a scalable solution because I have global data to feed into the workbench.

Does anyone know why I can't reproject my points into AZMED or AZMEA?


texts: What are technical skills that a GIS graduate should develop?

texts:
I am a GIS Analyst with one years experience and no idea where to begin with professional development. I'm competent with the ArcGIS platform, but I am not too familiar with its administration processes. I can work with basic Python programming. I'm taking a course in SQL next week, and I'm pretty terrible with web services.

Where/How/What should i focus on learning to step up to the next level in a GIS career.


texts: How can I determine the layer name in an edit session?

texts:
I have a problem with an editor extension in an Add-In for ArcGIS 10.
The editor extension listens on different events occuring on shapefile layers via registering listeners:

mEditEvents.OnChangeFeature += OnChangeFeature;
mEditEvents.OnCreateFeature += OnCreateFeature;
mEditEvents.OnDeleteFeature += OnDeleteFeature;

All these callbacks get are a parameter of type IObject, which can be cast to IFeature.
With this IFeature I can get the feature class and its alias name (IFeature.Class.AliasName) which seems to be the name of the shapefile associated with this feature, but I want the name of the layer associated with this feature which can be different from the shapefile's name.

Is there any way how to get this layer name?

Thanks in advance!

Update:

I found a solution (see my comment below) - thanks for the help!


texts: select random rows with python in arcgis

texts:
I am new to using Python in ArcGIS, and I wonder: how do I select random lines (or every two lines) of an attribute table? I must select these lines and calculate them, so they are different from other.

I actually have some conditions. For example: I will have to calculate the area of a polygon and the area is less than 30 for example, I have to select one random point within that polygon and calculate his column so that it is different the other.

Part of calculating the area is already ready, lacking only get select random points.


texts: Finding distance between buildings and nearest road using QGIS

texts:
I am very new to working with GIS data. I am currently using the QGIS tool to work with the data.
I have two datasets:

Buildings - polygon shapefile with information about different buildings in a city. There is a unique column called "propertyid".
Roads - line shapefile with information about different roads in the city. There is no unique column, but we do know that each line item refers to one road, so we have created an index column based on "@row_no".

Objective: To find the distance between a given building and the nearest road.
What I have already tried:

Converted buildings layer to points(centroids) and lines layer to points(extracted vertices) and run a distance matrix using these two layers. But this doesn't give me the desired result for some reason.
Since I am very new, I could be doing something wrong.
Using the same point layers mentioned above, ran the distance to the nearest hub algorithm. Again, I may have committed an error in doing this as well.

It would be great if someone could guide me in the right direction and tell me if

what I'm trying to achieve is possible, and
if it is, then how would you go about doing it.


texts: Can I interpolate but follow the trend of another surface?

texts:
I have a point data set will a fairly low sample density and a raster that should loosely correlate with the values of the points. For example, this could be elevation and temperature readings.
I want to create a surface that interpolates the points, preserving the values at the points where I have data, but in between follows the trend from the existing raster. A bit like is this very bad drawing.

Has anyone done this before and know how to do it. I suppose I could interpolate the points to start with (spline will preserve my values so I'll use that. then add/multiply version of my raster that has been normalized so that its 1 at each point then + or - the size of difference from the average at each point. I'm not sure if this would work, if it would be the best method, or how to produce the normalized raster. Does anyone have any ideas?


texts: Importing existing .gpkg files to create new Geopackage layer using QGIS

texts:
I'm a new QGIS user running 3.14.0 on a Mac with Mojave 10.14.6. While following along with this tutorial (https://www.youtube.com/watch?v=kCnNWyl9qSE), I reached around the 22 minute mark and encountered an issue.
When attempting to create a "New GeoPackage Layer" by importing an existing .gpkg file I've already created, I click "..." next to Database expecting to be able to browse my computer and import this file. However, while I can see the file I want it is grayed out/not able to be selected. The window that pops up says "Save As..."; it seems like my only option is to create and save a brand new GeoPackage layer rather than create a new one using the existing .gpkg file like the tutorial does. What I need to be able to do is select "MyFirstData.gpkg"  (image below), but it won't let me.
Especially since I'm so new to QGIS, I feel clueless about how to fix this problem and I can't continue the tutorial unless I can manage to successfully import and existing .gpkg file.
Why does my system work differently than the one in the tutorial and why can I only "Save As" rather than upload existing files?


texts: OpenLayers 3 insert layer below other layers

texts:
I have this JSFiddle that can turn the OSM data on and off, the only issue being when it adds it back in the layer is added on top of everything else.

I know this is a relatively simple fix using the insertAt from the maps layer collection but I can't seem to get it working with this example.Perhaps i'm placing the map.insertAt in the wrong place?


texts: Running buffer algorithm from processing on each feature from one layer in pyqgis?

texts:
In a standalone outside script for QGIS 2.18, i have a layer which contain multiple features to bufferize :

vlayer = QgsVectorLayer(out_dir + "selected_shape.shp", "blabla", "ogr")

I try to run processing for each features

from qgis.core import *
from processing.core.Processing import Processing
Processing.initialize()
from processing.tools import *
import processing

buff_res = processing.runalg("qgis:fixeddistancebuffer", vlayer, 0.2, 99, False, None,)

get_output= processing.getObject(buff_res['OUTPUT'])

for f in get_output.getFeatures():
    print f.attributes()

The qgis:fixeddistancebuffer compute the buffer only on the first feature from my layer. How can i compute the buffer on each feature on this vlayer ?


texts: Leaflet+GeoJSON: Why are lat-lng coordinates in the wrong order?

texts:
I'm coding an app with Phonegap, which requests user's current position. As an overlay I have loaded a geojson layer. I want the app to automatically extract a value of said layer at the exact location, where the user currently is (without the user manually clicking on the layer).

I've tried Leaflet PIP extension, but couldn't get it to work with leaflet 1.0.3. I've also tried simulating a touch/mouseclick event at screen coordinates, but that didn't work either.

I've found the extension Leaflet.PointInPolygon, which can check, if a point is lying within a polygon. Therefore a transformation from geojson to an array with coords is required.

I wrote the following code to parse my geojson and write it into arrays to add polygon layers to my map (in order to check later, if map.getCenter() coords are within any of these):

var coord = [];
for (var i = 0; i < data.features.length; i++) {
  feat = data.features[i];

  for (var j = 0; j < feat.geometry.coordinates.length; j++) {
    coord = feat.geometry.coordinates[j];

    var polygon = L.polygon(coord, { color: "red" }).addTo(map);
  }
}

I'm facing the prolem now  that lat and lon are in the wrong order:

[[10.77932197941255,52.03584196145097],[10.77932197941255,52.03575906909681],[10.779636389761993,52.03575906909681],[10.779636389761993,52.03584196145097],[10.77932197941255,52.03584196145097]]

How can I fix this?

Is there maybe a leaflet extension or javascript framework, which can directly check, if given coordinates are within a geojson object?


texts:
I am trying to install cartopy on OSGeo4W.  Even though all required packages  already downloaded and updated, the installer won't recognize them.  I keep getting the following error:

setup.py:168: UserWarning: Unable to determine GEOS version. Ensure you have 3.3 or later installed, or installation may fail.
setup.py:224: UserWarning: Unable to determine Proj4 version. Ensure you have 4.8.0 or later installed, or installation may fail.
Proj4 version 0.0.0 is installed, but cartopy requires at least version 4.8.0.

I also tried downloading an older version of cartopy and attempted to build it against the installed shapely and GEOS.  I received a long dump of output, including the error:

*LINK : fatal error LNK1181: cannot open input file 'proj.lib'
error: command "C:\Users\tug28727\AppData\Local\Programs\Common\Microsoft\Visual
C++ for Python\9.0\VC\Bin\amd64\link.exe"' failed with exit status
1181*

Any ideas?


texts: Importing .geotiffs and wavelet denoise plug-in - GIMP

texts: I am trying to run the wavelet denoise plug-in tool on GIMP to smooth a .geotiff. Unfortunately when I import the .tiff, the spatial refernece is lost and the image appears black. My goal was to smooth away uncorrelated roughness while preserving breaklines. Does anyone have any experience with this plugin or importing .tiff's in GIMP.


texts: I am trying to reconcile and post the edits. I would like to know if I can post the versions to two different parents. Could some one help me with c# coding on how to do this.


texts:
I try to use the processing 'qgis:extractbylocation' in QGIS 2.8.2 from the QGIS algorithm provider

  #Récupérer tous les départements qui s'intersectent avec la zone d'étude
    res=processing.runalg('qgis:extractbylocation', departementLayer, layerZone, u'intersects', 0,None)
    #Créer une couche vector issue de l'extraction des départements dans la zone
    departementLayer=QgsVectorLayer(res['OUTPUT'], "Departement(s) dans la zone d'etude", "ogr")
    # QgsMapLayerRegistry.instance().addMapLayer(departementLayer)
    #Définir un ecodage Windows compatible avec le français
    departementLayer.setProviderEncoding(u'iso-8859-1')
    #Modifier l'encoage
    departementLayer.dataProvider().setEncoding(u'iso-8859-1')
    #Vérifier que le couche est valide sinon afficher que la couche est erronée
    if not departementLayer.isValid():
        print "Layer failed to load!"
    #QgsMapLayerRegistry.instance().addMapLayer(departementLayer)
    #Modifier la barre de procession à 11%
    self.dlg.progressBar.setValue(11)

however when, i run the application , it gives not just the districts(departemets) inside the layerzone but also those that are next to the layerzone . My question is : what processissing can i use to have just the ditricts(departement)inside the layer zone ?


texts: 'use estimated metadata' doesn't work in QGIS 2.14, 2.16; metadata queries of death

texts:
I have a few huge PostGIS tables, and despite setting 'use estimated metadata' when creating the underlying datasource, it STILL issues these incredibly frustrating 'queries of death' at both layer addition and loading saved projects causing full table scans:

SELECT count(distinct ("block_id"))=count(("block_id")) FROM "christest"."block_th_view"
SELECT st_extent("the_geom") FROM "christest"."block_th_view"

Can't we just add a bounding box metadata tag on a layer definition, and trust that it is indeed using a unique primary key, and optionally have QGIS just blindly trust the project file on start-up?

Consider a well-indexed table with billions of GPS points with a free running counter as its unique key, that only allows display up to a zoom scale of 1:1000. Good GIST indexing and fast RAID disks make it very usable ... once given a chance... BUT..not one, but TWO (!!!) full table scans are required each time to load the layer, that yields the same obvious results: one billion unique keys and the entire world. It is an utter waste of both the user's and the database's precious resources to recalculate this, much as it would be to recalculate pi to a million places every time you wanted to use that function on your calculator.

feature request:

Add a new tag "boundingBox" or "extents" child XML node under a maplayer, within QGIS project files, that is blindly trusted when a global preference is set, rather than "figuring it out" on the fly by performing a full table scan with SELECT st_extent() during load, and optionally skip the unique key check.  Furthermore, add an additional dialog box to optionally  set this manually when the layer is first added, as is done in GRASS GIS.

This will make QGIS -far- more useable for very large database layers.


texts: Adding Feature Class to Terrain using ArcPy gives ERROR 000800?

texts:
I'm writing a script to convert lidar data (in .las format) into a terrain without going through all the clicking and so I can send it out to neighboring states with similar programs (I'm using ArcGIS 10.1 btw). I've gotten almost there with the script but issues arise when I try adding the Multipoint to the Terrain in the Feature Dataset.  I've done this step in model builder and exported that to a python script and am using the "in_features" that's output from the model.  Incidentally, the model runs just fine in arcmap, but when I try to run the exported python script I get errors.

import os
import arcpy, math, glob
from arcpy import env

arcpy.CheckOutExtension("3D")
arcpy.env.overwriteOutput = True
#make workspace in memory
arcpy.env.workspace = "in_memory"

#create a feature dataset for export
#change geodatabase to a getparameter
geodatabase = r'C:\Users\peter.a.metzger\Desktop\Projects\LiDAR\test_data\lidar_test.gdb'
outgeodatabase = 'test_lidar_output'
projection = r'C:\Users\peter.a.metzger\AppData\Roaming\ESRI\Desktop10.1\ArcMap\Coordinate Systems\NAD83_UTM_zone_14N.prj'
created_FC_tool = arcpy.CreateFeatureDataset_management(geodatabase, outgeodatabase, projection)
created_FC = created_FC_tool.getOutput(0)
#Create a Terrain

terrain_name = 'test_Terrain' #Change to a getparameter
arcpy.CreateTerrain_3d (created_FC, terrain_name, 1.24605815948, 50000, '', 'WINDOWSIZE', 'ZMIN', 'NONE', 1)

output_terrain = os.path.join (created_FC, terrain_name)
arcpy.AddTerrainPyramidLevel_3d (output_terrain, 'WINDOWSIZE', '2.5 1000; 5 2000; 10 4000; 20 8000; 40 16000')
#process the .las file and add to the feature dataset and terrain

os.chdir (r'C:\Users\peter.a.metzger\Desktop\Projects\lidar_model_test_data\R06305328')
for LASfiles in glob.glob("*.las"):
    print LASfiles
    pointinfo = "point_file_information"
    #Get Pointfile information
    arcpy.PointFileInformation_3d(LASfiles, pointinfo, "LAS", "", projection)

    #get average point spacing variable
    pointspacing = arcpy.SearchCursor (pointinfo)
    field = "Pt_Spacing"
    for row in pointspacing:
        avg_space = (row.getValue(field))
    print avg_space

    #convert las  to multipoint

    #create lable for multipoint file
    listlable = list(LASfiles)
    del listlable[-4]
    multipointlable = ''.join(listlable)
    print multipointlable

    print created_FC
    outputPlable = os.path.join (created_FC, multipointlable)
    #class_code, change to a getparameter
    class_code = 2

    print outputPlable

    outputPlable_parameters = outputPlable + ' Shape Mass_Points 1 0 40 true false R06305328las_embed <None> false'
    print outputPlable_parameters
    arcpy.AddFeatureClassToTerrain_3d (output_terrain, outputPlable_parameters)

The error I'm getting is:

  "Traceback (most recent call last):
    File "C:\Users\peter.a.metzger\Desktop\Projects\LiDAR\las2terrain.py", line 64, in
      arcpy.AddFeatureClassToTerrain_3d (output_terrain, outputPlable_parameters)
    File "C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\ddd.py", line 1843, in AddFeatureClassToTerrain
      raise e
  ExecuteError: Failed to execute. Parameters are not valid.
  ERROR 000800: The value is not a member of .
  Failed to execute (AddFeatureClassToTerrain)."

And it's from the parameter that I set of: "outputPlable_parameters = outputPlable + ' Shape Mass_Points 1 0 40 true false R06305328las_embed  false'"

These values are what the model builder is using, yet it's not working.


texts: Getting QGIS "zoom to layer" to work with PostGIS data in multiple projections

texts:
So I have a .qgs file that contains multiple postgis layers, some in EPSG 4326 and some in EPSG 26910, with the .qgs file in 4326.  When I use the "zoom to layer" functionality on the layer in 4326, the view area zoom completely out to the entire projected globe and I get a sequence of 6 errors in the CRS Log Messages saying
Warning Simplify transform error caught: forward transform of (inf, inf)  Error: Fallback transform failed.
I attempted to alter the .qgs XML to find/replace instances of 'srid=4326' with 'srid=26910' but that only makes the opposite problem happen (layers in EPGS:26917 not working).  I have a feeling it may be related to https://github.com/qgis/QGIS/issues/34983, but not sure.
Does anyone know if there is a way to configure the .qgs file to work with layers in multiple projections, such that the "zoom to layer" functionality works across the different layers?


texts: Loaded shapefile into SDE without spatial proj creates problem

texts:
Via a Python script - various shapefiles were loaded into SDE with no spatial reference set.
This causes a strange problem - when you view the featureclass in ArcCatalog (geography view) - it shows the following dialog:

When you look at the Resolution tab of the featureclass it has an invalid resolution value:

It should normally show a small value - like as 0.000000001.  The resolution is effecively a mulitplier or divider that ESRI GDB storage uses to convert real numbers of the vertex coordinates into integers for storage, so in this context, ‘0’ is invalid as you can’t divide by zero.

Is there any way to change this value?, as it is certainly not editable in the dialog and
Calcualate Default XY Tolerance in ArcToolbox does not work.  Can it be altered by ArcObjects  - or could that break things.


texts: Finding out the extent of the default zoom levels of OpenLayers

texts:
I'd like to reproduce the zoom levels of OpenLayers in PostGIS. I am not used to working with JavaScript and am having difficulties to find out the exact extent in meters of the default zoom levels of OpenLayers. I can't find it anywhere on the internet.

With the function ST_Expand I can create a shape similar to a zoom level. I would like to use that function and recreate the default zoom levels 12, 13 and 14 of OpenLayers. How do I do that?

I am thinking of finding out the distance of the extent from OpenLayers, and using that as an input for the function ST_Expand. But how do I find that distance?


texts: Feature templates do not change when using LYR file and PostgreSQL?

texts:
We just started migrating from SQL Server to a PostgreSQL Server for production.  Traditionally, we load a LYR file containing all of our symbology, then point it to our new GDB on the SQL Server.  If we needed to modify an existing feature template, we could do so without any trouble.

Now, with PostgreSQL up and running, things are working a little different.  Following the same process, we load in the LYR file, point to our GDB, and start editing.  Now when we try to modify a feature template, specifically a field with a domain, the changes are not applied to features created from the template.  We can open the template, and our changes are visible, but new features created with the template do not reflect changes.

Also, if we set up a new map document, add data with the same symbology as the LYR file (manually, not using the LYR file), there are no troubles with the feature templates.

Are there any known issues/solutions related to my situation?

EDIT: Software versions

ArcGIS 10.4.1

PostgreSQL 9.4.21

SQL Server 2014

EDIT 2: Further testing

I've updated the PG cluster to 9.4.21, but the problem persists.  I've also tested using LYR files for other schemas we use, and there are no issues there.  I tried opening the LYR file in a fresh MXD, connecting to GDB, saving the MXD, and then running MXD Doctor to see if something was corrupted in the Featur Layer.  No luck there either.

EDIT 3: ESRI Support

ESRI Support has been chewing on the issue for 2 months now with no resolution.  They have been able to replicate the issue with ArcMap 10.4 - 10.6, and PostgreSQL 9.4.5 - 10.3.  They are convinced that this is a bug.  I'm surprised no one else has run into this issue before now.  After digging further into our production workflows, I've been able to reproduce the issue with different gdb scehams using MXDs and LYRs.


texts: Pyintaller cannot find pyqt4 plugin directories for a PyQGis program

texts:
I have a PyQt application that I am trying to compile into one exe using PyInstaller 3.1. It does compile into an exe, however, the GUI doesn't show up upon launch (clicking the exe). One of the compilation message is:

ImportError: No module named site
5576 ERROR: Cannot find PyQt4 plugin directories
5576 WARNING: Hidden import 'PyQT4._qt' not found (probably old hook)

I think this is why the UI wouldn't show up. I've already add all of the paths that has PyQt4 in the pathex attribute for the Analysis object to no avail. Below is my spec file:

# -*- mode: python -*-

block_cipher = None

a = Analysis(['main32.py'],
             pathex=['C:/OSGeo4W/apps/Python27/Lib/site-packages/PyQt4','./src'],
             binaries=None,
             datas=None,
             hiddenimports=[],
             hookspath=None,
             runtime_hooks=None,
             excludes=None,
             win_no_prefer_redirects=None,
             win_private_assemblies=None,
             cipher=block_cipher)
pyz = PYZ(a.pure, a.zipped_data,
             cipher=block_cipher)
exe = EXE(pyz,
          a.scripts,
          a.binaries,
          a.zipfiles,
          a.datas,
          name='MyExe',
          debug=False,
          strip=None,
          upx=True,
          console=False )

Any ideas what's next to try?


texts: AttributeError: 'QgsRasterLayer' object has no attribute 'pendingFields'

texts:
Remember i am using Python Version 3.6.2 QGIS Version 3.2

while running my code i am getting this error:

fields = selectedLayer.pendingFields()
AttributeError: 'QgsRasterLayer' object has no attribute 'pendingFields

if i am using fields() instead of pendingFields() because of qgis version 3x

getting this error:

fields = selectedLayer.fields()
AttributeError: 'QgsRasterLayer' object has no attribute 'fields

i am new in this field so please tell me properly what should i do

here is my code:

    filename = self.dlg.lineEdit.text()
    output_file = open(filename, 'w')

    selectedLayerIndex = self.dlg.comboBox.currentIndex()
    selectedLayer = layers[selectedLayerIndex]
    fields = selectedLayer.pendingFields()
    fieldnames = [field.name() for field in fields]

    for f in selectedLayer.getFeatures():
        line = ','.join(unicode(f[x]) for x in fieldnames) + '\n'
        unicode_line = line.encode('utf-8')
        output_file.write(unicode_line)
    output_file.close()


texts: Rendering Google Terrain in QGIS Layout Manager

texts:
I am using the Google Terrain map (via XYZ Tiles).

This seems to make the system very slow, and when making a .jpeg through the Layout Manager, the Google Terrain will not display (render) when zooming in.

I am using QGIS 3.4.2-Madeira. Can anyone suggest something to resolve this issue?


texts: Creating unique points at the start of multipart polylines in ArcGIS Pro

texts:
I have polylines from a network analysis representing public transport trips. These are multipart polylines.

I need to identify the starting point for each trip. However, using vertices to points to create starting points, I get starting points for every part of the line.
To get only one point, I tried transforming the lines with multipart to single part and then unsplit, as suggested here, as well as dissolve but none of them worked.

Fyi, the lines have IDs.


texts:
I am running QGIS 2.14 and PostGIS 2.3 for PostgreSQL 9.5.

Is there any way, whether in QGIS or PostGIS (or somewhere else), to automatically make polygons which are not contiguous to be contiguous, let's say by meeting halfway were they do not touch ?

I didn't find any consistent source of information on the web allowing this to be done, it even seems to me that the topology validation tools cannot consider the case as being a topological error, but I must say I didn't find these tools to be really ergonomic, so if it is possible I just didn't manage to make them work properly.

To illustrate, here comes a picture showing what my problem is:


texts:
Pan and Zoom toolbar greyed out on ArcGIS Desktop.

I recently installed ArcGIS 10.4.
I have checked the data frame properties and the Extent is set to Automatic.
I am using the data view not layout view.
Normally I could pan and zoom in and out with ArcGIS 10.3 on this view but now I can't select these options.

Can anyone advise me how to activate this toolbar again?


texts: How to import multiple textfiles in one step

texts: Is ist possible to import multiple textfiles (csv, tsv etc.) as layers in one step or do i have to do it file by file?


texts:
We have a web map application that uses a Geocortex viewer and references feature services on AGOL. Users get into the site with Active Directory credentials. When a user tries to update a record they receive the following error:

There was a workflow error running activity: Exception has been thrown by the target of an invocation.
Workflow 'Update Manhole Inspections Workflow' failed
Unhandled exception: 'One or more errors occurred.' in activity '1.396: Query for Feature'.
One or more errors occurred.
Error processing task results.
The request was aborted: Could not create SSL/TLS secure channel.

We've patched over our ArcGIS desktop installations with the TLS 1.2 patch but haven't published anything to AGOL since. Any ideas why this is happening?


texts: Creating a raster based on values of other rasters

texts:
I have three single band rasters, raster_a, raster_b and raster_c with pixel values 1 and 0.

I am trying to burn the values of raster_a to raster_b and create a "new_raster_b" where  if: raster_a == 1, raster_b becomes 1 and if raster_a ==0 raster_b becomes  0. Special condition: if raster_c==0 raster_b remains 0 even if raster_a ==1.

I am reading all the 3 rasters as gdal datasets and converting them to array. How can I apply the above mentioned conditions to create an array for new_raster_b using these arrays?

raster_a_ds = gdal.Open(raster_a_path)
raster_b_ds = gdal.Open(raster_b_path)
raster_c_ds = gdal.Open(raster_c_path)

print('Reading all rasters as arrays')
dcl_array = raster_a_ds.ReadAsArray()
tcm_array = raster_b_ds.ReadAsArray()
tcd_array = raster_c_ds.ReadAsArray()


texts: Checking via ArcPy if ArcMap is in edit session?

texts:
I've created a Python add-in button that helps speed up my coworkers workflow by copying one feature class attribute to another. It uses the arcpy.UpdateCursor function to update a row in the target feature class. As it exists now, this button script can be run regardless of the editing mode. Obviously when its run in an edit session, the user can choose to stop editing and not save changes, but this is not the case when the script runs outside of an edit session.

How can I add a check to the script that will stop the script from running if ArcMap is not currently in an edit session?

This concerns ArcMap 10 & 10.1

I also want to check with other ArcMap users to verify that updates to tables is not normally allowed without being in an edit session.

So how is this script running outside of an edit session?

This script also brings up another question about the seemingly serendipitous selection order ArcMap performs that just happens to work for me when I update the 2nd feature class table from a list, but that's for another day.

Here's the script as it works now (without any 10.1 editor implementation):

How to add a check to ensure the user is in an edit session?

def onClick(self):
    #Reference mxd
    mxd = arcpy.mapping.MapDocument("CURRENT")
    #Reference the main Data frame
    mm = arcpy.mapping.ListDataFrames(mxd, "MainMap")[0]
    #Reference the Water System Valve feature class
    waterValves = arcpy.mapping.ListLayers(mxd, "Water System Valve", mm)[0]
    #Reference the fire hydrant feature class
    fireHydrants = arcpy.mapping.ListLayers(mxd, "Water Hydrant", mm)[0]

    #Use the extent of the main DF to select all valves in the current view
    dfAsFeature = arcpy.Polygon(arcpy.Array([mm.extent.lowerLeft, mm.extent.lowerRight, mm.extent.upperRight, mm.extent.upperLeft]), mm.spatialReference)
    arcpy.SelectLayerByLocation_management(waterValves, "WITHIN", dfAsFeature,"", "NEW_SELECTION")

    arcpy.SelectLayerByAttribute_management(waterValves, "SUBSET_SELECTION", "LOCATIONID IS NULL")

    fields = ["LOCATIONID"]

    row, rows = None, None
    rows = arcpy.UpdateCursor(waterValves,fields)
    row = rows.next()
    valveList = []
    append = valveList.append

    #Loop through the valves table to update LocationID
    while row:
        builder = str(row.QSNO)+"-"+ str(row.VALVESEQNO)
        row.setValue("LOCATIONID", builder)
        append(builder)
        rows.updateRow(row)
        row = rows.next()

    del row, rows

    #New selection for fire hydrants
    arcpy.SelectLayerByLocation_management(fireHydrants, "WITHIN", dfAsFeature,"", "NEW_SELECTION")
    arcpy.SelectLayerByAttribute_management(fireHydrants, "SUBSET_SELECTION", "LOCATIONID IS NULL")

    row, rows = None, None
    rows = arcpy.UpdateCursor(fireHydrants,fields)
    row = rows.next()

    #Loop through fire hydrant table to update LocationID
    while row:
        for locID in valveList:
            construct = str(locID) + "-FH"
            #print construct
            row.setValue("LOCATIONID", construct)
            rows.updateRow(row)
            row = rows.next()

    del row, rows, valveList, mxd


texts: Can QGIS create Vector Polygons from Raster Data?

texts:
I normally use Esri, am new to QGIS, and don't know a lot of terminology and definitions.

I am working with georeferenced geologic maps. Each formation is coded with a different color or color pattern.

Is it possible for QGIS 2.18.12 to analyze these data, identify areas with that particular spectral signature, and then use the analysis to create vector polygons?

Hand digitizing is too time consuming.


texts: Create a polygon around a point with interpolated distances from point in ArcGIS using Python

texts:
I have a shapefile of points with 16 fields associated with each point - these four fields are Q1, Q2, Q3, and Q4 (quadrants). Each quadrant field contains a value that represents a distance from each point in nautical miles.

For example I have a lat/lon point with values:
Q1 = 10,
Q2 = 20,
Q3 = 30,
Q4 = 20,

So one point, with 4 rays (one in each quadrant) a Q# distance away from the point, and I want to calculate a round-ish polygon around the point based on these 4 distances.

I need to do this programmatically hopefully using arcpy, numpy, or some other python module.

I was thinking to use the degree-distance tool on arcmap, but I beleive that can only be done by "point and click" and not by an arcpy tool.


texts: Calculating overlap of several lines in QGIS

texts:
I've traced several lines that represent people's paths. Some of these lines overlap several times. I would like to get a view showing for each segment of my lines the number of line overlaps.

The methods I have used until now are:

convert the lines into points (via Qchainge plug-in) and then create a heatmaps.
using a procedure that shows the result in vector format (but very long):

Starting from the original layer "orig":

From "orig" create a layer of segments called "seg" by using Explode lines (QGIS).
From "seg" extract the center points of each segment, called "centr" using Point along lines with parameter 0.5 (QGIS).
from "centr" delete duplicates, call "centruniq" using Delete duplicate geometries (QGIS)
From "centruniq" create buffers of 50 cm called buff using Buffer (QGIS).
Run Count points analysis in polygon (QGIS) with input "buffer" and "centr", call output "count".
To "seg" element duplicates, call "seguniq" using Delete duplicate geometries (QGIS)
Join attributes from location of "seguniq" with count call final using intersects: segment vector with total number of passages attribute ("orig" elements)
style by category

I would like to ask if there was a faster solution that would allow me to display the number of line overlaps?


texts: Otsu thresholding and classic image processing algorithms on radiometrically calibrated Sentinel-1 SAR images

texts:
Is it possible to apply the Otsu thresholding (or any other image processing algorithm dedicated to image digital values in the interval [0, 255]) to a Radiometrically calibrated Sentinel-1 SAR product for which the pixel values (sigma0) are signed real (negative or positive decimal) and are not normalized (logarithmic scale)?

If it's feasible, how should the intensity values be rescaled to the byte interval ([0, 255]) prior to the processing?


texts: Selected Feature how to zoom in ArcMap viewer window

texts:
I will selected one feature in ArcMap that feature zoom in viewer window in ArcMap i tried in ArcMap VBA but no result

Public Sub zoomtool()

Dim pDoc As IMxDocument
Set pDoc = ThisDocument
Dim pMap As IMap
Set pMap = pDoc.FocusMap

-----------------------------------------------Here how to call viewer window

Set act = pMap

Dim pLayer As IFeatureLayer
Dim pFSel As IFeatureSelection
For i = 0 To pMap.LayerCount - 1

Set pLayer = pMap.Layer(i)
Set pFSel = pLayer

'Get the selected features
Dim pSelSet As ISelectionSet
Set pSelSet = pFSel.SelectionSet

Dim pEnumGeom As IEnumGeometry
Dim pEnumGeomBind As IEnumGeometryBind

Set pEnumGeom = New EnumFeatureGeometry
Set pEnumGeomBind = pEnumGeom
pEnumGeomBind.BindGeometrySource Nothing, pSelSet

Dim pGeomFactory As IGeometryFactory
Set pGeomFactory = New GeometryEnvironment

Dim pGeom As IGeometry
Set pGeom = pGeomFactory.CreateGeometryFromEnumerator(pEnumGeom)

pDoc.ActiveView.Extent = pGeom.Envelope
pDoc.ActiveView.Refresh
Next i

End Sub


texts: Using overviews of individual bands of a virtual raster

texts:
If I create a virtual multi-band raster (using gdalbuildvrt -separate to combine three layers into an RGB colour image) and each of the separate bands already has its own .ovr external overview, do I still need to create an overview for the RGB .vrt?
In other words, will applications (such as QGIS) perform just the same in either case? Or will most applications fail to utilise (and combine on the fly) the overviews of the individual bands?
If I do want to create a single RGB overview .vrt.ovr, can gdaladdo still use the overviews of individual bands (to avoid having to revisit the source pixels directly)? Is there an option I can specify to insist on whether it takes this computational shortcut or else regenerates them from scratch?


texts: Difference in lengh calculations for ESPG:32632 and ESPG:4326

texts:
I have a shapefile with CRS ESPG:32632 where are shown some power lines of a northern Italy area.

I need to calculate the length of the lines and add it, as a new field, in the attribute table. I am also asked to provide the shapefile with with CRS ESPG:4326.

If I save the file with the CRS ESPG:4326, reopen it and then perform the length calculation (function $length of qgis field calculator) the results appear incorrect.

If instead I perform the length calculation operation with the original CRS and only then save the file with CRS ESPG:4326 the result seems ok.

Can someone explain to me why?


texts: How can I export a Numpy array made from MODIS data for use in ArcMap?

texts:
I have an array made from MODIS Terra Daily snow data that I would like to export for use in ArcMap. How would I output the 'probablility_of_snow array' array near the bottom of this post to an image to use in ArcMap?

Here are some of the steps I followed to compute the array:

 from pyhdf.SD import SD, SDC
 import numpy as np
 import matplotlib.pyplot as plt
 import gdal

 #Get 2001 data
 file = SD(MOD10A1_A2001001_h09v04_006_2016091183527_HEGOUT.hdf], SDC.READ)
 datasets_dic = file.datasets()
 print ('datasets_dictonary\n', datasets_dic)
 print ('file info', file.info())

datasets_dictonary:

 {'NDSI_Snow_Cover': (('YDim:MOD_Grid_Snow_500m_16', 'XDim:MOD_Grid_Snow_500m_16'), (527, 1038), 21, 0), 'NDSI_Snow_Cover_Basic_QA': (('YDim:MOD_Grid_Snow_500m_16', 'XDim:MOD_Grid_Snow_500m_16'), (527, 1038), 21, 1), 'NDSI_Snow_Cover_Algorithm_Flags_QA': (('YDim:MOD_Grid_Snow_500m_16', 'XDim:MOD_Grid_Snow_500m_16'), (527, 1038), 21, 2), 'NDSI': (('YDim:MOD_Grid_Snow_500m_16', 'XDim:MOD_Grid_Snow_500m_16'), (527, 1038), 22, 3), 'Snow_Albedo_Daily_Tile': (('YDim:MOD_Grid_Snow_500m_16', 'XDim:MOD_Grid_Snow_500m_16'), (527, 1038), 21, 4), 'orbit_pnt': (('YDim:MOD_Grid_Snow_500m_16', 'XDim:MOD_Grid_Snow_500m_16'), (527, 1038), 20, 5), 'granule_pnt': (('YDim:MOD_Grid_Snow_500m_16', 'XDim:MOD_Grid_Snow_500m_16'), (527, 1038), 21, 6)}

File.info

 file info (7, 5)

I get the data for several years and I assign each year's data to an array:

 sds_obj = file.select('NDSI_Snow_Cover')
 Array_2001 = np.array(sds_obj.get())
 print (Array_2001)

The data looks like this for each year:

 [[ 84 85 80 ..., 250 250 250]
 [ 82 82 85 ..., 250 250 250]
 [ 78 82 82 ..., 250 250 250]
 ...,
 [250 250 10 ..., 250 250 250]
 [250 250 250 ..., 250 250 250]
 [250 250 250 ..., 250 250 250]

I do some data cleaning and math (not all shown here):

 probablility _of_snow = np.divide(snow_days,sample_size)

This is the array I want to export to make a map in ArcMap:

 probablility_of_snow array  :

  [[ 0.75        0.75        0.71428571 ...,  0.71428571  0.71428571
  0.71428571]
  [ 0.625       0.625       0.625      ...,  0.71428571  0.71428571
  0.71428571]
  ...,
  [ 0.42857143  0.5         0.5        ...,  0.71428571  0.71428571
  0.71428571]
  [ 0.42857143  0.5         0.375      ...,  0.71428571  0.71428571
  0.71428571]]

Here is some info from the MODIS file:

 #Get transformation and projection info
 width = 1038
 height = 527
 file = SD(MOD10A1_A2001001_h09v04_006_2016091183527_HEGOUT.hdf, SDC.READ)
 sds = gdal.Open(file, gdal.GA_ReadOnly).GetSubDatasets()
 geoInf = gdal.Open(sds[0][0])
 geoT = geoInf.GetGeoTransform()
 proj = geoInf.GetProjection()

 print(geoT)

 (-118.1513669444455, 0.0045054006101485754, 0.0, 44.71875, 0.0, -0.004505401117436153)

 print (proj)

 GEOGCS["Unknown datum based upon the Clarke 1866 ellipsoid",DATUM["Not specified (based on Clarke 1866 spheroid)",SPHEROID["Clarke 1866",6378206.4,294.9786982138982,AUTHORITY["EPSG","7008"]]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]

I don't understand what to do next to create and export an image to use in ArcMap.


texts: Dynamic property inside the wkt:// on Geoserver

texts:
I try to create dynamic symbolizers with WKT. This link says a CQL expression can be used inside a WellKnownName or a link:href. I want my symbolizers to have a dynamic size depending on my data on postgis.

Here's my SLD :

<?xml version="1.0" encoding="ISO-8859-1"?>
<StyledLayerDescriptor version="1.0.0" xmlns="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc"
  xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
  <NamedLayer>
    <Name>Real Size</Name>
    <UserStyle>
      <Title>Real Size</Title>
      <Abstract>Real Size</Abstract>
  <FeatureTypeStyle>
    <Rule>
      <Name>container</Name>
   <ogc:Filter>
     <ogc:PropertyIsEqualTo>
       <ogc:PropertyName>ct</ogc:PropertyName>
       <ogc:Literal>container</ogc:Literal>
     </ogc:PropertyIsEqualTo>
   </ogc:Filter>
<PointSymbolizer uom="http://www.opengeospatial.org/se/units/metre">
 <Graphic>
    <Mark>
      <WellKnownName>wkt://POLYGON((0 0,${length}/1.5 0,${length} -${width}/2,${length}/1.5 -${width},0 -${width},0 0))</WellKnownName>
      <Fill>
        <CssParameter name="fill">#0000ff</CssParameter>
      </Fill>
      <Stroke>
        <CssParameter name="stroke">#0000ff</CssParameter>
        <CssParameter name="stroke-width">100</CssParameter>
      </Stroke>
    </Mark>
    <Size>600</Size>
  </Graphic>

The error I've got is :

  Rendering process failed The specified mark "wkt://POLYGON((0
  0,${length}/1.5 0,${length} -${width}/2,${length}/1.5 -${width},0
  -${width},0 0))" was not found!

What am I doing wrong? Is there another way to pass dynamic attributes into my graphics? (like SVG?)


texts:
I am reposting a question that i asked in stackoverflow about leaflet and custom renderers.

Basically what i want to achieve is create a custom renderer that alters the shape of CircleMarker or extend the Path constructor to create custom shaped (vector) Markers.

What i've done so far is this:

import {
  CircleMarker as LeafletCircleMarker,
  Canvas as LeafletCanvas,
} from "leaflet";
import { withLeaflet, Path } from "react-leaflet";

const myRenderer = (layer) => LeafletCanvas({ padding: 0.5 });

class CustomMarker extends Path {
  createLeafletElement(props) {
    const enhancedProps = {
      ...props,
      renderer: myRenderer
    };

    const el = new LeafletCircleMarker(
      props.center,
      this.getOptions(enhancedProps)
    );
    this.contextValue = { ...props.leaflet, popupContainer: el };
    debugger;
    return el;
  }

  updateLeafletElement(fromProps, toProps) {
    if (toProps.center !== fromProps.center) {
      this.leafletElement.setLatLng(toProps.center);
    }

    if (toProps.radius !== fromProps.radius) {
      this.leafletElement.setRadius(toProps.radius);
    }
  }
}

export default withLeaflet(CustomMarker);

I am just setting the Canvas constructor as a renderer like the examples in the Leaflet docs but i am getting an Error: The provided object is not a Layer in layer/Layer.js:159 addLayer function.

Can someone please instruct me in how to implement a renderer?


texts: level of spatial aggregation of polygons

texts:
I am looking for a method that allows me to obtain the spatial aggregation of polygons for each cell of a grid with 10kmx10km resolution.

In each cell (10kmx10km), I want to be able to consider that big polygons that share a line segment as well as big polygons within certain distances (e.g., 1km), represent areas with high aggregation. Smaller polygons with the same spatial distribution represent a medium aggregation cell. While small disperse polygons represent areas with low aggregation.

What is the best way to do this (in ArcGIS, QGIS or R)?


texts:
I tried to convert tiff image to pdf using this command in OsgeoShell

gdal_translate -of PDF -a_srs EPSG:32643 bangalore.tif bangalore.pdf

It gets converted into pdf but the TerraGo plugin is not getting enabled for drawing polygons or lines in Adobe Reader.So what mistake i have done and How to do this ??


texts: Can style be saved with Qpackage and how can the new .qgs use the new .shp locations automatically?

texts:
I've used the Qpackage plugin to move all the layers of a project to a new location. At first my try was to move the source files, but it seems to have made new files, which is no problem.
But now it seems the .qgs on the new location doesn't automaticly use these new layer files. When I load in one of the files the style is also not the same as from the original file. Can I use Qpackage to save/move style files as well? And is there a way to make the new .qgs automaticly use the newly made files for the layers as well? Or do I have to change the source file for all the layers manualy.
I'm used to the packaging mechanics of Adobe Indesign, so the QGIS way of packaging isn't very logical to me.


texts: Using PostGIS relationships (M:N) in QGIS without loading all related tables in layertree

texts:
I am working with a PostgreSQL/PostGIS (version 12.3/3.0.1) database and I am using QGIS (version 3.10) to visualize and work with my data. In the database I have two schemas:

schema for spatial data tables, e. g. gebaeude
schema for documents, e. g. dokument, dok_gebaeude (without spatial information)

I set a m:n relation between gebaeude and dokument using the table dok_gebaeude.
My tables are:

dokument: id, bezeichnung, dateityp, pfad

dok_gebaeude: id, fk_dokument, fk_gebaeude --> fk_dokument references dokument(id) and fk_gebaeude references gebaeude(id)

gebaeude: id, geom, nummer, nutzung

In QGIS I loaded the relationship with the project properties>relation manager and I can use it to connect documents (stored in table dokument) with gebaeude-objects. All work fine. BUT...
In QGIS I only need the table gebaeude and I don´t want to load the tables dokument and dok_gebaeude in my layers (the latter are only technical tables). Unfortunately it is neccessary to load all three tables (dokument, gebaeude and dok_gebaeude) in my layers in QGIS, otherwise I am not able to add the relationship in my project properties>relations.
Is there any possibility to hide or remove or not even load tables (dokument and dok_gebaeude) in my layertree without losing the relationsships?


texts: Why is the following CSV file not being loaded as a points layer in QGIS?

texts:
I've been trying to add a certain CSV file as a points layer to QGIS for a while now, without much success. It can be viewed over here. (It was adapted from this .xls data file). In the online displayer, it does not render well. That's probably because information was lost when I saved the original .xlsx file as a CSV.
However, when I add it as a delimited text file in QGIS with ; as the regular expression delimiter, the sample looks good. Still, the point coordinates are not shown on the map when I choose longNieuw and latNieuw as the x and y coordinates, respectively.
I'm not sure what goes wrong, as I have earlier successfully added other point layers in QGIS.
Question: why can the linked CSV file above not be added properly as a points layer file in QGIS?
Note: here's an example of a coordinate tuple of the csv file I'm trying to upload:
28,69733300601 ; -100,51019685

The CRS of my own project is EPSG: 3857 . I've tried adding the points layer in both this CRS and in EPSG: 4326.


texts: How can I convert a DataTable to a GeoDatabaseTable

texts:
I have a DataTable full of "stops" and I want to convert it into a table that can be GeoCoded. However, I am getting an error when I try to execute it.

The Message: Error HRESULT E_FAIL has been returned from a call to a COM component.

The Error Code: -2147467259

public static void DataTableToGeoDataBaseTable(DataTable inTable, string outPath, string outName)
    {

        DataTable dTable = new DataTable();
        dTable = inTable;

        Geoprocessor GP = new Geoprocessor();

        // Create the tool process object.
        ESRI.ArcGIS.DataManagementTools.CopyRows copyrowsTool = new
            ESRI.ArcGIS.DataManagementTools.CopyRows();

        // Set parameter values.
        copyrowsTool.in_rows = dTable;
        copyrowsTool.out_table = outPath + "\\" + outName;    //"Z:\\Projects\\Project\\NewTable"

        // Execute the tool.
        GP.Execute(copyrowsTool, null);  //ERROR HRESULT E_FAIL HAPPENS HERE
    }

Where am I going wrong?

I have also tried this version

public static void DataTableToGeoDataBaseTable(DataTable inTable, string outPath, string outName)
    {

        DataTable dTable = new DataTable();
        dTable = inTable;

        ESRI.ArcGIS.ConversionTools.TableToTable tableToTable = new ESRI.ArcGIS.ConversionTools.TableToTable();
        tableToTable.in_rows = dTable;
        tableToTable.out_name = outName;
        tableToTable.out_path = outPath;

        Geoprocessor gp = new Geoprocessor();
        gp.Execute(tableToTable, null);

    }


texts:
I have installed the following before I installed ArcSDE 10.2.1:

ArcGIS for Desktop 10.2.1(standard)
ArcGIS for server Enterprise 10.2.1(standard)
MS SQL server 2008
create enterprise geodatabase successfully

After I installed ArcSDE 10.2.1, when I tried to perform arcsde post installation" I cannot find it in "Start -> All Programs -> ArcSDE -> arcsde for microsoft sql server post installation", I only got "Command References" under ArcSDE.

Anything I missed?


texts: Representing Roads as Corridors using ArcGIS Desktop?

texts:
I am building a map of historic roads with ArcGIS Desktop 10.7.1.  Because of the poor accuracy of old maps and the fact that historic roads changed alignment often - sometimes seasonally - a better symbolic representation of such a road is to show a central alignment as a line, but then to show buffer on either side of that central line which together would represents a "corridor".  The width of the corridor can vary by road segment depending upon the accuracy [or not] of the source data.

I believe I have seen this done before, but I cannot find a reference to it, probably because I'm using the wrong terminology.


texts: What is the maximum number of characters for a label in QGIS 3?

texts:
I am trying to determine the maximum number of characters for a label in QGIS 3. I am at over 1,000 characters and it does not display. Is there another way to display text on a map with a large number of characters - text box, etc.?

City    Latitude    Longitude

Some City   51.500210   -0.126280
Second City 48.852421   2.351008
Third City  48.211534   16.384351
Fourth City Lorem ipsum dolor sit amet, consectetuer adipiscing elit. Aenean commodo ligula eget dolor. Aenean massa. Cum sociis natoque penatibus et magnis dis parturient montes, nascetur ridiculus mus. Donec quam felis, ultricies nec, pellentesque eu, pretium quis, sem. Nulla consequat massa quis enim. Donec pede justo, fringilla vel, aliquet nec, vulputate eget, arcu. In enim justo, rhoncus ut, imperdiet a, venenatis vitae, justo. Nullam dictum felis eu pede mollis pretium. Integer tincidunt. Cras dapibus. Vivamus elementum semper nisi. Aenean vulputate eleifend tellus. Aenean leo ligula, porttitor eu, consequat vitae, eleifend ac, enim. Aliquam lorem ante, dapibus in, viverra quis, feugiat a, tellus. Phasellus viverra nulla ut metus varius laoreet. Quisque rutrum. Aenean imperdiet. Etiam ultricies nisi vel augue. Curabitur ullamcorper ultricies nisi. Nam eget dui. Etiam rhoncus. Maecenas tempus, tellus eget condimentum rhoncus, sem quam semper libero, sit amet adipiscing sem neque sed ipsum. N    45.438171   12.334102
Florence    43.769019   11.250152


texts: Gdal ReadAsArray() different values than gdal band statistics

texts: I am confused about why if I open a band with GDAL I get, for example, band.GetMaximun(): 1000, but if I do ReadAsArray I get myArray.max(): 1200 (the same if I do np.max(myArray)). Does someone knows why this difference?


texts: QGIS equivalent to "Features to Line" ArcGIS tool

texts: I have a python script in Arc10 that uses the features to line tool. I however don't have the license to run it. I was wondering if there is a Qgis tool or process that would essentially do the same thing?


texts: How to split two interscting polygons by the average linestring of the shared area (QGIS / python shapely)

texts:
I have two polygons,

  POLYGON((-5 -5,-5 5,5 5,5 -5,-5 -5))
  POLYGON((1 1, 1 7, 7 7, 7 1, 1 1))

and wish to split them evenly.
The result should be a split along the red marked line.

I'm using qgis/ python shapely and wish to do it using python.