texts: Creating grid index (fishnet) features with indexing according to columns

texts:
I'm working with west to east river features and I need to index a grid using only numbers and accordinly to columns and not lines. The grid index features tool gives me this numerical sequence of index, for a grid with 3 lines and 6 columns:

1;2;3;4;5;6;

7;8;9;10;11;12;

13;14;15;16;17;18

But I really need this one:

1;4;7;10;13;16

2;5;8;11;14;17

3;6;9;12;15;18

Any tips?


texts:
I'm having a problem importing a .pbf file in PostgreSQL using the osm2pgsql.
I've downloaded the italy-latest.osm.pbf and as I do every year I'm going to upload it in PostgreSQL.
I did it for other country files and all worked well, but when I try with this one the process starts and after 4 hours after having read all nodes, ways and relations it simply stops.

I tried with the -C option augmenting the cache size, but nothing change.
Now they are 3 days I'm trying but I'm not an expert in these kinds of things.

Do you know if there is a more frindly way to import these data in Windows?


texts: OpenLayers Show message (Information message or Alert )(or HTML content)on the map

texts:
J'utilise OpenLayers et j'aiemrais afficher une fenetre affichant des messages sur la carte openLayers

Voici un exemple de ce que je voudrais faire : http://www.annuaires.justice.gouv.fr/annuaires-12162/annuaire-des-juridictions-21798.html

j'aimerais afficher un cadre comme la fenetre où il y écrit:
"Affiner votre recherche..."

J'ai 'essayé de faire :

Test

et la j'ai bien une petite fenêtre qui s’affiche sur la carte mais y a un décalage en haut de la carte.

savez vous comment afficher du texte ou du contenu html sur la carte ?

English Translation (based on Google Translate. Improvements welcome!)

I'm using OpenLayers and I would like to display a window that displays messages on the OpenLayers map.

Here is an example of what I would do: http://www.annuaires.justice.gouv.fr/annuaires-12162/annuaire-des-juridictions-21798.html

I would like to display a frame as the window where it's written:
"Affiner votre recherche ..."

I'm trying to do:

Test

and I have a small window that appears on the map but there is a gap at the top of the map.

Do you know how to display text or html content on the map ?


texts: Complex combination of raster bands in Geoserver

texts:
I have a TIFF file containing various bands. At the moment I am using SLD to assign red/green/blue channels to bands. I would like to be able to perform more complex combinations on those bands, eg add 2 bands and assign it to the RED channel for instance.

Is this currently possible with Geoserver?


texts: Finding the closest point to a polygon in a shapefile using GeoPandas

texts:
I have the following map of statistical regions of Belgium:

The blue dots correspond to coordinates for which I have temperature data. Surrounded by the red circle is Brussels. As it can be seen, there are no points inside Brussels and therefore I don't have temperature info for it.
I would like to assign the temperature in the closest blue point to Brussels.
How could I find the closest point of the grid to the Brussels polygon using GeoPandas?


texts: Editing attribute values in multiple fields using ArcMap

texts:
I have about 500 shapefiles that contain polygons. One polygon in each shapefile. I want to make a new field in their attribute tables called name.
To give you an example the name of each shp has the same formula like "Ikaria_20_09_16". Ikaria is the place and the rest are date. Each shapefile has different name of place and of course different date.
How can I perform an edit where I can take the first word of the name and put it in the field?
In addition I want to make new fields about date but in separate ways. One for day, one for month and one for year. I want also to take the from the name of the shp what I need and put it in field.
Is there any way to do that?


texts: Bounding box generation, retaining fields?

texts: I generated a bounding box by calculating the extent of each feature using "Polygon from layer extent..." tool located in "Vector/Research Tools" which works great but I lost the fields, especially the name of each polygon. Does anyone know how to keep the fields in the process?


texts: Changing default color palette from RGB to CMYK in ArcMap?

texts:
Is there a way to change the default colour palette in ArcMap? Currently, for each symbol, I have to change RGB to CMYK to choose my colour (see screenshot). This is pretty annoying since all my colours are in CMYK.

I also had a look at the utility tool AdvancedArcMapSettings.exe (more details in ArcMap Help) but it doesn't provide an option.


texts: Creating empty copy of active layer in PyQGIS?

texts:
I'm trying to create an empty copy of a certain layer, so no objects only tablestructure, crs etc.

from qgis.utils import *
layer = qgis.utils.iface.activeLayer()
print layer.name()
provider = layer.dataProvider()
print provider.fields()[0].name()
writer1 = QgsVectorFileWriter("d:/temp/qgis/test3.shp",
                              "CP1250",
                              provider.fields(),
                              provider.geometryType(),
                              provider.crs(),
                              "ESRI Shapefile")

When I open the created file (test3) there is no tablestructure, crs. etc.

Running the script the second time with test3 still open, it works.

Obviously I miss something .

An other good working solution I found later:
http://gis.stackexchange.com/q/156096/67332/

it loads it into a memory layer.


texts:
I've a huge number of shapefiles and I would like to create a VRT as I do with gdalbuilvrt for raster.
I find here https://github.com/dwtkns/gdal-cheat-sheet that it should be possible and the format of VRT for OGR vector files exists and is described here : http://www.gdal.org/drv_vrt.html.

Where can I find documentation or a utility program to do that?


texts: Why does QGIS 2.2.0 displaying triangles instead of rectangles?

texts:
When displaying the demographic datas from INSEE (http://www.insee.fr/fr/themes/detail.asp?reg_id=0&ref_id=donnees-carroyees&page=donnees-detaillees/donnees-carroyees/donnees-carroyees-200m.htm), the results under qgis 2.2.0 is quite surprising :

instead of rectangles I just see triangles in the bottom right corner (so covering half the surfaces)
the same shape file opened with MapInfo 10.0 seems correct with all the rectangles correctly displayed and filled.

I can't see the reason of such a difference, it seems for me to be a bug of Qgis ..


texts: I have downloaded Landsat collection 1 Level 1 imagery (Landsat 4,5,7 and 8) from https://earthexplorer.usgs.gov/. I would like to know that do we have to apply some conversions or corrections before calculating any index for example ndwi and ndvi or is this analysis ready data?


texts: GeoServer limits number of SQL View parameters?

texts:
I tried to create a sql view with 53 params, but when I clicked at "Guess parameters from SQL" it showed only 25 params.

If you have had this issue, is there a workaround?

The featuretype.xml file contains all 53 parameters, the problem seems to be on the sql view interface.


texts: Identify polygon touching common polygon

texts:
In the case of the image towards left, the red polygon is overlapping the blue polygons.

In the case of the image towards right, two red polygons are overlapping blue polygons. The two red polygons are overlapping on one common polygon, i.e. polygon no. 2.

Is there any way to identify the polygons, i.e. particularly those polygons in which 2 or more red polygons are overlapping on a common blue polygon?


texts: Why is the OGC Symbology Encoding standard not being implemented and is SLD an alternative?

texts:
I searched exhaustively over the internet as well as at this Exporting Symbology from ArcGIS for Desktop to Microstation or QGIS? but didn't find a satisfactory answer. Hope to get the answers from the experts over here.

Suppose I have styled a polygon shapefile (e.g.: Border width = 0.26, border color = green, red colored simple fill, label displayed, etc ). When I save the style in QGIS, it is saved in XML format which is specific only to QGIS. Same is the case with other software which they store it in their own format. Which means, I cannot exchange styling information between different software.

Though the OGC SE standard is in place, I was just confused why this is not being implemented by any software. QGIS 1.8 supports SLD but then what is the purpose of storing the style as QGIS specific .qml file ?

Also, since SLD standard is used for web-services, is it OK to use this standard in this case (for exchanging data between different software)?


texts: Joining field works in ArcMap but not in ArcPy?

texts:
I am trying to join a field of values found within a DBF to a multi-polyline shapefile. They share a common field name (filename) and the values found within each are identical. I was able to run Join Field without a hitch in ArcMap, but I'm receiving two errors when I try to run it in Python with identical parameters.  Here's the block of code:

merge_features = r"C:/Users/kellyj/Desktop/Projects/CostAnalysis/merge_features"
import arcpy
from arcpy import env

arcpy.env.workspace = merge_features
arcpy.AddJoin_management("merged_lines.shp", "filename", "merged_table.dbf", "filename")

These are the errors I receive:

  ExecuteError: Failed to execute. Parameters are not valid.
  The value cannot be a feature class
  ERROR 000840: The value is not a Raster Layer.
  ERROR 000840: The value is not a Raster Catalog Layer.
  ERROR 000840: The value is not a Mosaic Layer.
  WARNING 000970: The join field filename in the join table merged_lines is not indexed. To improve performance, we recommend that an index be created for the join field in the join table.

Any thoughts as to what's happening here?


texts: Is it possible to add a command to the ArcMap Find dialog context menu?

texts:
I think I already know the answer to this... which is "no", but I wrote this up so I'll ask anyways.

Is it possible using ArcObjects to extend/customize the ArcMap Find tool/dialog by adding a command to the context menu that appears when you right-click a search result? The command would also need to know the location of the item that was right-clicked.

I see that there is an IFinder interface one can implement to create a new find dialog but I just want to add a menu item to the existing one.

My goal is to provide an address lookup using a user-specified Geocoding service (i.e. what the "Locations" tab does), and directly run a custom identify command at the location of the user's chosen search result.

If I can't use the existing functionality to do this then it's not worth reinventing the wheel as it would just eliminate an extra step or two for the user. I'm also doing this within an ArcMap add-in, so implementing interfaces, registering component categories, etc. is not an option.


texts: Display Leaflet WMS Layers from same source on different panes

texts:
I am using the Leaflet WMS plugin to add in layers to my map via a WMS.

I have the following code below that adds all of the layers using the same source but as separate layers.

However, I want to set each of the layers to its own pane so that when I add/remove them to/from the map they appear in the correct order as defined by the list of panes, rather than the order they were added to the map as currently happens without the panes.

var source = L.WMS.source("http://IP ADDRESS:PORT/geoserver/wms", {
format: 'image/png',
transparent: true,
tiled: true,
opacity: 0.6,
minZoom: 12,
tileSize: 512,
attribution: "Attribution Here"
});
var layer1 = source.getLayer("WORKSPACE:layer_1");
var layer2 = source.getLayer("WORKSPACE:layer_2");
var layer3 = source.getLayer("WORKSPACE:layer_3");
var layer4 = source.getLayer("WORKSPACE:layer_4");


texts: How to check how many line intersects a intersecting point in same geometry

texts:
I want to know that how many line intersects a point. In figure below it can be seen that 3 lines crosses a point. 2 red lines and 1 blue line. I want to know in the same geometry that many lines crosses an intersecting point.

I have tried it with st_intersection but its not working. Suppose if I have a point similar to that in the picture the query should return 3.

Is there any way in to check how many lines crosses an intersection point?

what I have tried is:

select Count(Distinct a2.id) as no_of_line_intersecting,
st_intersection(a1.geometry,a2.geometry)
 as geom from Area a1,Area a1 where
ST_Touches(a1.geoemtry, a2.geometry)

But I am not getting the right results. Does somebody have any idea about this? I am using spatialite


texts: Why does pyQGIS QgsVectorLayer().isValid() return FALSE with Oracle views?

texts:
This works perfectly fine with an Oracle TABLE but not with a VIEW. Why?

uri = QgsDataSourceURI()
uri.setConnection(host, port, db, username, pw)
uri.setWkbType(QGis.WKBPolygon)
uri.setDataSource(schema, table, geomcol)
vlayer = QgsVectorLayer(uri.uri(), table, 'oracle')
print vlayer.isValid()

..and returns TRUE. However, the same code does not work with a VIEW. Both TABLE and VIEW have correct entry in USER_SDO_GEOM_METADATA. All spatial tables in the database load correctly; none of the views do.
All spatial tables and views load correctly in QGIS using Layer->Add Layer->Add Oracle Spatial Layer. Is it meant to work only with tables or am I doing something wrong?


texts: Is there a tool in QGIS similar to ArcGIS Combine?

texts:
I am looking for a tool in QGIS 2.18 similar to the Combine tool in ArcGIS for rasters. I have four different raster layers that I would like to combine into a single raster that has values that represent the unique combinations of the source rasters. My source raster layers have 2-6 different classes. E.g. geology types 1-6, topographic position class 1-3, Aspect class 1-2. So they have values that range from 1-6. I would like to use this to stratify my study area by these unique combinations so that I can then randomly select sample sites.

This other question (Seeking tool in QGIS similar to ArcGIS Combine?) finds a solution that tells you how much of the area has different combinations, but I would like the result to be either a raster or vector layer so that I can then use it for other spatial analyses.

I can probably do it by converting them to vector layers by doing some unions and then some field calculations. But I didn't want to have to stuff around with invalid geometries that always seem to crop up when working with vectors and this seems a clunky way to go about it.


texts: CSRS to NAD83 datum shift with ogr2ogr using grid shift file fails

texts:
I am trying to do a datum shift transformation using ogr2ogr on the Ubuntu Desktop 16.04 command line. All files, including the grid shift file, are in the same directory. I run the following command:

  ogr2ogr -f "ESRI Shapefile" -s_srs "+proj=utm +zone=17 +ellps=GRS80
  +towgs84=0,0,0,0,0,0,0 +units=m +nadgrids=ON83CSv1.gsb" -t_srs EPSG:26917 survey_marker_NAD83.shp
  survey_marker_csrs.shp

I get the following error:

  ERROR 1: failed to load datum shift file
  Failed to reproject feature 0 (geometry probably out of source or destination SRS).
  ERROR 1: Terminating translation prematurely after failed
  translation of layer survey_marker_csrs (use -skipfailures to skip errors)

What am I doing wrong?


texts: MGRS(Military Grid Reference System) reverse geocode in ArcMap?

texts:
How can we reverse geocode MGRS (military grid reference system) strings in ArcMap?

MGRS strings use alphanumeric notation to identify a location on the earth.

For example, I have 34TCK9320739642. How can I get the address (city, state, name, country, etc.)?


texts:
I'm trying to calculate values for a raster I have based on other raster values.

The problem is that I get the error:

  -0.054.multiply is not a function

and I believe it is because I'm writing something wrong, but can't find the problem.

This is how I try to calculate the new values for the raster:

var NewRas=(2.12).subtract((-0.054).multiply(RAS1))+((0.12).multiply(RAS2)).add((3.14).multiply(RAS3));

My end goal is to calculate the new pixel values for my new raster base on this equation:

EDIT: I have managed to calculate it this way, as demonstrated in GEE tutorial:

But still I don't understand why what I wrote is wrong


texts:
I want to import GeoJSON data into MySQL to compare performance with other databases. Is there any way I can import following data into MySQL. And if it possible, how to create table for such records.

{"type":"FeatureCollection","features":[{"type":"Feature","geometry":
{"type":"Point","coordinates":[-122.085823,37.386337]},"properties":
{"id":"1","name":"Frankie Johnnie & Luigo Too","address":"939 W El Camino Real, Mountain View, CA"}},
{"type":"Feature","geometry":{"type":"Point","coordinates":[-122.083237,37.387138]},
"properties":{"id":"2","name":"Amici's East Coast Pizzeria","address":"790 Castro St, Mountain View, CA"}},
{"type":"Feature","geometry":{},"properties":{"id":""}},
{"type":"Feature","geometry":{},"properties":{"id":""}}]}


texts: Closing QGIS custom form programmatically?

texts:
I have a custom form for entering and editing attribute values of a point shapefile. I also added a button for deleting the actual feature. Deletion of the feature works fine and then I want to close the form.

I tried:

MyDialog.close()

but it only empties the dialog (all widgets gone) but leaves an empty window.

How can I close the custom form programmatically?

System is Windows 8.1, QGSI 2.8.1


texts: Can you use LEN in SelectLayerByAttribute_managment in Arcpy?

texts:
I am trying to use LEN to select all strings in a Column that are a length of a phone number, and am unable to enter the correct expression, is it possible to use LEN in SelectLayerByAttribute_management ?

arcpy.SelectLayerByAttribute_management("Contacts_tbl", "NEW_SELECTION", "LEN(Contact_Phone_Work) = 12")


texts:
I have an issue with creating an address locator on ArcGIS 10.2.2.
I always used an address locator somebody else made, it's made on ArcGIS 9.2. When I look into its properties, it says

  Description: US Streets, Style: N/A

.
The thing is I made a locator by my own, using the US address dual range style, which I think it's the same as the previous locator. I use all the required fields.

But the new locator takes forever to complete the job, and most of the times it ends up crashing, whereas the old locator keeps working smoothly.

Anybody has a clue on why this could be happening?

SOLVED:

I ran into this explanation:
http://blogs.esri.com/esri/supportcenter/2010/10/29/haunted-by-geocoding-brew-up-something-new-with-9-3-1-locator-styles-in-10/

So I downloaded the locator styles from 9.3 and followed the instructions to add them.
Then I created a 9.3 style locator, and it worked as a charm. I guess it's a 10.2.2 bug.


texts: Repaint the map after updating a feature in QGIS

texts:
I update a feature in a vector layer using the following Python code:

layers = QgsMapLayerRegistry.instance().mapLayers()
layer = layers.get('my_vector20140614125509945')

fid = 1
attrs = { 0 : "hello", 1 : 100 }
layer.dataProvider().changeAttributeValues({ fid : attrs })

The feature is updated but the map is not redrawn.  I am using catagorized style.  The value of 100 will be shown with a different color.


texts: Is it possible to create print templates in 2.6 and 2.8?

texts:
Is this possible? I can only find the option to 'Add items from template'. When using this option, the page setup defaults to A4 landscape and you must manually change these settings using the composition tab within the print composer.

Version 2.4 and earlier did not need this extra step. Applying a pre-saved print composer not only loaded in all the objects previously saved within the template but also replicated the page set-up

Am I missing something really obvious?


texts: How to implement a shoelace formula into Python 2.7.8

texts: Here's the image of the equation I'm using, it is a shoelace formula. I need help to point out what each variable represents and how to apply this in python script, math is not my strong-suit! From what I understand the X and Y are coordinate inputs and "A" is Area. I don't know what "i" is. The end goal is to create centroids.


texts: QGIS 3.19 Development version installment

texts:
I'm trying to install the qgis development version 3.19, as I need it to load a map from Maptiler.
While trying to install I keep getting the same error:
I do get the 3.19 on my device, but when opening i get this error:

----------------------------------------
Does anyone know how to fix this or what I'm doing wrong?


texts: How to let ogr/gdal CreateLayer() create a `geom` field instead of `wkb_geometry`?

texts:
When creating a new table (PostGIS) using OGR/GDAL (2.1+), is there a way to select a specific name for the geometry field (such as geom), rather than using the default name?

I created a new table via python3 ogr/grdal liek this:

ds = gdal.OpenEx(connection_string, gdal.OF_VECTOR | gdal.OF_UPDATE)
lyr = ds.CreateLayer( "mylayer", XXX, ogr.wkbLineString)

The problem is that I can't find a parameter to specify the name of the geometry field, only its type. In the table created, the geometry field seems to be called wkb_geometry by default:

wkb_geometry geometry(LineString,26945)

Is there a parameter or option to let CreateLayer() use a different name such as geom?

Related question: Use Postgis to convert wkb_geometry data type to geom datatype


texts: Setting extent while interpolating in ArcGIS Desktop?

texts:
I have a set of lines representing level curves. I have also a polygon.

I want to interpolate them for creating a hypsometric map. I use the Topo to Raster command but the problem is that even if I set the extent of the map with that polygon layer, my map is created bigger. I tried to clip it but it remains the same.

I want that my map to be only as that polygon.

Is there a clipping method or a property that I can set while I interpolate the curves in order to generate a map with that extent?


texts: I am planning to do an advanced program in spatial science and engineering. It states as a requirement in the program, that the student should have a through understanding in math and statistics. So far I have done only up to undergraduate level (i.e. only A/Ls) in both the fields and am baffled as to how much of understanding I should have in the fields. Also I my BS degree is a special in Geology. As the program is basically an engineering graduate level one, what do you recommend? Some good books, would be highly appreciated.


texts: Creating Custom Layers using ArcGIS Pro SDK for .NET?

texts:
We have a Custom Extension developed for ArcMap, which includes the Creation of a Custom Layer. The Code for the Layer is similar to what you see in this Answer to Custom Layer not drawing using ArcObjects?. This has been developed in ArcObjects & C#.

I'm now evaluating migrating the Extension to ArcGIS Pro, but can't find any documentation for a Custom Layer in ArcGIS Pro. The Only thing I find is an ArcGIS Idea for Custom Layer in ArcGIS Pro.

Can we create Custom Layers using the ArcGIS Pro SDK for .NET ?


texts: Difference between Geostatistical Analyst and Spatial Analyst toolboxes for interpolation?

texts: I did interpolation on point data using "Geostatistical Analyst" tool. It gives me surface layer. Then I exported this layer into raster image (by right click on this layer=>export data into raster) for convenience of raster calculation. But few minutes ago I found another approach on internet to create raster image from point data. That is to use "Spatial Analyst tool". So I'd like to ask there is any difference with output of further calculation between the two raster images created from 2 difference interpolation approach in ArcGIS for Desktop.


texts: I am trying to use the tobler's hiking function to estimate the distance/time from an origin to destination points. The contour interval is 500 but I do not know how to interpret the results. Does it mean that the distance between each of the contour line is 1 hour because I used aniso_hr1 and where does the contour interval come in the interpretation?


texts: Storing large number of mapsets with many raster maps

texts:
OS X El Capitan, GRASS GIS 7.0.1

I have a location with about 530 read-only mapsets with each about 2600 raster maps (plus PERMANENT and one link to a read-write mapset), overall 11.255.481 raster maps and about 181 GB.

What is the best way of storing these in a space efficient way and to
have them accessible read-only?

I have them at the moment in a grass
database and thought of putting these into a squashfs which should
work (tried with a subset and it was nice - mounting fast, space efficient), but after two days of compressing (and only having compressed
30% - slow speed presumably due to huge number of files) I gave up.

I would like to avoid having to install postgresql as I tried it before
(several years ago) and given up due to it's complexity.

I will analyse the maps in R, so any storage backend which I can use
From R and GRASS would be OK.

Any suggestions?

ADDED INFO

Number of Files   : 11.604.844  (find . -type f | wc -l)

Size of all Files : 181G   (du -sch ./*)

One timing example:

01:39:24 /Volumes/simASM$ time mksquashfs ./simASM.Outeniqua.Random.Outeniqua_2_0.Outeniqua.41719.* ~/tmp/simASM.squashfs
Parallel mksquashfs: Using 8 processors
Creating 4.0 filesystem on /Users/rainerkrug/tmp/simASM.squashfs, block size 131072.
[=======================================================================================================================================================================================|] 92263/92263 100%

Exportable Squashfs 4.0 filesystem, gzip compressed, data block size 131072
    compressed data, compressed metadata, compressed fragments, no xattrs
    duplicates are removed
Filesystem size 365271.48 Kbytes (356.71 Mbytes)
    50.25% of uncompressed filesystem size (726846.67 Kbytes)
Inode table size 801639 bytes (782.85 Kbytes)
    20.23% of uncompressed inode table size (3962401 bytes)
Directory table size 617819 bytes (603.34 Kbytes)
    20.52% of uncompressed directory table size (3010422 bytes)
Number of duplicate files found 84524
Number of inodes 122969
Number of files 110768
Number of fragments 1775
Number of symbolic links  0
Number of device nodes 0
Number of fifo nodes 0
Number of socket nodes 0
Number of directories 12201
Number of ids (unique uids + gids) 2
Number of uids 1
    rainerkrug (501)
Number of gids 1
    staff (20)

real    2m2.877s
user    1m6.972s
sys 1m57.861s
01:41:31 /Volumes/simASM$


texts: Calculate average number of direct neighbors in cities in PostGIS

texts:
In need to calculate the average number of direct neighbors in pre-selected cities in PostGIS.
So far, I have only found solutions to select neighbors based on distance (e.g. all neighbors within 1km from my focal point).
What I need are really the direct neighbors (so no other points in between) irrespective of the distance.
I have a table where I can select coordinates by cities, e.g. select lat, lng, geom from mytable where city = 'Munich' (I have also created a geom column with gist index).
Ideally, I want to calculate for each row in the table (selected by cities) the number of direct neighbors that are also located within the same city.


texts: pgRouting not returning node with correct z-value

texts:
have an issue with pgRouting that I've modified to consider multi-floor indoor areas. I've created a mock floor plan of any typical hotel. I copied the floor plan and assigned a different z-value to essentially mirror the floor plan, with one positioned directly above the other.  I then made vertical connections between the floors for elevators and stairs to create a complete 3D network.

Using the SQL code found here https://github.com/mapcentia/geocloud2/blob/master/app/scripts/sql/pgr_createTopology3D.sql, I created the modified pgr_createTopology3d and pgr_pointToIdZ functions to build topology on the 3D network. Following that, I followed normal procedure for building the SQL View to serve the queried results as a layer to GeoServer

When constructing the SQL View, I used a modified version of pgr_fromAtoB to consider z-values.  I modified as follows (see @underdark's answer to pgRouting with multiple floors):

-- Find nearest node
EXECUTE 'SELECT id::integer FROM ways_vertices_pgr
    ORDER BY the_geom <-> ST_GeometryFromText(''POINT('
    || x1 || ' ' || y1 || ' ' || z1 || ')'',4326) LIMIT 1' INTO rec;
source := rec.id;

EXECUTE 'SELECT id::integer FROM ways_vertices_pgr
    ORDER BY the_geom <-> ST_GeometryFromText(''POINT('
    || x2 || ' ' || y2 || ' ' || z2 || ')'',4326) LIMIT 1' INTO rec;
target := rec.id;

However, it doesn't appear that the z-value is being considered.  For example, I have nodes lying exactly on top of each other with the same lat/long, but with differing z-values. I can issue the ST_GeometryFromText command and limit the results to 2 and see the both of the nodes are being returned at the lat/long I've specified.  However, the z-value for the bottom node is 0 and the upper node is 1.  No matter what z-value I specify in the request, the upper node is always returned as the top result.

Does anyone know why?

How can I get it to return the bottom node if I specify the z-value as 0?


texts: How to rotate (svg) symbols in QGIS map composer to match map rotation?

texts:
Because of my thesis layout requirements, I have to rotate 270° some of my maps in composer.
There are (many) svg smbols (icons), and they are rotated as well.
I would like to (non-)rotate them in order to keep them upright, so one doesn't have to continuosly rotate the page or the neck ;)
Is there a way to do that? Note that I would like to do it on the composer stage since some of my maps are normally placed north-south.
Thanks for the attention!


texts: Formatting h-axis of time series chart in Google Earth Engine

texts:
I am trying to make a time series NDVI based on mm-yy. Since I get monthly images to a single image (in map_m function), I lose 'system:time_start' from properties. So, I use the loop's time to set a new 'system:time_start' property. But when I create the time series chart using this property, I cannot get time mm-yy format in my chart. I guess that's because GEE is expecting original format of 'system:time_start' value of sentinel images.
How to get 'mm-yy' in h-axis of chart using the 'system:time_start' property I created?

Link to GEE code:

Here is my code:

// Point for filterbound
var roi = ee.Geometry.Point([92.15091245970632,21.205549008619254])

// Center the display to ROI.
Map.centerObject(roi, 11);

// Function to get NDVI.
var getNDVI = function(image){
var NDVI = image.normalizedDifference(['B8', 'B4']).rename('NDVI')
var clouds = image.select('QA60').gte(1024).unmask().not().rename('CloudMask')
return image.addBands(NDVI)
};

// Gets the current time.
var eeNow = ee.Date(Date.now());

// Gets image collection of Observation Period
var series = ee.ImageCollection('COPERNICUS/S2')
          .filterBounds(roi)
          .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 10))
          .map(getNDVI)
          .select(['NDVI']);

// List of years and months
var years = ee.List
 .sequence(2016, ee.Number.parse(eeNow.format('Y'))); // Gets latest year
          // from eeNow as string then converts string to a number
var months = ee.List.sequence(1, 12);

// Filters images annualy-monthly
var map_m = function(i) {
i = ee.Number(i);
var filtered_col = years.map(function(j) {
var filtered = series
    .filter(ee.Filter.calendarRange(i, i, 'month'))
    .filter(ee.Filter.calendarRange(j, j, 'year'))
var max =filtered.reduce(ee.Reducer.max())
var time = ee.Date.fromYMD(j,i,1)
time = ee.Number.parse(time.format('YYYYMMDD'))
return max.set('system:time_start', time)
        .aaddBands(ee.Image(time).rename('time'))
});
return filtered_col;
 };

 var img_col = ee.ImageCollection(months.map(map_m).flatten());
 print(img_col,'image collection');

 // Removes null images
 var ndvi = img_col
 .map(function(image) {
  return image.set('count', image.bandNames().length())
 })
 .filter(ee.Filter.gt('count', 1))
 print('final images: ', ndvi)

 var title = {
  title: 'NDVI anomaly over time',
  hAxis: {title: 'Time', format: 'mm-yy'},
  vAxis: {title: 'NDVI'},
  };
print(ui.Chart.image.series(ndvi.select('NDVI_max'),
  roi,
  ee.Reducer.max(),
  30,
 'system:time_start')
   .setOptions(title));


texts: Making Fishnet that is specific size in square meters on ground?

texts: I have a raster image in GCS WGS 1984. For sampling I want to take a grid of 500 m^2 in my area of interest and pick the centroid of this square. How do I make sure that the grid created with Fishnet tool is 500 m^2 on the ground? I tried with the Cell width parameter but ended up with very large grids. I can't specify the number of rows and columns because the numbers can vary as I proceed with my project.


texts: Although I have seen many examples about doing this at ArcGIS Desktop, I am still trying to find a how to using QGIS. What I want to do is something like the image below.


texts:
Is there a way to fill while editing in ArcGIS?

You know similar to the popular paint tool, where an enclosed area is filled with a colour. I'd like to fill an enclosed area (enclosed by poly lines or polygons or both).

I'd like to be able to do this in an edit session, a batch process or tool would not provide the functionality that I am looking for.

Here is an image to better explain the situation:

I need to fill the area between the High Water Mark (the dark blue polyline) and the lake (the blue polygon). I've highlighted the area I need to fill with a red scratch. Now there is a lot of gaps and inaccuracies in the data, so I need to manually check each one as I do it. It takes forever to trace first the watermark and the lake. I don't see why there isn't an ArcGIS editing tool to fill in gaps. If Microsoft can create one anyone can.


texts:
I am looking for a DTM for the city of Agadir in Morocco.

Could anyone help me to find it?


texts:
I have a set of points with datatype as geometry(MultiPoint, 32643), I want each point to be converted in short line segments. The purpose of doing this is to use "split by lines" tool in QGIS which splits a line layer with another line layer.

The suggestion in one of the answer is as below, but I could not find the geometry by expression tool. Also, I tried to implement the below command in PostgreSQL but it didn't work.

Convert your points into short line segments with the geometry by expression tool.
Use an expression like this:

make_line($geometry, translate($geometry, 1, 1))


texts: Should I install all plugins to my new installation

texts: I am new to QGIS and I am preparing for first time use. Should I install all plugins that are not installed?


texts: addParameter to raster in Python QGIS plugin

texts:
I create a plugin for processing two rasters. I want to set two input layers and one output layer.

class CombineRaster:

"""QGIS Plugin Implementation."""

    Input = 'Input raster'
    Input2 = 'Input second raster'
    Output = 'Output'

def __init__(self, iface):
     self.addParameter(QgsProcessingParameterRasterLayer(
         self.Input,
         self.tr("Input raster"), None, False))
     self.addParameter(QgsProcessingParameterRasterLayer(
         self.Input2,
         self.tr("Input second raster"), None, False))
     self.addParameter(QgsProcessingParameterRasterDestination(
         self.Output,
         self.tr("Output"),
         None, False))

When loading the plugin, I get an AttributeError error: 'CombineRaster' (my plugin) object has no attribute 'addParameter'. What is wrong?


texts: How to include altitude using QGIS for field surveys capturing lat/long for utility locations?

texts:
I am using QGIS with a USB GPS to capture lat/long of utility locations.
I have created a .csv file with a few fields e.g. address, cross st, id number, etc. and made it "editable".
When I locate a utility item I utilize the the live GPS "add point" feature.
QGIS adds the Lat/long (X/Y) to the first two columns of the file and "prompts" for the information from my .csv file.
I can then save the file as a .csv and I have the attributes with lat/long.

Is there an "easy" way to have QGIS add/include the altitude "Z" along with the lat/long? Would this require a special script, etc.?


texts: Arcpy Script to find Field Values Begining with a Certain Character

texts:
I'm trying to find and update a list of all the field values for alternate pin numbers (syntax being: XX-XX-XXX-XXX).  I want to cycle through the pin numbers to find all of the pins beginning with 14 and 34 (14-XX-XXX-XXX, 34-XX-XXX-XXX).

I am needing to update those parcel pins to have a township number in front of the pins.  (EX: 01-14-XX-XXX-XXX and 01-34-XX-XXX-XXX)

This is what I have so far:

import arcpy
#Set workspace
arcpy.env.workspace = "C:\prj\IL\Bond\geodb\Bond_Master.mdb"

#Define input features and shapefile tables
ParcelPolygon = r"C:\prj\IL\Bond\geodb\Bond_Master.mdb\Cadastral\Parcel_Poly"

#Define field for cursor
altpin = "ALTPin"

#Create update cursor
with arcpy.da.SearchCursor(ParcelPolygon, (altpin)) as cursor:
for row in cursor:
    try:
        first_index_value = int(row[0].split("-")[0])

        if first_index_value == 14:
            print "{0} - This parcel is in Tamalco township".format(row[0])
        else:
            print "{0} - This parcel is NOT in Tamalco township.".format(row[0])

    except Exception as error:
        print "Value: {0} is invalid\nError: {1}".format(row[0], error)


texts: QGIS 3.4 distance matrix returns the wrong point

texts:
I'm using Distance Matrix to find the closest point between one set (4 points) and another set (4 points), detailed below. All results look okay except one, where the returned point is not the closest. Any ideas what is it that I'm doing wrong?

Set 1 (input)

Set 2 (target)

The output layer shows the following. However, TargetID for 1 should be 3, and not 2. By the way, I did another distance matrix layer where the "Use only the nearest (k) target point" was 0, giving an output of all distances between all points, and indeed the distance between 1 (bilbao) and 3 (oviedo) was smaller than the distance between 1 (bilbao) and 2 (alcala de henares), which suggests that qgis knows the correct distance between points but returned a wrong one when I'm looking for the nearest point only.

Distance Matrix layer


texts: Difference between GCS_European_1950 and GCS_International_1924

texts:
In ArcGIS, there are two datums they look like same except their names; GCS_European_1950 and GCS_International_1924. Where I started to work, I found out that they are using both datums. Since they are two different datums, ARCGIS asks for Geographic Transformation Parameters and there are no predefined any parameters. But I realized somebody had defined a custom transformation parameters(with Coordinate_Frame Method) in the past and its parameters were defined as 0. Is it OK to do transformation this way between these two datum from GIS point of view?

It seems like use of GCS_International_1924 is not common and I am quite confuse about the difference between those two geographic coordinates?

GCS_European_1950
WKID: 4230 Authority: EPSG

Angular Unit: Degree (0,0174532925199433)
Prime Meridian: Greenwich (0,0)
Datum: D_European_1950
  Spheroid: International_1924
    Semimajor Axis: 6378388,0
    Semiminor Axis: 6356911,946127947
    Inverse Flattening: 297,0

GCS_International_1924
WKID: 4022 Authority: EPSG

Angular Unit: Degree (0,0174532925199433)
Prime Meridian: Greenwich (0,0)
Datum: D_International_1924
  Spheroid: International_1924
    Semimajor Axis: 6378388,0
    Semiminor Axis: 6356911,946127947
    Inverse Flattening: 297,0


texts: Joining data without having a common ID or same spatial location

texts:
I have two polygon shapefiles, which don't have same spatial locations, neither have common attributes. Is it anyhow possible to join these files in ArcGIS 10.3?

I know the ArcGIS supports two kinds of joining,  with common attribute or with same spatial location. But in this case what I actually need is to generate common feature Id for the two files. Is it possible to generate?

Not in a conventional way, I have generated labels for the polygons then created annotations with them. They brought the desired results for many, but not for all the polygons. They got same new Object IDs but few sequences doesn't match. Can't find the reason.


texts: Exporting all domains to table in ArcGIS Pro using ArcPy?

texts:
I want to export all domains from gdb to table. I know that i can do run domain to table tool but it will only export values of that domain.
I also know that there is a batch option but even for that you need to know domain names.

I used to have a python script for that but that  was in ArcPy not ArcGIS pro python. I was searching through different blogs but havent found anything yet.

Has anyone tried it on ArcGIS Pro?


texts:
Not sure if this question belongs (Please tell me if it isn't)

I have an AutoCAD drawing that is in centimeters, and every time I insert a dimension, it appears in centimeters (which makes sense since the unit is set to centimeters in "Drawing Utilities") but is there a way to make the dimensions appear in meters while keeping the unit system in the drawing in centimeters?

*I Use AutoCAD 2014


texts: Creating a view in ArcCatalog - ORA-01722

texts:
I'm creating a view in ArcCatalog 10.6 on an Oracle DB. Never created a view before but with some Googling, I've managed the SQL below. The view works but only partially. In ArcCatalog, preview shows the first 2000 records. When I click "move to end of table", it gives ORA-01722 error. This happens because the first two characters are not always numeric, some are alpha. I thought using ELSE would calculate them to '9999' but I guess this isn't the case. I need some help where if any of the first two characters are alpha, calculate the value as is, no CONCAT.

SELECT GRANTOR, SKETCHNUM, LIBER, FOLIO, CASE
WHEN SUBSTR(CONTRACT_, 1, 2) BETWEEN 00 AND 20
THEN CONCAT('20', GISEDITOR.EASEMENTS_ROW.CONTRACT_)
WHEN SUBSTR(CONTRACT_, 1, 2) BETWEEN 21 AND 99
THEN CONCAT('19', GISEDITOR.EASEMENTS_ROW.CONTRACT_)
ELSE '9999'END AS ENGPLANDWG
FROM GISEDITOR.EASEMENTS_ROW


texts: Generating sequential numbering for duplicate values in a field sorted by ObjectID

texts:
I am trying to figure out a solution to an field calculation problem. What I'm trying to do is automatically (using either ArcPy or the Field Calculator) generate sequential numbers, beginning at 1, in a new blank integer field (called 'Point_ID') for every recurrence of a value in a second field (called 'Line_ID'). The sorting of the sequential values in the second field ('Line_ID') will be based on the order in a third field ('FID'). Can anyone help me do this? I am very green when it comes to advanced field calculations and ArcPy (and Python in general, for that matter). So your explicitness is much appreciated.

Another way to put it: there are duplicate values in the Line_ID field. I would like to create a new field that counts up for every duplicate occurrence in the Line_ID field, with the sort order based on the FID field. So, if there are nine values in the Line_ID field that have the value "A," the new Point_ID field will go from 1 to 9, with the order based on the sort on FID.


texts:
I am working on a series of footpath routes and have been asked if it possible to calculate the total ascent that a walker would have climbed if the travelled along each route?

I know that it is simply not a case of max-min because the footpath is not just one single climb but has many stages where the path ascends and descends.

Presently, I have the data as a series of polylines and an underlying dtm in ArcGIS but do not have a spatial analyst extension.

I have GRASS GIS and QGIS too so quite happy to consider external methods if this will help answer the question.


texts:
Where do I find WorldView satelite imagery that can I download that are hard to find around ?

There are few that I know of are:

EarthExplorer

GLCF

Glovis


texts:
I am working with the image below and I am trying to find how many pixels are in the grey corners of the image, and I want to assign them a probability of 5%. The bottom layer is my DEM and the top layer is a raster of distances. I know how many pixels I have on my raster layer, however I am not sure of how I can get the pixels outside of the circles and how to assign those pixels a probability. I would appreciate any help I can get!

Image Link


texts: Graphical modeler calculator missing in QGIS3

texts: There used to be a calculator in the modeler-only tools in QGIS graphical modeler. But I notice that it is missing in QGIS3. How do I do numerical calculations as part of the model? For example, take input in km and convert it to meters for certain algorithm.


texts:
I have a fresh install of QGIS 2.4. When trying to open a .shp file I get the error message below.

I've followed their advice and have done several re-installs and also cleared all qgis files from my system (Parallels on Mac OS X 10.6) before those re-installs. I've tried QGIS 2.2 which caused the same problem.

Funnily all worked for a day, the following day and since not anymore...

Has anyone come across this before?


texts: Buffer with dissolve - Geopandas - unary_union multipolygon

texts:
I am trying to buffer points which are quite close together, and then dissolve the overlapping buffer polygons into larger polygons. After doing the buffer is where I have the issue. I have tried using dissolve('geometry') but that gives this error.

KeyError: 'geometry'

I think the proper way to do this is using unary_union but my output is a multipolygon and I cant use it.

gdf2 = gpd.read_file(SELECT)
gdf2b = gdf2.buffer(gridDist)
gdfu = gdf2b.unary_union
gdfu.plot()

AttributeError: 'MultiPolygon' object has no attribute 'plot'

Or if I try to send the unary_union to a file

gpd.GeoSeries.to_file(gdf3u, BUFF1, driver="ESRI Shapefile")

AttributeError: 'MultiPolygon' object has no attribute 'index'

So my question is how can I either dissolve these polygons based on their self intersecting geometry and have the result be a geodataframe or a geoseries?
Or how can I transform the multipolygon I have from into a geodataframe?


texts:
I have a huge dataset (millions of points) and a irregular grid (quadtree) in shapefile format and I need to convert the points in countes over the grid provided. I managed to make it work, using shapely Polygon.contains and a naive double-loop that scales as O(n**2).

Obviously, the script is very slow. It takes already very long on thousands of points. Is there any other way to greatly reduce running time, without writing a lot more code?

What's the typical procedure to convert so many points to a grid system, and how is it called, so that I can search further information about it?


texts: Geotools: Unexpected result reprojecting bounding box to EPSG:3035

texts:
I have a bounding box given in EPSG:4326 and I want to reproject it into EPSG:3035. How should I be doing it?

Here’s what I tried:

try {
    CoordinateReferenceSystem sourceCRS = CRS.decode("EPSG:4326");
    CoordinateReferenceSystem targetCRS = CRS.decode("EPSG:3035");
    MathTransform trans4326to3035 = CRS.findMathTransform(sourceCRS, targetCRS);
    ReferencedEnvelope bounds4326 = new ReferencedEnvelope(5.0, 15.0, 47.0, 55.0, sourceCRS);
    Envelope bounds3035 = JTS.transform(bounds4326, trans4326to3035);
    System.out.println(bounds3035);
} catch (Exception e) {
    e.printStackTrace(System.err);
}

Result:

Env[-917569.8210525657 : 536626.7649395498, 8365819.661985753 : 9500037.05372034]

But these coordinate values are not what I expected. Cross-checking in QGIS, the values should be approximately 3900000 : 4700000, 2600000 : 3600000. This makes sense, because the center of projection of EPSG:3035 is 10°E, 52°N and the false easting/northing is 4321000/3210000, so the reprojected bounding box should be around that.

My Geotools version is 15.1.


texts: Server side clustering - cluster (feature) geometry

texts:
I am trying to implement server side clustering and show it using OpenLayers. I am using this good material as an example. As a clustering result, a get some list of  features, some of them are clusters. Like this (it is an output of previous example):

Array (
    [0] => Array (
        [0] => Array ( [id] => marker_3 [lat] => 59.431602 [lon] => 24.757563 )
        [1] => Array ( [id] => marker_4 [lat] => 59.437843 [lon] => 24.765759 )
        [2] => Array ( [id] => marker_6 [lat] => 59.434776 [lon] => 24.756681 )
    )
    [1] => Array ( [id] => marker_5 [lat] => 59.439644 [lon] => 24.779041 )
    [2] => Array ( [id] => marker_2 [lat] => 59.432365 [lon] => 24.742992 )
    [3] => Array ( [id] => marker_1 [lat] => 59.441193 [lon] => 24.729494 )
)

where 0 element is cluster containing 3 features.
But how can I compute this clusters geometry? Or how can I visualize it using OpenLayers?


texts: Does anyone success installed Mask R-CNN IN GRASS GIS?

texts:
As title, does anyone success installed Mask R-CNN in GRASS GIS?
Also, I found a script for updating GRASS GIS in order to use MASK R-CNN as below:
To use environment setting flags like –overwrite, the user has to update his GRASS GIS with the following patch:
lib/python/script/core.py:
    elif var.startswith (b'opt_'):
        options[var[4:]] = val
    elif var in [b'GRASS_OVERWRITE', b'GRASS_VERBOSE']:
-       os.environ[var] = val
+       os.environ[var.decode("utf-8")] = val.decode("utf-8")
    else:
        raise SyntaxError("invalid output from g.parser: %s" %line)

How can I update? I'm new to GRASS GIS.


texts: Convert multiple lines into a polygon in QGIS

texts:
I am working on a huge CAD-map that I am supposed to convert into a GIS map. Currently I am trying to convert the symbols used for survey points (see below) into polygons, of which I could calculate a central point as the real representation in GIS.

I got rid of all the unnecessary lines, and wanted to convert the remaining line octagon into a polygon by use of the "Line to Polygon tool", but that does not work, apparently, because every edge is a single line (see below). However the lines appear to be correctly snapped to each others end points.

I tried several workarounds, but nothing really worked. Any ideas?


texts:
I have done a lot of research and cannot seem to get the following expression right. I have three (3) attribute columns - ID, Diameter, and Material, I would like to label sewer lines with the unique id, diameter, and material. E.g. GS1945 12" VCP. Ultimately the script would leave out any information the is NULL.

def FindLabel ([ID],[Diameter],[Material])
     if IsNull([Diameter]) and IsNull([Material]) then
          return [ID]
     elif IsNull([Diameter]) then
          return ([ID] + " " + [Material])
     elif IsNull([Material]) then
          return ([ID] + " " + [Diameter] + "''" )
     else:
          return ([ID] + " " + [Diameter] + "''" + [Material])


texts: QGIS Custom Python Function to Return Current Extent QGIS 2.18LTR

texts:
Related Question Charts and Graphs that only display features in current extent
Trying to write a custom function to return the current view extent. The below code is what I wrote based on the QGIS Cookbook https://docs.qgis.org/testing/pdf/en/QGIS-testing-PyQGISDeveloperCookbook-en.pdf. The return from the Output Preview is
<built-in method xMinimum of QgsRectangle object at 0x000000...
Any suggestions on where I am going wrong?

Function Code:

import math
from qgis.core import *
from qgis.gui import *
from PyQt4.QtCore import QVariant
from PyQt4.QtGui import QAction, QMainWindow
canvas=QgsMapCanvas()
e=canvas.extent()
@qgsfunction(args=0, group='Custom', usesgeometry=True)
def CurrentExtent(value1, feature, parent):
    return str(e.xMinimum)+" "+str(e.yMaximum)+", "+str(e.xMaximum)+" "+str(e.yMaximum)+", "+str(e.xMaximum)+" "+str(e.yMinimum)+","+str(e.xMinimum)+" "+str(e.yMinimum)+"))"


texts: Some labels not displaying in ArcGIS Desktop?

texts:
I am using ArcGIS Desktop 10 and trying to label line features. It labels 75% of the features normally, while others do not display any label.

There is no obvious geographic or label correlation between the malfunctioning features and they are all part of the same shapefile. If I alter the label slightly (for example: changing 12" to 12'') I can often get the label to display, but this isn't a practical solution.


texts: Getting precipitation data over each point in a Feature Collection in Google Earth Engine?

texts:
I want to generate a time series chart of precipitation data using CHIRPS in Google Earth Engine. I have some twenty points over which I want the data individually, i.e, I want twenty charts.

So far, I have managed to get a single chart over all the points altogether. How can I modify the following code as per my requirements?

Link to the code: https://code.earthengine.google.com/f38eca09db9e303f69d5a69e9d863a0a

//Filter image collection
var precip = chirps.filterDate('2000-01-01','2008-12-31');

//Create and print rainfall chart
print(ui.Chart.image.series(precip.limit(5000),tunga,ee.Reducer.mean(),1000).setOptions({

title: 'gp_1: PPT time series',
hAxis: {title: 'Day'},
vAxis:{title:'rainfall (mm/day)'}
}));


texts: Area measurement validation: ellipsoid vs geoid

texts:
I try to validate area measurements on orthomosaics using objects of known size on the ground. I'm working on centimeter level measurements.
As far as i understand an orthomosaic in WGS 84 is projected to the reference ellipsoid surface.
I set up a square of 1 m² on the ground which is level. Because it is level it means that it is parallel to the geoid (or more precicely to the local surface normal vector of the geoid).
By defintion ellipsoids are only an approximation of the geoid, which means the ellipsoid plane is not necessarily parallel to the water-level. This in turn would result in a error in area measurement.
Because RTK GPS systems usually use a bubble level for setup i would assume that this error is small. But how small is it? Are there locations on earth were this is a problem?


texts:
With the following script I was able to determine the maximum values of a column [POINT_Z] for certain groups from [ID_1].

    arcpy.MakeFeatureLayer_management(inLAS[:-4] + '_POLE_POINTS.SHP', "POLES")
    reset_layers("POLES")

    with arcpy.da.SearchCursor(inLAS[:-4] + '_POLE_POINTS.SHP','ID_1') as max_cursor:
        for row in max_cursor:
            mast_value = sorted({row[0] for row in max_cursor})

    #iterate through the list to find highest value in field [pointz], but only in features which have the value of the element of the list "mast_value"
    field = "ID_1"
    for mast in mast_value:
        expression = field + " = " + str(mast)
        with arcpy.da.SearchCursor(inLAS[:-4] + '_POLE_POINTS.SHP',('POINT_Z', 'FID'), where_clause=expression) as mast_cursor:
            for row1 in mast_cursor:
                print('ID ' + str(mast) + u': max altitude ' + str(max(mast_cursor)))
                arcpy.SelectLayerByAttribute_management("POLES", "NEW_SELECTION", "\"FID\" = " + str(max(mast_cursor)[1]))

With the output I can display the results. A manual check showed that the result is correct.

ID 1: max altitude (477.56, 5997) ID 7: max altitude (523.27, 3401) ID 16: max altitude (495.69, 3203) ID 38: max altitude (495.69, 5574) ID 39: max altitude (470.16, 5986) ID 44: max altitude (503.93, 4713) ID 63: max altitude (503.97, 4185) ID 73: max altitude (463.67, 4205) ID 75: max altitude (518.59, 2880) ID 77: max altitude (463.63, 5877) ID 79: max altitude (485.91, 5594) ID 87: max altitude (482.06, 853) ID 88: max altitude (475.2, 1836) ID 91: max altitude (492.0, 1074) ID 97: max altitude (468.77, 1323) ID 98: max altitude (462.52, 1546) ID 99: max altitude (463.81, 1294)

To continue working (selection), I need the corrosponding value [FID] (second value in the brackets) of the maximum value of [POINT_Z] as a variable. My attempt to write the result into a variable

myMastList = str(max(mast_cursor))

fails with an error message.

File "E:\Documents\LAS\Las_Wire.py", line 264, in masten
    myMastList = str(max(mast_cursor))
ValueError: max() arg is an empty sequence

How can I write the value of [FID] from the highest value of [POINT_Z] (or both [FID] and [POINT_Z]) into a variable?


texts:
I have a raster file that is projected not detailed enough. Every pixel in the raster should be replaced, just a bit.
With ArcGIS I calculated the offset per pixel. I couldn't find a way to do this offsetting in ArcGIS, so I switched to python. I found out that I can get the current envelope with ogr with :

env = geom.GetEnvelope()
print "minX: %d, minY: %d, maxX: %d, maxY: %d" %(env[0],env[2],env[1],env[3])

And now I have to change these to the (already calculated) new origin.

Here is some sample code:

    ## From gdal info page
wkt = "LINESTRING (1181866.263593049 615654.4222507705, 1205917.1207499576 623979.7189589312, 1227192.8790041457 643405.4112779726, 1224880.2965852122 665143.6860159477)"
geom = ogr.CreateGeometryFromWkt(wkt)
env = geom.GetEnvelope()

I have to say that I don't entirely understand the envelope function. Is this a function that returns the minX etc from the projection? Or does it just calculates this all from the geometry?

Does anybody now how to change this or knows another workaround?


texts: pyProj Coordinate Transformation Incorrect

texts:
I'm using pyProj to transform extents between two different geographical systems. In a case when both my systems are the same I expect to get back the same extents as inputted but this isn't happening for the below mentioned case:

import pyproj
inProj = pyproj.Proj("+init=EPSG:{0}".format(3857)) # Mercator
outProj = pyproj.Proj("+init=EPSG:{0}".format(3857))
x =  (-20037507.0672, 20037509.6184)
y =  (-1467048.29156, 8625918.8737)
print pyproj.transform(inProj,outProj,x,y)

The output comes out to be : ((-20037507.0672, -20037507.067178484), (-1467048.2915600014, 8625918.873699998))

Which is incorrect as my X-Coordinate extents are same for xMax and xMin.

If I change my co-ordinate systems

import pyproj
inProj = pyproj.Proj("+init=EPSG:{0}".format(3857)) # Mercator
outProj = pyproj.Proj("+init=EPSG:{0}".format(4326)) # WGS84
x =  (-20037507.0672, 20037509.6184)
y =  (-1467048.29156, 8625918.8737)
print pyproj.transform(inProj,outProj,x,y)

My X-coordinate output comes out to be : (-179.99998854118687, -179.99998854099357)

Which is again incorrect. Is this due to a domain issue in pyProj or is there any other reason ?


texts: Label QGIS diagram sector with percentage etc

texts:
I'm working with QGIS 2.14.3, using the diagrams function for vector layers. Is it possible to label each single diagram sector with an attribute / calculated percentage? In the layer properties tab / labels I only see an option to label the whole feature. In the diagrams menu, I can't find any labelling option.

It would be very nice to label for example every single pillar in this diagram with an attribute / calculated percentage. The displayed label is the label of the whole feature.


texts:
Using ArcMap 10.3, I have a shapefile with dissolved roads and I have 10 records with the Road ID field (they are all unique). I want to join an Excel file with Road-ID with 15 records and I want to append Geometry. for example:

So I want to have Road ID 10  with values 5 and 10 after join in the dissolved layer.
I guess it is  Many to one join. So Actually I want to copy the record, 2 times in case Road Id 10 is two times in the Excel file and so on. Right join is not working, here.
Should I use iterators?


texts:
According to Using QGIS3 Processing algorithms from standalone PyQGIS scripts (outside of GUI) I try to implement the existing "Split Lines with Lines" algorithm into a QGIS Plugin.

To do this I need the path, where the algorithm file is stored, but I do not manage to find it. Can you tell me where it is stored and if there is an easy way to find out, where algorithms are stored?


texts: R-Script is not working from QGIS processing toolbox (error massage)

texts:
I would like to use some r-scripts from processing toolbox but I get this error massage:

The downloaded source packages are in
â/tmp/RtmpdP9CQf/downloaded_packagesâ
Warning message:
In install.packages("rgdal", dependencies = TRUE) :
installation of package ârgdalâ had non-zero exit status
tryCatch(find.package("raster"), error=function(e) install.packages("raster", dependencies=TRUE))
[1] "/usr/local/lib/R/site-library/raster"
library("raster")
Loading required package: sp
library("rgdal")
Error in library("rgdal") : there is no package called ârgdalâ
Execution halted

I installed r-base and r-base-dev 3.4 from CRAN-repository.
R is activated as provider in my QGIS 2.18.9 processing toolbox.

EDIT: I also tryed to install RGDAL via R and got a similar error:

downloaded 1.6 MB
* installing *source* package ‘rgdal’ ...
** package ‘rgdal’ successfully unpacked and MD5 sums checked
configure: CC: gcc -std=gnu99
configure: CXX: g++
configure: rgdal: 1.2-7
checking for /usr/bin/svnversion... no
configure: svn revision: 660
checking for gdal-config... no
no
configure: error: gdal-config not found or not executable.
ERROR: configuration failed for package ‘rgdal’
* removing ‘/usr/local/lib/R/site-library/rgdal’

The downloaded source packages are in
    ‘/tmp/Rtmp3y3vA3/downloaded_packages’
Warning message:
In install.packages("rgdal") :
  installation of package ‘rgdal’ had non-zero exit status
Error in library(rgdal) : there is no package called ‘rgdal’


texts: Displaying point data that is close together using ArcGIS Pro

texts:
I am using ArcGIS Pro to display point data on a map. Each of these points need to be displayed with clearly visible symbols. However, some of them sit very close together, and when displaying on a printed map, they appear as a collective black blob, rather than individual symbols. When I make them small enough that they are separate points, they appear only as tiny dots, not as unique symbols.

How can I best display these points?

Is annotating the symbol itself an option? (as an aside, I do have two inset maps to display the busiest areas, but this is still a problem on the rest of the map).


texts: Working with DEMs: which should I do first calculate slope or reproject?

texts:
I have some 10 m digital elevation model rasters in UTM, and I need to generate 30 m resolution slope rasters in Albers.

Which should I do first: reproject the rasters, or calculate slope?

In case it matters, I am planning to use gdalwarp and gdaldem to do this.


texts: What should QGIS course for advanced ArcGIS Desktop users cover?

texts:
I will build a 16 hour (2 day) course on QGIS for ArcGIS users. The target users that breathe and sleep ArcGIS. As we're trying to reduce our costs, we're thinking of switching some people to QGIS and having a floater ArcGIS license rather than have everyone with it. I've used QGIS extensively and very familiar with it (versions 0.7 to 1.8)

Are there any course materials that are specific to "QGIS to ArcGIS" or that target advanced operations in QGIS?


texts: How can I extract the ground from a DTM?

texts: I need to create a shapefile that only includes the vegetation in the scene. I was thinking I could take the DTM and the DSM and subtract the DSM from the DTM. Is there a way to do this in ArcMAP?


texts: Rename bands in a multiband raster in ArcGIS 10.x

texts:
I have imported a multiband raster in .tif format into a geodatabase using ArcGIS 10.7 and want to rename the bands to something more descriptive. However, I can't find out how to do this easily.
Does anyone know how to rename bands in ArcGIS 10.7?


texts: Checking Shapefile type (Point, Line, Polygon) in ModelBuilder using python script?

texts:
I am creating a model and would like the input to branch off to different tools within the model based on the shapefile type (point, line, polygon). I haven't found a tool in ModelBuilder that will do this and am not very skilled in python.

Is there a tool I'm not seeing or a script that exists for this?


texts: Storing entire GeoJSON document in PostgreSQL/PostGIS?

texts:
I am trying to insert GeoJSON document in my PostgreSQL, I cannot just simply use the data type Geometry and use the expression ST_GEOMFROMGEOJSON because my geoJSON looks like this, and I need to store it all, what do you suggest?

GeoJSON Document

{
"type": "FeatureCollection",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::26916" } },

"features": [
{ "type": "Feature", "properties": { "h": -0.400000, "mean": null }, "geometry": { "type": "MultiPolygon", "coordinates": [ [ [ [ 337120.0, 4659412.460953020490706 ], [ 337121.0, 4659412.271286088041961 ], [ 337121.32677364279516, 4659412.0 ], [ 337121.581764082075097, 4659411.0 ], [ 337121.0, 4659410.470878113992512 ], [ 337120.0, 4659410.360407529398799 ], [ 337119.0, 4659410.916399526409805 ], [ 337118.937826662964653, 4659411.0 ], [ 337119.0, 4659411.488602316007018 ], [ 337119.127736570662819, 4659412.0 ], [ 337120.0, 4659412.460953020490706 ] ] ] ] } }
]}

There are more lines but I've cut it down. I have tried inserting it as text only, but then when I wanted to export the result it offered me CSV, XHTML,tabbed or XML and they all start putting  random " marks. so I am guessing Text data type is not the way. I am sure someone out there tried inserting complete GeoJSON not only the geometry part.

bytea the .geojson file?


texts: Feature Class to Feature Class renames field, but not alias

texts:
I have the following model:

In my two FC to FC tools shown above, I rename many fields. For example, table 1 has field "MEAN" and I rename it "approval_rate." Table two has "approval_rate" and I rename it "approval_rate_old." I join those two tables with renamed fields, and then run calculations on the fields by their new names (starting with Calculate Field (2)).

However, Calculate Field (2) fails, even though I have the syntax and field names correct in my simple VB expression. I think it is failing because FC to FC (2) is not changing the alias, just the field name. So when I open the resulting attributes table, it still shows all the previous names. In the example above, this leads to two fields with the alias "approval_rate" in the same table. So when Calculate Field (2) goes to run an expression using "approval_rate" it does't know which field to use and thus fails with "General Exception Error."

Does my hypothesis seem correct? Why is one FC to FC tool changing the name and alias, and one is only changing the name, not the alias?

Here is the screenshot of my Calculate Field (2) tool.

And here are screenshots of the relevant columns in the input table.


texts:
I have a table map_building_data and I create view all_incidents_unions

CREATE OR REPLACE FORCE VIEW ALL_INCIDENTS_UNIONS
(
   ORA_GEOMETRY
)
AS
SELECT sdo_aggr_union (MDSYS.sdoaggrtype (ora_geometry, 0.005))
             ora_geometry
       FROM (  SELECT sdo_aggr_union (MDSYS.sdoaggrtype (ora_geometry, 0.005))
                         ora_geometry
                 FROM (  SELECT sdo_aggr_union (
                                   MDSYS.sdoaggrtype (ora_geometry, 0.005))
                                   ora_geometry
                           FROM (  SELECT sdo_aggr_union (
                                             MDSYS.sdoaggrtype (ora_geometry, 0.005))
                                             ora_geometry
                                     FROM (  SELECT sdo_aggr_union (
                                                       MDSYS.sdoaggrtype (ora_geometry,
                                                                          0.005))
                                                       ora_geometry
                                               FROM (  SELECT sdo_aggr_union (
                                                                 MDSYS.sdoaggrtype (
                                                                    ora_geometry,
                                                                    0.005))
                                                                 ora_geometry
                                                         FROM (  SELECT SDO_AGGR_UNION (
                                                                           SDOAGGRTYPE (
                                                                              COM.GEOM,
                                                                              0.005))
                                                                           ora_geometry
                                                                   FROM MAP_BUILDING_DATA COM
                                                                  WHERE ID =9307433
                                                               GROUP BY MOD (ROWNUM, 128))
                                                     GROUP BY MOD (ROWNUM, 64))
                                           GROUP BY MOD (ROWNUM, 32))
                                 GROUP BY MOD (ROWNUM, 16))
                       GROUP BY MOD (ROWNUM, 8))
             GROUP BY MOD (ROWNUM, 4))
   GROUP BY MOD (ROWNUM, 2));

After create view I am inserting data on USER_SDO_GEOM_METADATA

INSERT INTO USER_SDO_GEOM_METADATA (TABLE_NAME,COLUMN_NAME,DIMINFO,SRID)VALUES ('ALL_INCIDENTS_UNIONS','ORA_GEOMETRY',MDSYS.SDO_DIM_ARRAY(MDSYS.SDO_DIM_ELEMENT('X',8255000,8754700,0.005),MDSYS.SDO_DIM_ELEMENT('Y',4261000,4614700,0.005)),'2400000');

After this when I create layer on geoserver and want see result, this return error

Rendering process failed java.io.IOExceptionORA-13226: interface not supported without a spatial index ORA-06512: at "MDSYS.MD", line 1723 ORA-06512: at "MDSYS.MDERR", line 8 ORA-06512: at "MDSYS.SDO_3GL", line 1224

please help ..........


texts: Geoserver/geotools DB2 extension support

texts:
I'm evaluating the migration of a postgis/postgres database to a DB2 spatial db.
A Geoserver instance is connected to the postgis datasources, so I know there is an extension available for making the same connections from Geoserver to the DB2 datasources.

I've been looking for an up-to-date page or documentation where I can find the supported versions of DB2 or the DB2 drivers for each gt-jdbc-db2 release (the library behind the extension AFAIK), but I haven't been able to find any info about this.

In the end, what I need to know is the versions of DB2 supported by the extension.


texts: Displaying rectified raster (topos) in Google Earth

texts:
Looking for tool or an easy way (limited programming skills here) to display geotiffs in Google Earth. Specifically Canadian topos like this

This Geogratis site has the same topos in regular tif and pdf formats but can't seem to locate a way to display them in Google Earth. I found this script but don't know how to implement it.
Any help would be appreciated!


texts: Using a spatial query to filter the atlas pages in QGIS

texts:
I would need to create an atlas of all the map sheets (red rectangles) that refer to a specific municipal territory (for example the selected yellow polygon) as in the image. The atlas coverage layer is the map sheet layer.
I try to understand if there is a possibility to use a simpler and more automatic way than using an expression like

"Map_Sheet_number" IN ('459090', '459100', 459130, '459140')

as in the example, to filter the atlas coverage layer.

Is it possible (and how) to use a spatial query to filter the atlas coverage layer (map sheets) using another layer (municipal territory)?


texts: How to find Nearest point on line from a point cordinate in SQL server

texts:
I want to find out the nearest point on my road network which is a multiline string from a point with lat lon value.

The below code I have tried for finding buffer distance line from a point.

DECLARE @a geometry;
SET @a = geometry::Point(73.0127941495,19.26388944035,4326);
    SELECT top 1 * FROM Routes
    WHERE geom.STIntersects(@a.STBuffer(0.002)) = 1 order by @a.STDistance(geom);

I want a solution for finding Nearest point on line from a point cordinate


texts: Overlay .mxd file to ArcGIS in Android Studio

texts:
I'm trying to overlay a floor plan of a building to the app I've built using the site:

https://developers.arcgis.com/android/10-2/guide/add-geocoding-to-your-app.htm

I have path segments and nodes drawn inside the image that I would like to be able to use my phone to navigate from my GPS location to whatever room number I specify (which has a node). I can't publish the file to any ArcGIS server because I have signed an NDA concerning the CAD files for the floor plans.

In short, how can I overlay the .mxd file onto the ArcGIS map in my Android application?


texts:
I would like to use log4net as the logging framework for my addin. However, now I'm having trouble in enabling logging in debug mode and in the released esriAddin file. It works in either one, not both at the same time.

I have a dedicated configuration file for log4net (called log4net.xml). When the Build Action of the file is set to Content, the logging works in the debug mode but not in the released addin. When Build Action is AddinContent, it's the opposite.

The reason for this is that when in debug mode, the solution uses the bin folder as the base directory, but when the addin is compiled the folder is the "esriaddin" file.

So when in debug mode and Build Action=Content" the log4net.xml is at C:\Users\%username%\AppData\Local\ESRI\Desktop10.0\AssemblyCache\{9FFCCA4E-6DD6-81B5-17C3-7683C725F70D}. When I change Build Action to "Addin Content" (while still in debug mode), the log4net.xml is not copied to the directory!

My problem is how can I configure my solution so that the log4net config file would be in a static location regardless of the Build Action? I need to do this because I want to automate the building process and changing the Build Action manually is not an option.

I've checked this thread but its solution is to edit the ArcMap.exe.config file which is not an option for me (many users for the addin, the file may not be empty).

Any suggestions?


texts: How to find the coordinates (longitude, latitude) of a point in an aerial image?

texts:
I will be receiving drone images from an associate. It'll have longitudes/latitudes values from the corners of the image.  I don't have access to his image processing S/W.  It'll be in a relatively small area 1,000,000m².  Not sure if I would need the accuracy of something like the Haversine formula.

I need to find the coordinates (longitude, latitude) of a given point (trucks for example) in the image.  This point is ever changing image to image.
I'm seeking an equation/Java code/Python code etc.


texts: Can QGIS Time Manager be used with millions of data points?

texts:
I have several months of tracking data for researchers driving across north america and i would like to animate their progress.  The Time Manager seems to work fine for maybe one weeks worth of one researchers data but if I use a much larger data set it freezes Qgis.  Is there a way to run the plugin through something other than QGis desktop like a command prompt so I can monitor progress and spot any data errors?

Maybe an alternative that can handle several million data points?


texts: Viewing pseudo node and dangle nodes using ArcPy?

texts: We know that we require production mapping extension and ARCINFO license to view the pseude node and dangle nodes.We have to see these nodes without these functionality through python scripting,but I do not have to use any in-built tool ,I have to make my own tool through python which can helps me viewing the nodes


texts: Is it possible to add values+symbols to an imported symbology in ArcMap?

texts:
I am trying to use an imported symbology from a layer file in ArcMap 10.2. However, the shapefile that I want to make a symbology for has two or three more values than are present in the imported legend. Whenever I want to add the new values and manually edit their symbols, the entire imported legend disappears and I end up with all new symbols.

Is it possible to add values and symbols to an imported legend?

I am using layer properties->categories->unique values in the layer properties window, same as in the layer file.


texts: Finding average elevations from a Raster DEM, applying those to Vector Grid created with GME

texts:
I'm using ArcMap 10.3, and here's what I have:
I have an angled vector grid created through GME's GenVecGrid tool which has been clipped to an irregular size and angled relative to the xy axes to represent a study area. When I say it's irregular, I mean to say while each cell is full (no partial cells), portions of the rectangular grid have been removed to reflect the study area. Here's an image to show what I mean:

I also have a raster DEM of the area in question (also pictured in the above image link).

Here's what I'm hoping to do:
For each cell in the vector grid, I want to find an average elevation from the raster DEM and incorporate this into the vector grid itself (in the attribute table). Ideally, I'd like to use an interpolation method like kriging to find an intelligent average, so to speak.

In searching forums, I have seen recommendations for using the Add Surface Information tool (ArcToolbox > 3D Analyst > Add Surface Information), as well as Interpolate Shape (ArcToolbox > 3D Analyst > Interpolate Shape), but these were unsuccessful.

Can anyone guide me in the right direction?


texts: How to generate grid shapefile with grid squares that align with landmasses in R

texts:
I have a hard time trying to describe the shapefile exactly, but I am trying to generate a grid over the Gulf of St. Lawrence, such that only the gridsquares which lie over water are included. I would like this to extend to the narrow southern strait at (-60, 47)

I'm really new to GIS in R, and below is my attempt at brute forcing it by generating rows of gridsquares and specifying the bounding box for each line individually. I used st_make_grid from the sf library to make the 5 rows.
Each gridsquare is 0.25x0.25 degrees.

I feel like there must be a more elegant solution to this... if anyone knows one.


texts: QGIS 3.0.3 topology checker to find gaps between polygons

texts:
When using QGIS 2.x I have used Topology checker to find gaps between polygons, crossing lines and duplicated points/vertex. I also liked that it shoved a list of errors where I could zoom into each error to fix it.

In version 3.0.3 I see that topology checker is running automatically and tells me at once when I put two point on top of each other.

How can I make it check for gaps? Is there somewhere to set the topology rules?


texts: I would like to set up an OSM routing server on my LAMP stack. I am interested in doing a lot of of routing queries.  I have a search function that is closely related to getting the driving distance between lat/long locations and will thus be handling a lot of routing queries, even in a development environment.  I know how to use the Mapquest API but would like to set up an OSM routing server before I sign up for a business licence with a proprietary mapping service.  Where would I start this project? I have looked at OSRM and osm4routing but am still a bit bewildered and can't find anything very well documented out there.


texts: Using ogr2ogr w/ -s_srs to combine inconsistent EPSG files into PostGIS?

texts:
I'm having trouble combining shapefiles from TIGER and different local governments into a single postgis db to query legislative locations.  I've imported all data to multipolygon geom, but found some datasets' polygon shapes returning really large floats values lat/lng.  Looking more closely, I noticed inconsistent EPSG coordinate systems, based on searching prj2epsg.org with contents of my .prj files:

#Census Data (EPSG:4269 - GCS_North_American_1983)
GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]

#Texas Data (EPSG:3081 - NAD_1983_Lambert_Conformal_Conic)
PROJCS["NAD_1983_Lambert_Conformal_Conic",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",1000000.0],PARAMETER["False_Northing",1000000.0],PARAMETER["Central_Meridian",-100.0],PARAMETER["Standard_Parallel_1",34.91666666666666],PARAMETER["Standard_Parallel_2",27.41666666666667],PARAMETER["Latitude_Of_Origin",31.16666666666667],UNIT["Meter",1.0]]

I can't seem to get ogr2ogr to explicitly transform the target coordinate systems with -a_srs EPSG:3081 -s_srs EPSG:3081 -t_srs EPSG:4269:

ogr2ogr -f PostgreSQL PG:"host=..." file.shp -nlt MULTIPOLYGON -s_srs EPSG:3081 -a_srs EPSG:3081 -s_srs EPSG:3081 -t_srs EPSG:4269 -lco GEOMETRY_NAME=the_geom -nln table_name.

I've tried a number of different variations on the above, but at best I get ERROR 1: AddGeometryColumn failed for layer table_name, layer creation has failed.

What am I missing here?

I'm still a bit of a novice on work with shapefiles.

Example source files:

From Texas: ftp://ftpgis1.tlc.state.tx.us/DistrictViewer/House/PLANH309.zip
From Tiger: ftp://ftp2.census.gov/geo/tiger/TIGER2013/CD/tl_2013_us_cd113.zip


texts: Access SpatiaLite connections with PyQGIS QgsProviderMetadata (QGIS 3.10.7 LTR)

texts:
I would like to list spatialite connections and access a specific connectionwith PyQGIS. For postgres provider it's pretty straightforward:
>>> QgsProviderRegistry.instance().providerMetadata('postgres').connections()
{'alkis (authcfg)': <qgis._core.QgsAbstractDatabaseProviderConnection object at 0x000001D663232CA8>, 'baumradar (authcfg)': <qgis._core.QgsAbstractDatabaseProviderConnection object at 0x000001D663232948>, 'eb77 (authcfg)': <qgis._core.QgsAbstractDatabaseProviderConnection object at 0x000001D663232E58>,...
>>> QgsProviderRegistry.instance().providerMetadata('postgres').connections()['alkis (authcfg)']
<qgis._core.QgsAbstractDatabaseProviderConnection object at 0x000001D663232CA8>

But the same thing for spatialite provider returns a traceback:
>>> QgsProviderRegistry.instance().providerMetadata('spatialite').connections()
Traceback (most recent call last):
  File "C:\OSGEO4~1\apps\Python37\lib\code.py", line 90, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
_core.QgsProviderConnectionException: Provider spatialite has no connections method

Same problem with e.g. wms provider, both are available:
>>> QgsProviderRegistry.instance().providerList()
['DB2', 'WFS', 'arcgisfeatureserver', 'arcgismapserver', 'delimitedtext', 'gdal', 'geonode', 'gpx', 'mdal', 'memory', 'mesh_memory', 'mssql', 'ogr', 'ows', 'postgres', 'spatialite', 'virtual', 'wcs', 'wms']

Shouldn't all providers accessed through QgsProviderMetadata()? How to list and access spatialite (and other) connections using this class?


texts: ArcGIS Geometric Network: match downstream segments with their ultimate upstream source

texts:
I have a geometric network that represents a combined sewer system. The network functions as expected - I can set flags and trace upstream and downstream.

My goal is to identify all of the ultimate sources (line-segments) that flow to each of several downstream weirs (also line segments). Because it is a complex branching network, there can be multiple source segments per weir and, theoretically, a source segment could terminate at multiple weirs. A final deliverable would be a table with two columns: Weirs and sources:

| Weir ID  | Source ID |
------------------------
|    1     |     A     |
|    1     |     B     |
|    2     |     C     |
|    3     |     D     |
|    4     |     A     |
|    1     |     E     |
------------------------

What I've figured out so far:

I can create a layer of points to represent the 90 or so weirs in the network.
I can use the point layer as flags to execute the "Trace Geometric Network" tool from ArcToolbox.
I can save out the traced flowpaths from source to weir as a new layer.

If I only had to do this for a single weir, it would be easy. I'd look at the traced network and write down the IDs of any source line segments. Or, if it's a particularly complex trace, I'd run a query on the saved-out flowpath and select any segments that don't have a segment upstream. Something like:

SELECT * FROM saved_flowpath sf1 WHERE NOT EXISTS
    (SELECT 1 FROM saved_flowpath sf2
     WHERE sf1.upstream_node_id = sf2.downstream_node_id)

However, with ~90 weirs in this particular network, and the likelihood of having to do this again for other networks, I want to find a way to automate it. But without being able to write recursive queries in the ArcGIS query builder, I'm at a loss for how to proceed. Anyone have any thoughts?


texts: Couldn't load PyQGIS...on Windows custom build

texts:
I have built QGIS 3.2.0 with the release code downloaded from GitHub, using CMake and Visual Studio 2017. Basically I followed the instructions in INSTALL file. When I run QGIS.exe, it shows:

Couldn't load PyQGIS.
Python support will be disabled.

Traceback (most recent call last):
  File "", line 1, in
ModuleNotFoundError: No module named 'qgis'

Python version:
3.6.0 (v3.6.0:41df79263a11, Dec 23 2016, 08:06:12) [MSC v.1900 64 bit (AMD64)]

QGIS version:
3.2.0-Bonn 'Bonn', exported

Python path:
['I:/OSGeo4W64/apps/qgis/./python', 'H:/Users/Aurus/AppData/Roaming/QGIS/QGIS3\\profiles\\default/python', 'H:/Users/Aurus/AppData/Roaming/QGIS/QGIS3\\profiles\\default/python/plugins', 'I:/OSGeo4W64/apps/qgis/./python/plugins', 'I:\\OSGeo4W64\\bin\\python36.zip', 'I:\\OSGeo4W64\\apps\\Python36\\DLLs', 'I:\\OSGeo4W64\\apps\\Python36\\lib', 'I:\\qgis\\bin', 'I:\\OSGeo4W64\\apps\\Python36', 'I:\\OSGeo4W64\\apps\\Python36\\lib\\site-packages']

And:

An error occurred during execution of following code:
qgis.utils.uninstallErrorHook()

Traceback (most recent call last):
  File "", line 1, in
NameError: name 'qgis' is not defined

Python version:
3.6.0 (v3.6.0:41df79263a11, Dec 23 2016, 08:06:12) [MSC v.1900 64 bit (AMD64)]

QGIS version:
3.2.0-Bonn 'Bonn', exported

Python path:
['I:/OSGeo4W64/apps/qgis/./python', 'H:/Users/Aurus/AppData/Roaming/QGIS/QGIS3\\profiles\\default/python', 'H:/Users/Aurus/AppData/Roaming/QGIS/QGIS3\\profiles\\default/python/plugins', 'I:/OSGeo4W64/apps/qgis/./python/plugins', 'I:\\OSGeo4W64\\bin\\python36.zip', 'I:\\OSGeo4W64\\apps\\Python36\\DLLs', 'I:\\OSGeo4W64\\apps\\Python36\\lib', 'I:\\qgis\\bin', 'I:\\OSGeo4W64\\apps\\Python36', 'I:\\OSGeo4W64\\apps\\Python36\\lib\\site-packages']

I place a qgis.env on the binary directory with following contents, instead of modifying system environment variables:

PATH=I:\OSGeo4W64\apps\qgis\bin;I:\OSGeo4W64\apps\Python36;I:\OSGeo4W64\apps\Python36\Scripts;I:\OSGeo4W64\apps\qt5\bin;I:\OSGeo4W64\apps\Python27\Scripts;I:\OSGeo4W64\bin;H:\WINDOWS\system32;H:\WINDOWS;H:\WINDOWS\system32\WBem
GDAL_DATA=I:\OSGeo4W64\share\gdal
GDAL_DRIVER_PATH=I:\OSGeo4W64\bin\gdalplugins
GDAL_FILENAME_IS_UTF8=YES
GEOTIFF_CSV=I:\OSGeo4W64\share\epsg_csv
JPEGMEM=1000000
OSGEO4W_ROOT=I:\OSGeo4W64
PROJ_LIB=I:\OSGeo4W64\share\proj
PYTHONHOME=I:\OSGeo4W64\apps\Python36
QGIS_PREFIX_PATH=I:/OSGeo4W64/apps/qgis
QT_PLUGIN_PATH=I:\OSGeo4W64\apps\qgis\qtplugins;I:\OSGeo4W64\apps\qt5\plugins
VSI_CACHE=TRUE
VSI_CACHE_SIZE=1000000
O4W_QT_PREFIX=I:/OSGeo4W64/apps/Qt5
O4W_QT_BINARIES=I:/OSGeo4W64/apps/Qt5/bin
O4W_QT_PLUGINS=I:/OSGeo4W64/apps/Qt5/plugins
O4W_QT_LIBRARIES=I:/OSGeo4W64/apps/Qt5/lib
O4W_QT_TRANSLATIONS=I:/OSGeo4W64/apps/Qt5/translations
O4W_QT_HEADERS=I:/OSGeo4W64/apps/Qt5/include

Hope it will be solved.


texts: Varying Graduated Symbology Visibility Scale in published service

texts:
One of the things I love about ArcGIS Pro is the ability to vary the graduated symbology of a layer by scale. However, when a map is published to Portal with this functionality, it loses it.
Is there a way to replicate this functionality on a published service?


texts: How to use QGIS3 Authentication Manager to authenticate hundreds of users?

texts: Has anyone used QGIS3 Authentication Manager to handle hundreds of users?  What is the best way to implement role based access for QGIS3?  We are loading layers from PostgreSQL 10.


texts: Spatial Reference and Area/Length Attributes relationship

texts:
I currently allow users to add a shapefile in. From the shapefile, allow the user to generate a new layer base on the polygons from the shapefile's polygon they imported.

Generating the new polygon and putting them into my own layer is no longer a problem, however, I see the Area/Area_Length attribute to be different between my new and old polygon.

Only difference is I projected the Map's spatial reference to the new polygon, instead of using the layer's spatial reference.

Do you guys know if that is the reason behind the different attributes?

Code:

foreach (IPolygon polygon in polygons)
{
   //Need to project to the map spatial reference
   polygon.Project(IMap.SpatialReference);

   //Generate a new polygon
   newPolygon = new Polygon() as IPolygon;

   ((Polygon)newPolygon).AddPointCollection((IPointCollection)polygon);

   //Add the feature
   subcatchmentFeature = new SubcatchmentFeature(newPolygon);
   AddFeature(subcatchmentFeature);
}


texts: NoData cells appeared in segemented raster

texts:
I am trying to classify an aerial image of a dike into grass and weeds. I used an aerial RGB image of 10 cm resolution being a TIF file. This aerial image is cleaned and contains no gaps. Only outside the aerial image raster NoData is present.

I want to perform a object based supervised classification on this image. Before classifying the image, the image was segmented using the Segment Mean Shift tool with the values 19 for spectral detail, 10 for spatial detail and 10 for minimum segment size.

After the segmentation it turned out that the segmented raster has holes in it with NoData cells for the red band at some locations.. These locations are visualised with light blue in the attached figure. The remaining colours are from the segments. The second figure shows approximately the same part of the study area on the RGB image.

Do you have any tips or suggestions why these NoData gaps, shown in light blue, appear in my segmented dataset?


texts: ReconcileVersions tool yields ERROR 000301

texts:
I am wanting to create a version, step into the version, make edits, and then post and reconcile that version.  When I get to the ReconcileVersions_management tool, it is giving me

ERROR 000301: The workspace is of the wrong type

I think that might mean "the version is locked for editing" or "the connection to the db is still open--we can't operate on it".  I'm posting my long-ish code, looking for answers.

import arcpy
import os
from datetime import datetime

# make a connection to SDE.DEFAULT
try:
    filename1 = "fn1"
    foldername1 = os.path.join(os.environ['APPDATA'],"ESRI","Desktop10.2","ArcCatalog")
    servername1 = "gissvc"
    servicename1 = "sde:oracle11g:GISDEVR"
    dbname1 = ""
    auth1 = ""
    un1 = "MyUser"
    pass1 = "MyPass"
    userinfo1 = ""
    ver1 = "SDE.DEFAULT"
    sv1 = "SAVE_VERSION"

    if os.path.isfile(os.path.join(os.environ['APPDATA'],"ESRI","Desktop10.2","ArcCatalog", filename1 + ".sde")):
        os.remove(os.path.join(os.environ['APPDATA'],"ESRI","Desktop10.2","ArcCatalog", filename1 + ".sde"))
    print "creating connection to " + ver1
    arcpy.CreateArcSDEConnectionFile_management(foldername1, filename1, servername1,
                                                    servicename1, dbname1, auth1,
                                                    un1, pass1,  userinfo1, ver1, sv1)
    ws1 = os.path.join(os.environ['APPDATA'],"ESRI","Desktop10.2","ArcCatalog", filename1 + ".sde")
    ver2 = "net_edits"

# create a version

    for version in arcpy.da.ListVersions(ws1):
        if ver2 in version.name:
            print ver2 + " exists, deleting"
            arcpy.DeleteVersion_management(ws1, version.name)
    print "creating new version, " + ver2
    arcpy.CreateVersion_management(ws1, "SDE.DEFAULT", ver2, "PUBLIC")

# connect to new version

    filename2 = "fn2"
    foldername2 = os.path.join(os.environ['APPDATA'],"ESRI","Desktop10.2","ArcCatalog")
    servername2 = "gissvc"
    servicename2 = "sde:oracle11g:GISDEVR"
    dbname2 = ""
    auth2 = ""
    un2 = "MyUser"
    pass2 = "MyPass"
    userinfo2 = ""
    #ver2 = "net_edits"
    sv2 = "SAVE_VERSION"

    if os.path.isfile(os.path.join(os.environ['APPDATA'],"ESRI","Desktop10.2","ArcCatalog", filename2 + ".sde")):
        os.remove(os.path.join(os.environ['APPDATA'],"ESRI","Desktop10.2","ArcCatalog", filename2 + ".sde"))
    print "creating connection to " + ver2
    arcpy.CreateArcSDEConnectionFile_management(foldername2, filename2, servername2,
                                                    servicename2, dbname2, auth2,
                                                    un2, pass2,  userinfo2, ver2, sv2)
    ws2 = os.path.join(os.environ['APPDATA'],"ESRI","Desktop10.2","ArcCatalog", filename2 + ".sde")

# edit some things -- select where status = 5, change to "tacos"

    tbl_vw = "t"
    where = 'STATUS = 5'
    arcpy.MakeTableView_management(os.path.join(ws2, "fooditems"), tbl_vw, where, "", "")
    arcpy.CalculateField_management(tbl_vw, "STATUS", '"tacos"')

# close connections

    arcpy.ClearWorkspaceCache_management()
    print arcpy.GetMessages()

# reconcile with post

    reclogname = "reconcile_log_" + str(datetime.now().strftime('%Y-%m-%d')) + ".txt"
    reclog = os.path.join("F:\PythonScripts","ErrorLogFiles",reclogname)
    if os.path.isfile(reclog):
        os.remove(reclog)
    print "reconciling with post and delete_version options set"
    arcpy.ReconcileVersions_management(ws2, "ALL_VERSIONS", "SDE.Default", ver2,
                                   "LOCK_ACQUIRED", "ABORT_CONFLICTS", "BY_OBJECT",
                                   "FAVOR_TARGET_VERSION", "POST", "DELETE_VERSION", reclog)
except Exception as e:
    print e.message


texts: The proper way to add markers to a layer group in Leaflet

texts:
I always add markers to a Leaflet layer group like this:

someMarker.addTo(someLayerGroup);

But looking at the docs I just realized that you can also do:

someLayerGroup.addLayer(someMarker);

Are there benefits of using one way over the other?


texts: What kinds of line segments/edges require high accuracy in a true surface-of-the-ellipsoid representation?

texts:
I've been musing (and doing prototype coding) for a 'projection free' geographic codebase with your basic point, line & polygon primitives.

Rather than dealing with all the sacrifices that come along with projecting to the plane, however, I'm writing algorithms that work directly on the surface of the ellipsoid.

One of the potential complications is that there are different kinds of "lines" out there:

(arcs of) great circles: shortest distance along the (constant-zero-elevation) surface between two points; should correspond exactly to line-of-sight paths.
rhumb lines: connect the two points with a path of constant direction -- for example, some state borders follow lines of latitude (which are not great circles).
curves: circular arcs (paths of constant distance from a given center-point); Bezier (not sure on correct re-interpretation in the context of a curved surface), etc.

Of the different kinds of paths (including ones I missed), which are important enough that they have an 'exact' representation, vs representing within an error bounds by short segments of a simpler path (e.g. short geodesic arc segments)?

Clarification edits: by 'exact' above, I mean parametric. In other words: computable to any desired accuracy, without a a densification-on-import step.

An edit, much later, to add a citation I've come across that closely parallels my own thoughts on the use of 3D unit vectors as a geographic primitive: A Non-singular Horizontal Position Representation (alt link). Best part? I didn't have to write it all out myself!


texts: Getting OpenLayers pixel value of DEM layer?

texts:
I am working in OpenLayers and need to get pixel value of DEM layer to get the altitude of the place.

Can you guide me?


texts: Trying to build qgis - set path to Bison and flex?

texts:
I'm trying to follow (probably very good) step by step instructions for building QGIS in a Windows environment hence I've got some difficulties to understand where I shall set the paths for Bison and flex in CMAKE and where do I find the option BINDINGS_GLOBAL_INSTALL?

I've uploaded a print screen from my CMAKE

The C compiler identification is unknown

The CXX compiler identification is unknown

CMake Error at CMakeLists.txt:10 (PROJECT):

No CMAKE_C_COMPILER could be found.

CMake Error at CMakeLists.txt:10 (PROJECT):

No CMAKE_CXX_COMPILER could be found.


texts: Does anyone know a good way to handle z ordering when using GraphicTracker in ArcEngine?

texts:
From the GraphicTracker Documentation:

  The first graphic added will appear
  underneath subsequent graphics if they
  have overlapping display areas.

This is not true when using ArcEngine 10.0 with polygon fill symbols on picture marker symbols - sometimes it works (1st time for sure) and sometimes it does not.

Has anyone figured out a proper way to manage z ordering of symbols when using ArcEngine and using GraphicTracker?


texts: Center And Zoom function using ArcGIS API for JavaScript?

texts:
I'm trying to insert x, and y coordinates manually through input text, and center the map and zoom it to that point. I have this code:

var x = $("#x").val();
var y = $("#y").val();
var point=new esri.geometry.Point(x, y, map.spatialReference);
map.centerAndZoom(point,6);

When executing this code, it gives me a blank page. What's wrong?
I also tried this code:

  map.graphics.clear();
  var xMin, yMin, xMax, yMax;
  var mp = new esri.geometry.Point;

  mp.x = $("#x").val();
  mp.y = $("#y").val();
  mp.spatialReference = map.spatialReference;
  var symbol = new esri.symbol.SimpleMarkerSymbol(esri.symbol.SimpleMarkerSymbol.STYLE_CROSS, 30, new esri.symbol.SimpleLineSymbol(esri.symbol.SimpleLineSymbol.STYLE_SOLID, new dojo.Color([255,0,0]), 1), new dojo.Color([0,255,0,0.25]));
  var graphic = new esri.Graphic(mp, symbol);
  xMin = mp.x-1000;
  yMin = mp.y-1000;
  xMax = mp.x+1000;
  yMax = mp.y+1000;

  var newExtent = new esri.geometry.Extent();
  newExtent.xmin = xMin;
  newExtent.ymin = yMin;
  newExtent.xmax = xMax;
  newExtent.ymax = yMax;

  newExtent.spatialReference = map.spatialReference;

  map.setExtent(newExtent);
  map.graphics.add(graphic);

Also does not work. Any help??