texts: Convert lat/long to arbitrary projection?

texts:
I'm writing Windows WPF application using ESRI's axMapControl, and I need to display points on the map, where the points are in lat/lon (WGS 1984), and the map projection is determined by the user.  (User selects which .mxd to display, and I cannot mandate which projection it uses.)

I've been searching around the web, looking for examples of how to project lat/lon, and this is what I've found:

double lon = -3.159875;
double lat = 51.465615;

var srf = Activator.CreateInstance(Type.GetTypeFromProgID("esriGeometry.SpatialReferenceEnvironment")) as ISpatialReferenceFactory3;

// Create Transformation from WGS84 to OSGB86
var geoTrans = srf.CreateGeoTransformation((int)esriSRGeoTransformationType.esriSRGeoTransformation_OSGB1936_To_WGS1984Petrol) as IGeoTransformation;

var wgs84GCS = srf.CreateGeographicCoordinateSystem((int)esriSRGeoCSType.esriSRGeoCS_WGS1984);
var bngPCS = srf.CreateProjectedCoordinateSystem((int)esriSRProjCSType.esriSRProjCS_BritishNationalGrid);

IPoint point = new Point();
point.PutCoords(lon, lat);

var geometry = point as IGeometry5;
geometry.SpatialReference = wgs84GCS;

geometry.ProjectEx(bngPCS, esriTransformDirection.esriTransformReverse, geoTrans, false, 0.0, 0.0);
point = geometry as IPoint;

double x = point.X;
double y = point.Y;

The above works fine for converting from WGS84 to OSGB86. My problem is that it depends upon a particular predefined geographic transformation, and there's no way for me to know which one to use.  I'll know that I'm converting from WGS84, but I'll want to convert to the coordinate system of the map.

In digging through the ESRI docs, I've found reference to IGeometry5 .Project(), which takes only a target spatial reference, which would seem to be what I need, except that it doesn't return the same values as ProjectEx(), in the example above, and the gives me this warning:

  A geometry is not densified before it is projected. This can lead to the output geometries not reflecting the 'true' shape in the new coordinate system. A straight line in one coordinate system is not necessarily a straight line in a different coordinate system. Use IGeometry2::ProjectEx if you want to densify the geometries while they are projected.

So, assuming I want to use ProjectEx(), instead of Project(), is there some way of finding the appropriate geographic transformation, given WGS84, and the spatial reference of a map?

Or is there some other way of projecting a geometry, that does densify, and can work with just the two spatial references, instead of a predefined geographic transform?


texts: Trying to split a complex shapefile by Australian states

texts:
I'm trying to work with a 1:1 million scale shapefile of the geology of Australia.

When the whole lot is loaded in to QGIS (over 200,000 polygons) it will freeze and go to slow to work with.

Ideally I would split the data by each state and just turn on which state I need at a time and that will hopefully run full speed.

However I can't figure out how to split the shapefile properly (700mb file),  I've tried a vector layer split and clip in the geoprocessing tools but I can't get them to work properly. The clip tool just generates a blank shapefile (I tried drawing a polygon shape as a vector file as the output layer so maybe that part was wrong).


texts: Sum line lengths for each grid cell in QGIS

texts: I'm trying to get data from a polyline layer into a grid (I've created this with regular points, points to polylines). I guess at first I have to cut all polyline segments at the grid lines and then I have to spatially join (which only works for the first feature found) all of them into each grid cell. I want to find out road network length differences in different areas.


texts:
I am trying to login to the http://myserver/arcgis/rest/admin page and when I log in with credentials I get no error message but I am returned to the same login page.  If I put in bad credentials I get an error message saying bad login.  I can login to the rest/services page and the ArcGIS Server Manager Site.

I found this but I am not sure if this applies http://forums.esri.com/Thread.asp?c=204&f=2388&t=257171

It looks like I am using a reverse proxy..from looking in the rest.config file under the rest directory.  Does anyone have a clue on what I should investigate?


texts:
Just wondering, how close is the public release of version 2?

I do realise that I can download 1.9, but was wondering about the public version.


texts: Command line tool to repair self-intersection in GeoJSON

texts:
I am working with this file. Kansas seems to have a self intersection:

TopologyException: Input geom 0 is invalid: Self-intersection at or near point -98.918045000000006 37.379078 at -98.918045000000006 37.379078
TopologyException: Input geom 0 is invalid: Self-intersection at or near point -98.918045000000006 37.379078 at -98.918045000000006 37.379078
TopologyException: Input geom 0 is invalid: Self-intersection at or near point -98.918045000000006 37.379078 at -98.918045000000006 37.379078

Is there a command line tool to repair this quickly and easily?

(I do not have access to QGIS.)


texts: Missing R scripts in QGIS 3 after install

texts:
After install Processing R Provider plugin (v.1.0.1) in QGIS 3.4 I only got test.rsx. I'm able to create new scritps but I can't acess scripts already installed in my OS (Linux Mint 19.1) trough R repository or download it (like previous version).

Should I change some configuration in QGIS or OS?


texts: Getting layer name from getFeaturesAtPixel in OL4

texts:
I'm trying to build a popup that only shows feature info from certain layers. How do I get layer name?

I've tried this:

const features = this.map.getFeaturesAtPixel(evt.pixel, {
  layerFilter: ((layerCandidate: layer.Layer) => {
    if (layerCandidate.get('name') === 'my layer') {
      return true;
    }
    return false;
  }
  )});`


texts: Raster calculator issue in QGIS 3.8 and 3.8.1

texts:
I am trying to perform a simple raster calculation on an elevation layer to occlude certain values. I have used the following expression:

( ( "clipped_elev@1">500) AND ("clipped_elev@1"<1000))*"clipped_elev@1"

On QGIS 3.8.0 and after the update to 3.8.1 I get pixels with either 0 or 1 which is not the expected result. I tried the raster calculator from both the Raster menu and from processing. I tried the same expression on the LTR version and got the expected raster result.

Has the syntax for the raster calculator change in 3.8 or is this a bug that needs to be reported?


texts: Infinite length lines coming off polygons in shapefile

texts:
I have been sent some OS MasterMap topography area shapefiles. (OS MasterMap is a national map product including pretty much every permanent structure and land parcel.)

A few of the polygons in each of them have what appear to be lines of infinite length attached to them.

The layers use the British National Grid CRS (EPSG 27700); the lines stretch back to the zero x and y coordinates in the west and south, but carry on infinitely to the north and east.

The three layers in question:

Context of lines against UK:

The lines are not selectable using select feature, and they disappear upon zooming in. GRASS GIS v.clean remove dangles tool removes the vertical ones, but not the horizontal ones.

As they are not selectable, I cannot identify the affected polygons to manually correct them.

I am using QGIS 2.18.0. This problem also occurs on my colleague's PC, with ArcGIS 10. On ArcGIS, the layers cause errors when attempting to get summaries or similar. On QGIS, the summary algorithms simply fail to complete.

What is the cause of these lines and how can I fix them?


texts:
I have a SQLite database, with point and linestring geometries. Between the points(blue) and the linestring(green) is a gap, which I want to close with another line(red). The line should have angle of 90 degres to the existing line and build a connection between the point and the line.
Which are the best SQLite spatial functions to do this?


texts: How to calculate distances between several points along network?

texts:
I have several points (bicycle share hubs), and I want to calculate the distance between each of them using a network. I want a CSV/XLS that has three headers: hub_id1, hub_id2, distance.

I've created a network dataset (after consultation), but how do I get the distance between each of the (blue) hubs pictured below?

I know a dollop of Python, so if this requires some scripting, I'm comfortable trying to apply what I know. I'll be the first to say, though, that I've not scripted with spatial data before.

FYI, both of these posts ask a similar question, but not using a network grid.

Would love some input! Thanks!


texts: How to calculate maximum distance to an interpolated value along a contour in QGIS?

texts: I'm planning to enter concentrations of groundwater contaminants collected from well locations into QGIS and then use interpolation functions to create contours of the groundwater plume.  Can QGIS calculate the maximum distance from a known point (e.g. a well) to an interpolated contour line (e.g. a contour line with concentration of 5 ppb)?


texts: How to use qgis:selectbylocation in PyQGIS?

texts:
I am trying to use the processing toolbox to perform a spatial query (points in polygon). I am not sure how to use the tool and I keep getting Error: Wrong number of parameters.

Here is my script:

polygon_path = r"D:\PythonTesting\SelectByLocation\mitte.shp"
points_path = r"D:\PythonTesting\SelectByLocation\punkte.shp"

polygon = QgsVectorLayer(polygon_path, 'Mitte', 'ogr')
points = QgsVectorLayer(points_path, 'Berlin Punkte', 'ogr')

processing.runalg('qgis:selectbylocation', points, polygon, False, False, False, 0)

I used the documentation, which indicates that I have to use 6 parameters.

However, if I call alghelp the following is returned, indicating 5 parameters:

ALGORITHM: Select by location
    INPUT <ParameterVector>
    INTERSECT <ParameterVector>
    PREDICATE <ParameterGeometryPredicate>
    METHOD <ParameterSelection>
    OUTPUT <OutputVector>

METHOD(Aktuelle Auswahl ändern mit)
    0 - creating new selection
    1 - adding to current selection
    2 - removing from current selection

Either way I am not getting it to run correctly.


texts: How do I use sys.exit(0) in an arcpy script to exit early without having an error message show up?

texts:
I have an arcpy tool script for ArcGIS 10.0 that has two major functional sections.  The user can choose whether or not to run the second section.  If the user chooses NOT to run the second section, I simply want to run a cleanup function and exit the script with a sys.exit(0) without having an error message posted to the tool results window.  There are two major threads here in GIS-SE about exiting arcpy scripts, but the solutions therein don't specifically address the error message.  The general structure of the code is as follows:

import sys
##import arcpy

def CleanUp():
    print 'Cleaning up ...\n'

def finish():
    CleanUp()
    print 'Exiting ...'
    sys.exit(0)

do_more = False  #or True ... input from user

#Section 1:  do some stuff
print 'Doing some stuff ...\n'

if not do_more:
    finish()

#Section 2:  do more stuff
print 'doing more stuff ...'
CleanUp()

If I run this code in the Python interpreter outside of ArcGIS/arcpy, it works as I'd expect, exiting gracefully with no error message; however, in my arcpy script with this same structure, the script exits, but a SystemExit error message is posted to the tool results window.  Is there a way to make the arcpy tool script eat the exception and bury the SystemExit error message?


texts: Add marker to projected coordinates from form-input Leaflet

texts:
I am working on a Leaflet map and trying to use an input form to insert my local coordinates that will be projected to WGS84 and displayed on map as a marker but when I try it the button seems not to work. I don't know if it's a Javascript syntax related issue or I am misusing Proj4.

Here is the code:
HTML
         X: <input type="number" id='latt'>
         Y: <input type="number" id='lngg'>
         <button onclick="addmarker()">addmarker</button>

Javascript
            // convert coordinates to Moroccan system
        var secondProjection = "+proj=lcc +lat_1=29.7 +lat_0=29.7 +lon_0=-5.4 +k_0=0.9996155960000001 +x_0=500000 +y_0=300000 +a=6378249.2 +b=6356515 +towgs84=31,146,47,0,0,0,0 +units=m +no_defs";
        var firstProjection ="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs "
        //add markers based on insertion
        function addmarker() {
            var inputlat = document.getElementById('latt').value;
            var inputlng = document.getElementById('lngg').value;
            var pj = proj4(secondProjection,firstProjection,[inputlat,inputlng]);
            var x = pj[0];
            var y = pj[1];
            var latlng = L.latLng(x,y);
            L.marker(latlng).addTo(map);
        };

EDIT: Type = "number" instead of "float"


texts: How do you test QgsDataSourceURI.setConnection is valid in a QGIS plugin?

texts:
I am currently trying to write my first Python plugin for QGIS and I'm using the the

uri = QgsDataSourceURI()
uri.setConnection(host,port,db,username,password)

syntax to create a connection to an Oracle database. This can then be passed to a QgsVectorLayer(uri.uri(), name, "oracle") command and the layer is added.

The trouble is that before calling this code I am asking the end user for the username and password. The details given are stored to save asking them again.

If they give incorrect details then QGIS will automatically ask them to correct the details but this is outside the scope of my Plugin and so the stored credentials would be wrong.

I want to check the credentials are valid after asking for them but before I try to add the layer.

Is there a way to use my constructed URI in python so that the database will tell me if they are valid?

Something like QgsCredentials but just returning true/false depending on the response of the server.


texts: Managing PostGIS default privileges for users in QGIS?

texts:
I have a PostGIS database for several projects, most of which we work with in QGIS. Each project has it's own schema. For each schema I have created a group role project-members, where project is replaced by the respective project name. What I want to achieve is that all members of a project have read rights for all objects that are created byother members within their project.

For that I thought of the follwing settings:

GRANT CREATE, USAGE ON SCHEMA *project*;

ALTER DEFAULT PRIVILEGES IN SCHEMA *project*
GRANT SELECT ON TABLES TO *project*_members;

Now each user establishes a connection to PostGIS in the QGIS browser and members of the same project should be able to see what the other members did and load layers. Unfortunately, this does not work with my setup. In this answer I learned why: The default privileges are only set for the group role and not for it's members. The proposed solution - to create tables only as the group role and not as individual users - works when I use the console but does not work for users that are connected to the database through QGIS because cannot use SET ROLE ...within QGIS. Also, I don't want to bother every QGIS user with PostgreSQL commands.

Does anybody have a hint, how I can solve this? Maybe I can change my setup in a way that works the way I would like it to.

With help of the answer from DavidP and this answer I could come up with the following solution: Whenever a table is created, the event trigger trg_create_table calls the function trg_create_table_func(). This function queries for all users, which have USAGE and/or CREATE privileges in the schema in which the table was created and grants them the SELECT privilege for the new table.

CREATE OR REPLACE FUNCTION trg_create_table_func()
RETURNS event_trigger
LANGUAGE plpgsql
AS $$
DECLARE
    obj record;
    usr record;
BEGIN
    FOR obj IN
    SELECT *
    FROM pg_event_trigger_ddl_commands()
    WHERE command_tag IN ('CREATE TABLE','CREATE TABLE AS')
    LOOP
        FOR usr IN
        SELECT split_part(unnest(nspacl)::text, '=', 1) AS user_name
        FROM pg_namespace
        WHERE nspname = split_part(obj.object_identity , '.', 1)
        LOOP
            EXECUTE 'GRANT SELECT ON TABLE ' || obj.object_identity || ' TO ' || usr.user_name;
        END LOOP;
    END LOOP;
END;
$$;

CREATE EVENT TRIGGER trg_create_table ON ddl_command_end
WHEN TAG IN ('CREATE TABLE','CREATE TABLE AS')
EXECUTE PROCEDURE trg_create_table_func();

Since this is my first time working with functions and event triggers in PostgreSQL, please feel free to improve if you know a better way.


texts: How to Calculate Utilization Distribution Overlap Index (UDOI) in ArcGIS

texts:
I'm interested in calculating the overlap of home ranges (=utilization distributions in my case) of a bird species using sighting coordinates. This article introduced UDOI as the best index for such purpose:

I have ArcGIS 9.3, and I am able to calculate kernel home ranges using the ArcMap plug-in "Home Range Tools" (HRT). I also have "Hawth's Tools" in my ArcMap. Now I am hoping to find a way to calculate UDOI to demonstrate the degree of overlap between individuals.

I see that UDOI can be calculated in R, yet I do not know how to use R. Can someone tell me if there's any Arc tool that does this? Thanks a lot in advance!

I also see that the user "MTGIS" has asked this same question before:

And @whuber explained the math of UDOI in detail in the wonderful reply. However, I am not super familiar with ArcGIS, and did not understand how to do those steps in Arc. Sorry about repeating the same question! @whuber and anyone familiar with ArcGIS, would you mind walking me through the Arc function that does this?

Thanks!!!


texts: Make a certain layer always appear on top of all the other layers?

texts: I'm finding a lot of material written about how to order features within the same layer, but I can't find anything about how to order the layers themselves in a WMS request. Please note that this is specifically about WMS.


texts: geojson or shapes of ADM3 in all countries in the world

texts:
Where or how can I get the ADM3 geojson or shapes for all the countries in the world?

I am new in Dataviz and I tried geonames.org but I only get points and not polygons.

ADM3 is the 3rd category of country division. So for instance in the case of USA, it is ADM0 = Country, ADM1 = State, ADM2 = County, etc..


texts: ArcGIS Server 10.1 - Cached tiles cut off features

texts:
I have a, ArcGIS MapServer with cached tiles. The problem I ran into is, that the features which are next to the MapServers extent are cut off.
The image below should explain the issue.
The blue line displays the x-min value of the extent. As you can see, only the half of the feature-icon is visible.
How can I tell my MapServer to cache the tiles with all features visible. Is there any way to specify some kind of padding to the extent?


texts: how to get a line that coincident with another line?

texts:
I have 2 linestrings:
Linestring 1 (-3 3, 6 8, 15 9, 22 13, 28.5 15.7) --> red color
Linestring 2 (0 4.66666666666667, 6 8, 10 8.44444444444444)  --> blue color

Actually both two linestrings are coincident

How do I get a part of linestring 1(red) that are not coincident with linestring 2 (blue) in PostGIS (as a multilinestring) ?


texts: Create a path following a road, interpolating from a set of points

texts:
I have several hundred vehicular routes within a city from which I wish to create shapefiles.  These take the form of directions such as "From Green Street, turn left into Brown Street, Right into Blue Street".

I first wrote a program to change them into the form "Corner of Green and Brown; Corner of Brown and Blue" and then geocoded them using the Google geocoding service.  I now have sets of coordinates as csv files, each file (route) has between 20 and 50 points.  I then plotted the route using www.gpsvizualiser.com and where the roads are straight it has worked very well.  However, roads are not straight, and I wish to constrain these paths to follow roads, I'm supposing by interpolating points between the coordinates I have.

I'm using QGis and Open Street Map data (I'd be prepared to try something else if this is not adequate).  Using the RoadGraph plugin I can find routes between two points, but it is very resource intensive, and as I have several hundred sets of points and thousands of pairs of points, I need to find a better way.

The question was asked and answered by a now-inactive user in Shortest Path using interpolated OSM points, however he provided his answer as a link to his blog, which I was disappointed to find is now defunct.

If there is no simple way to do this I will need to do it manually.  If that is the case, is there a simple way to "snap" a line that I draw in QGis to an underlying polyline layer (the OSM data)?  I may need to do some of them manually anyway, as I'm sure the data I have generated is not perfect.


texts:
I am trying to plot a point on top of a shapefile.
the Coordinates of the point are accurate. I verified it using the go to XY tool of ArcMap. but when I try to plot overlay using Python the location of both is distorted. following is my code.
loc =gpd.read_file(r'C:\Users\MJ\Desktop\populationFunc\rentalRtes\resapart\gurgaonGISdata\locGISd.shp')
loc = loc.to_crs(epsg = 32643)
loc.plot(figsize = (10,10))
lat =  28.44198794960
long = 77.07818660830
p1 = Point((long, lat))
df = pd.DataFrame({'a':[lat,long]})
po = gpd.GeoDataFrame(geometry = [p1], crs = loc.crs)
po = po.to_crs(epsg = 32643)
po.plot()
fig,ax = plt.subplots(figsize = (15,7))
loc.plot(ax=ax, color = 'none', edgecolor = 'black')
po.plot(ax= ax)`


texts: Moving marker in a Google Map on Android

texts:
I am trying to create something similar to
http://dev.syskall.com/map/

in my android phone application. Basically I have a List of GeoPoints and I need to make the car/marker move on the map according to the values of the geopoints. Putting the mapView code in a loop does not help at all! The map is invalidated till it reaches the last location in the list and then shows up with the latest location only so I can't see any movement of the marker.

If anyone has created this on Android - please do get back to me! Appreciate your help very much!


texts: Getting coordinates to display properly when imported from excel into ArcMap?

texts:
Both the base map and the table imported have been put to same coordinate system.
Is only showing a single point at a wrong location for all points.

When I click "identify" for the single point it shows all my points with correct lat and long but with different location than the lat and long entered from the table.


texts: Using ArcGIS Desktop to find data folder?

texts:
I am taking a GIS certification course on Courseera. I downloaded the ArcGIS for desktop 10.3.1. I am trying to work with ArcMap. I created a GIS folder on my desktop and followed the steps -given to me by my course instructor, to connect the folder to ArcMap. Once I do that, there is supposed to be an "IntroductionToArcMap" data folder that pops up. The folder does not pop up. I have uploaded the image below of what my screen is supposed to look like. The folder does not pop up.

Has anyone had this same issue?

Is there something that I need to do different or is it the software itself?


texts: Why can't I use SelectLayerByAttribute_management while a cursor is set

texts:
I've been having a problem with arcpy in that the simple line:

arcpy.SelectLayerByAttribute_management ("TEMP_LINES")

produces this response: Runtime error : ERROR 000358: Invalid expression Failed to execute (SelectLayerByAttribute).

I've worked out that this error only occurs when a cursor is set, like this:

lineCursor = arcpy.UpdateCursor("TEMP_LINES","","","combined_i")

Can I not use select by attribute and the cursor or do I have one of them set up wrong?

EDIT: whole code:

arcpy.MakeFeatureLayer_management (Lines,"TEMP_LINES")

lineCursor = arcpy.UpdateCursor("TEMP_LINES","","","combined_i")

for row in lineCursor:
    arcpy.SelectLayerByAttribute_management ("TEMP_LINES")
    arcpy.SelectLayerByLocation_management (Polygons, "INTERSECT", "TEMP_LINES")
    arcpy.FeatureClassToFeatureClass_conversion (Polygons, "temp_polys")
    polyCursor = arcpy.SearchCursor("temp_polys")
    blocks_list = []
    for row in cursor:
        blocks_list.append(row.getValue("BLOCK"))
    blocks_str = ', '.join(blocks_list)
    lineCursor.BLOCK = blocks_str


texts: ArcGIS 10.1 Python AddIn Using Multiprocessing crashes ArcMap?

texts:
I would like to run a multiprocessing task from a python add-in tool.  My issue is that the process keeps failing.  Basically crashes ArcMap.

Here is my basic code:

def function(startOID, endOID, fc):

    wrksp = r"c:\temp\mp_addintest\data\test_%s.txt" % (int(startOID) + int(endOID))
    # real logic removed to dumb it down
    with open(wrksp, 'w') as writer:
        writer.write("%s to %s from %s \n" % (startOID, endOID, fc))
    return wrksp
class btnMP(object):
    """Implementation for src_addin.MPButton (Button)"""
    def __init__(self):
        self.enabled = True
        self.checked = False
    def onClick(self):
        pool = None
        try:
            pythonExe = os.path.join(sys.exec_prefix, 'python.exe')
            multiprocessing.set_executable(pythonExe)
            pool = multiprocessing.Pool(4)
            results = []
            for i in xrange(4):
                results.append(pool.apply_async(function, [str(1),
                                      str(i),
                                      str("test")]))
            pool.close()
            pool.join()
            for result in results:
                print result.get()
        except:
            del pool
            print 'error'

If I run the code outside of ArcMap or from a toolbox, it works without a problem, but when I put the logic inside a button, it causes arcmap to crash.

My guess is that ArcMap is running in process for all python add-ins.  Is there a work around for this issue?

I've tried adding in the freeze_support() to the code as well, but that did nothing as well.


texts: Exporting All Streets in a City - Multiple Lines, With Latitude and Longitude

texts:
I installed Postgresql, Gis, SELECT postgis_full_version() reports,

 POSTGIS="2.2.2 r14797" GEOS="3.5.0-CAPI-1.9.0 r4084" PROJ="Rel.
 4.9.2, 08 September 2015" GDAL="GDAL 1.11.3, released 2015/09/16"
 LIBXML="2.9.3" LIBJSON="0.11.99"

I was able to import data for a single city, from a Geofabrik file, through osm2pgsql.

Now I'd like to export all pedestrian streets in this database in lat,lon,name format. Same street could show up in multiple rows, a street being composed of multiple line segments - that is fine.

I am able to use this query to get streets, and I believe they are line pieces, but how do I get the latitude / longitude associated with the beginning and end of that line?

SELECT name FROM planet_osm_line WHERE highway in ('pedestrian');

What would be the simplest query to do this?


texts:
I am trying to assign attributes from one point layer to another point layer based on closest location.

So I am trying to use a spatial join and generate the ArcPy code in order to ensure reproducibility.

I have tried two ways to do a spatial join.

First, "right click on the layer file in the table of contents" > joins and relates > join > ... . This works but does not generate Python code.

Second, Analysis toolbox > Overlay toolset > spatial join. This generates Python code but also an "Error 999999: An expected Field was not found or could not be retrieved properly. Failed to execute (SpatialJoin)." (the second sentence of the error message appears roughly 20 times). I played around with different "Field Maps", "Match options", "Search Radius" and "Distance Field name" without success.

Running the following python code has also generate the same error:

import arcpy

arcpy.env.workspace = "Path_to_workspace"
arcpy.SpatialJoin_analysis ("dataset1", "dataset2","Path_to_folder_with_output_file/output_file.shp","JOIN_ONE_TO_ONE", "KEEP_ALL", "#", "CLOSEST", "#", "#")

Do you know how I can generate Python code from the first approach or avoid the error in the second code?

Solution

I think the issue arose because "dataset1" and "dataset2" where generated from CSV files using coordinates and the MakeXYEventLayer_management command. So once I saved them as a shapefile and used the shapefile for the spatial join it worked. It is not obvious to me why this is the case, so feel free to add intuition.

I solved the issue using the following code:

import arcpy

arcpy.env.workspace = "Path_to_workspace"

# Define projection
spRef = "path_to_projection/projection.prj"

# Data set 1
arcpy.MakeXYEventLayer_management("path_to_table/table_dataset1.txt", 'x', 'y', 'dataset1',spRef)
arcpy.FeatureClassToFeatureClass_conversion('dataset1', 'Path_to_folder_with_input_file', 'dataset1.shp')
# Data set 2
arcpy.MakeXYEventLayer_management("path_to_table/table_dataset2.txt", 'x', 'y', 'dataset2',spRef)
arcpy.FeatureClassToFeatureClass_conversion('dataset2, 'Path_to_folder_with_input_file', 'dataset2.shp')

arcpy.SpatialJoin_analysis ("Path_to_folder_with_input_file/dataset1.shp", "Path_to_folder_with_input_file/dataset2.shp","Path_to_folder_with_output_file/output_file.shp","JOIN_ONE_TO_ONE", "KEEP_ALL", "#", "CLOSEST", "#", "#")


texts:
Doing a search of GIS.StackExchange, and all of StackOverflow, for the term "Luciad" will bring up zero relevant questions about the Luciad framework. From what I understand, it is comparable in many features to other libraries such as ESRI.

Why is it that no Luciad/LuciadMap questions are being asked or answered on any of these sites? Are there distinct disadvantages to the Luciad framework that prevent widespread adoption, or even discussion, of this framework?


texts:
Mac OSX 10.6.8 QGIS 1.8.0 : I'm mapping an area that's a rectangle slanted at 26 degrees, and neighbourhood maps are always printed so the rectangle is horizontal. I'm trying to make the labels readable when printed that way.

In advanced settings the placement 'around centroid' shows a rotation option, but it's greyed out and unresponsive; none of the other options allow specifying a label orientation. The move/rotate tools of ABC are also greyed out and unresponsive. (The 'free' setting puts labels everywhichway and is very awkward to read.) I've added x, y and rotation attribute fields (type real) and enabled editing.

Might there be a way of enabling the rotation field that I haven't spotted? Or, something else that will have a similar effect?


texts: QGIS loop/iterator with python in function editor (field calculator)

texts:
I am trying to create a loop to rename the fields of my attribute table. The fields should start at 1 and increment by one until the end of the table (e.g. 1,2,3,4,5,...,n)

my code on the python shell works fine, but not in QGIS, why?

Can anyone help me?

def loop (x):
    num = 1
    while num < x:
       print num
       num += 1

loop(188)


texts: I want to use my WMS layer as a background in OSM iD editor. I know it works in JOSM but I couldn't figure it out in iD editor. Does anyone know how to do this?


texts: "Levelling" elevation data coinciding with polygon?

texts:
I have a shapefile containing building geometries (polygons) and a elevation dataset. The elevation data doesn't perfectly match the buildings.

I want to "flatten" the elevation data so that wherever there's a building the elevation will be completely flat.

In short: For all pixels that intersect a polygon, alter their value to the average of the intersecting pixels.

How can I accomplish that?

This Question is more or less unrelated to a specific software, though I prefer using free software like QGIS, GDAL/OGR etc.


texts:
looking at rasters showing change in river channel over time. I have used the raster calculator to calculate the changes e.g. below is the attribute table when subtracting the 2016 raster from the 2015 raster.

-1 indicates areas of erosion whilst +1 shows areas of new floodplain. 0 is just areas of no change. The cell size is 1.2 so how would I convert the count to an area e.g. m^2?. All of the layers are under the BNG spatial reference


texts: Unable to execute WPS GET request in OpenLayers 3

texts:
I am working on a Web GIS development application and using GeoServer 2.11.1, OpenLayers 3 and WPS extension v1.0.0 on Windows 7 and Ubuntu Linux 14.04.

Objective: My objective is to use buffer for spatial analysis (means user will select feature after that user will assign  specific distance say 100m then output will display accordingly) using WPS extension v1.0.0

Development: I am using WPS HTTP GET request to get spatial result. For that I referring following:

http://geoprocessing.info/wpsdoc/1x0ExecuteGETEncoding

Also getting an error as following:

1."Null point exception error". Please check attached screenshot.

Also researching on my side as well as looking  alternative also.


texts: Performing Object Based Image Analysis with NAIP?

texts:
I use NAIP imagery (primarily in QGIS/GRASS) for agricultural analysis.  I want to separate out objects such as grain silos, farms, buildings, &c. using segmentation and classification.

Is this possible with NAIP's somewhat lackluster resolution?

I know with LANDSAT imagery I could pansharpen, but I don't have that luxury with quad-band.


texts: Using PostGIS for a large trajectory dataset

texts:
I have a spatial trajectory dataset which basically contains about 10E6 (lat,lon,time)-tuples, grouped together by an identifier so they form trajectories.

Currently, SQLite is used to store the raw data.

Lat/Lon are in WGS84 but span only small areas (i.e. European country size).

I would like to store it in a geodatabase for easy descriptive & exploratory analysis, eg. calculating the overall and average distance for trajectories, and basic data preprocessing, eg. dismissing all tuples that are outside a polygon, for example.

Thus, I thought PostGIS would be optimal since I have worked with it before.

Some questions arise though: Should I just store the tuples as POINTS with GEOMETRY data type, or should I use MULTIPOINTS for whole trajectories?

And what about performance? Performance is not critical since the data is used for a research project and not for a production application, but I don't want to wait 5 hours to calculate distances between a million point-pairs.

What about indexing? What type of indices should I use and on what type of data?

I know this question is kind of general but I would like to see if people have experience with storing trajectory data and maybe they can give me some advice.


texts: Translating ArcPy Commands to GDAL Equivalents?

texts:
I'm trying to figure out GDAL equivalencies to the following arcpy functions:

arcpy.env.overwriteOutput()
arcpy.MakeFeatureLayer_management()
arcpy.SpatialJoin_analysis()
arcpy.da.SearchCursor()
arcpy.Delete_management()
arcpy.AddMessage()
arcpy.AddError()

I did not see anything in the Python Script examples thread to help me translate the arcpy commands over to GDAL.  If I'm missing something obvious to an experienced GDAL / arcpy user would you be kind enough to point it out to me?


texts:
I've been using ArcGIS collector for field surveys and measurements. I would like the app to automatically populate a field with latest update time when I enter values in the fields I'm interested in.

Is that possible? I'd like to avoid using tracking, which seem to be a possible workaround when googling.


texts: Why does the Semi-Automatic Classification Plugin repeatedly claim that I need to restart qgis?

texts:
I have installed the Semi-Automatic Classification Plugin several times, both through the plugin manager, and via a manual installation.

But it repeatedly tells me I have to restart QGIS? But restarting QGIS does not change a thing.

I read on this site that installing the SEXTANT plugin was necessary, but that did not solve the problem either.

So far I tried all possible steps. I upgraded python to the latest version (3.6.1), upgraded QGIS to the latest version (2.18).

By the way, since I installed the Semi-Automatic Classification Plugin, I get a message that a mini-dump has been written, each time I close QGIS. I suspect this is related to the problem of getting the desired plugin running, but I have no clue what to do with it?

Ultimately I want to download Sentinel-2 imagery in QGIS, any suggestion on alternative plugins therefore are also welcome?


texts: Linking up GPS points taken from rail network

texts:
Bit of a newbie here, trying to work out whether any PostGIS functions can help with my task. Apologies in advance if my use of GIS terminology is not up to scratch yet.

I have a large database of GPS points, ie lat/long coordinate pairs (large being "millions" of points).
Currently they have no relational information, so they are not linked together.

However, they were collected over the UK rail network, typically every few metres. So all points map to rail tracks.

I need to process these points and link them together, so they map onto real railway lines.
I want to identify which GPS points are buffers (ie ends of track) and points (in the rail sense - ie where track splits) as well as those that are "just sections of track".
This needs to be a process that runs automatically with no user intervention.

There are many areas where multiple tracks run in parallel, or areas such as large stations with lots of points & sidings, where there are large numbers of GPS points in a small area.

Being rail, there is a limit to how much the line will curve, so that must be a key feature to aid analysis. As the GPS points are close together, it should be possible to see which form lines that could be track, and dismiss those that turn too sharply (which must be points from parallel tracks).

Initially I had been planning to write some java code that would pick GPS points from the DB that fall within a limited geographical area of each other, and look for GPS points within that region that continued in a straight-ish line. The add/update "previous node id" and "next node id" fields for each point.
Then repeat, using the last linked node as the centre of the next geographical region to analyse.
When the end of a track was reached, then I'd pick the next "untraversed" point in the DB and start again. Then repeat this cycle until every point in the DB had been traversed.

Initially the data is in a mySQL database table.
I'd added an extra field, in which I'd converted the lat/long pairs to a geo "Point" field.
Then used some rudimentary "MBRContains" queries to return points within a certain area (eg 100m2).
I was then going to write some java code that explored/traversed the points, with repeated calls to get points within small areas.

So my question(s) are:

Does PostGIS offers any functions that could assist with the task of identifying related points ?

Or are there any better approaches out there to achieve this ?

Maybe there are tools/routines to do just this kind of thing, but I don't know the terminology enough to know about it yet :)

Thanks for any advice.


texts: Significance of Pgr_trsp() Function returning -1 in id1 column

texts:
I want to understand what's the significance of -1 value in column id1 of pgr_trsp() function output. When I run the below query I was excepting the source and target edge to be snapped according to the floating value I provided as input, So that my final route is snapped precisely to my point of start and end. But pgr_trsp() gave the entire route without snapping. Instead it just had -1 value in column id1. This made me ask what the use of this -1 value. Is there another function we can use to snap the route based on this value. Below is query.

SELECT * FROM pgr_trsp('SELECT feat_id as id, source, target, cost_len as cost FROM test.table_route',
5,
0.5,
600,
0.2,
false, false) AS shortestPath
JOIN test.table_route nav ON shortestPath.id2 = nav.id order by seq;

Here 5 and 600 are the id of my start and end edges. 0.5 and 0.2 are my floating values to snap the start and end edges.


texts: ArcObjects - Interacting with existing ArcGIS Server layer

texts:
I have a bunch of users running a 10.2.1 VB.Net extension, which assumes they have a particular MXD loaded. That MXD contains an existing ArcGIS Server layer (ie. not added programmatically by the extension):

I can get an ILayer reference to that layer using IMap.Layer([idx]), but the returned COM object does not seem to support any other interface. Casting it to IMapServerLayer, IMapServerGroupLayer, IMapServerSubLayer etc all return Nothing.

What interfaces should this object support, and how do I get a more useful object than ILayer out of it? (eg. to be able to perform spatial queries on it?)

Edit: I had a crack at this sample:

Dim uidMapServerLayer As New ESRI.ArcGIS.esriSystem.UID()
uidMapServerLayer.Value = "{E9B56157-7EB7-4DB3-9958-AFBF3B5E1470}" ' IMapServerLayer GUID

Dim enLayers As IEnumLayer = pMap.Layers(uidMapServerLayer, True)
enLayers.Reset()
Dim eLayer As ILayer = enLayers.Next
While Not eLayer Is Nothing
    Dim eMSLLayer As IMapServerLayer = eLayer
    eLayer = enLayers.Next          ' <-- breakpoint here
End While

...and eMSLLayer is a valid object. But if I break it at the breakpoint line and execute CType(eLayer, IMapServerLayer) I get Cannot convert to 'Interface IMapServerLayer'...(this was what I was using to test if the ILayer object supported the various interfaces, so apparently I'm barking up the wrong tree)


texts: Import public online map layer into QGIS

texts:
I am collecting data about mobile phone coverage in my city. This is one of the public maps I found online. In the source code of the website owner I found the following ressource which is the visual layer of the LTE / 4G coverage:

https://netmap.vodafone.de/arcgis/rest/services/CoKart/netzabdeckung_mobilfunk/MapServer/export?dpi=300&transparent=true&format=png8&LANGUAGE=ger&layers=show%3A0%2C100%2C1&bbox=764766.9499904191%2C6664654.659215047%2C838146.497144092%2C6686018.808620726&bboxSR=3857&imageSR=3857&size=1920%2C800&f=image

I can play around with those parameters, I can call the image directly. But I can't get it into my QGIS as a layer.

Question: how do I have to set the "New WMS layer" source? Current approach is this which does not work.
https://netmap.vodafone.de/arcgis/rest/services/CoKart/netzabdeckung_mobilfunk/MapServer/

The QGIS version I am using is 2.14.21.


texts:
I want to bind ZoomBox function to a toolbar button to active and deactivate, its default activator is Shift key. Is there a way to activate without pressing the shift key?

In OpenLayers 2 I could add the ZoomBox function using:

self.zb = new OpenLayers.Control.ZoomBox({
        title: "Zoom box: zoom clicking and dragging",
        text: "Zoom"
});
self.map.addControl(self.zb);

and to active:

this.zb.activate();

I can't find the ZoomBox objects under ol.control class in Openlayers 4. Is there a ZoomBox equivalent or alternative in OpenLayers 4. If not how can I do this programmatically?


texts:
I am trying to do an unsupervised classification using GRASS GIS software. I have hyperspectral photo and now analyzing vegetation indices. Here is what I did:

Created 2000+ raster photos with different combinations of bands
Created raster of signatures for the same area using the same mask
Used i.group command to make group of all rasters that I have
Used i.cluster command - and here is my problem:

  WARNING: quantization file [gndvi_g40_n75] in mapset [xxx]
  missing
  WARNING: Unable to open header file for raster map
  WARNING: Unable to open raster map
  ERROR: Unable to open raster map

What should I do to not have those errors? It looks like i.group did not create everything properly...

So I cannot use i.maxlik without i.cluster


texts: How to create buffers around points in ArcMap when there is a lack of OID

texts:
I want to create buffers around points in ArcMap. Using the Buffer (Analyst) tool, I indicated 100m buffers around my points of interest.

I get the error:

ERROR 000339: Input price, ads (v8) - school.csv Events does not have OIDs
Failed to execute (Buffer).
Failed at Wed Nov 07 16:49:00 2018 (Elapsed Time: 0.01 seconds)

The points are in a CSV file. How do I correct this error in order to generate the buffers?


texts: GDAL python package with Anaconda and ElasticSearch 7

texts:
I'm trying to import feature data into elasticsearch v7.4 using GDAL to convert from an ArcGIS feature service layer into geoJSON for indexing.

We use anaconda 3 (latest version 2019.10) for our python environment. When you pull GDAL it seems to be pinned back at version 2.3.3 - despite being up to version  3.0.x on conda-forge.

Unfortunately only version >= 2.4 of GDAL works with ElasticSearch > 2.

I've tried switching to the conda-forge channel but incompatibilities make my environment unresolvable.

I've tried using pip instead but it's getting complicated as it requires MS Build tools as it's trying to build from source versions.

I've managed to make it work as a proof of concept using a GDAL docker image but I need to productionise this with Anaconda.

Does anyone have a workaround for this?


texts: Leaflet: Error when showing marker name within an info box/popup from GeoJSON data

texts:
I am trying to make it so when you click on the marker you get a popup with the name of the station within wales.
Live working of my code:
https://ces-web2.southwales.ac.uk/students/18018815/mitchtut/giscw2/leaflet.htm
I am currently getting this error:

leaflet.htm:167 Uncaught TypeError: Cannot read property 'properties' of undefined
at i.clickOnMarker1 (leaflet.htm:167)

code:


texts: MapInfo legend editing for right to left text

texts:
I work with some right to left text (Arabic, Hebrew) and I would like the legend to be from right to left instead of the current left to right.
Is there any way of doing this in MapInfo?


texts: Use of linear referencing on a route that span over 3 different UTM zones

texts:
I have a long route (~2400 KM) that span 3 different UTM zones. (N - 37, 38, 39)

On this route there is a sign every 1 km-approximately (km marker), these signs help in determining the exact location of assets along the road etc…

I have the actual location of these signs as a point and it turned out that the distance between 2 markers is between 1100 m to 1200 m in general so there is an error in these markers location caused by human mistakes. So I need to create new km markers that conserve the distance between the 2 markers (1000 m) using linear referencing and event tables on my road route.

But the problem is my route passes by 3 different UTM zones and all of my data are geographic WGS84 which is in degree, therefore, I can’t create km markers using linear referencing and event table in this case for each 1 km as my m measurements if the route was WGS 84.

And I need the rout to be 1 single record along the 3 UTM zones so I can’t divide the route based on the projection boundary. Even though I tried to do this and also had an error equal to 20 m on the boundaries of the projection zone.

What should my route projection be in this case? So that the linear referencing and event tables work without any distortion.

For sure the route should be in a projected coordinate system since we're using meters (1000 m) to dived the route or you think we can use the geographic coordinate system by transforming the m value to a decimal degree.


texts: Performing ellipsoid to spheroid conversion on a DTM using gdal

texts:
I have a DTM that is referenced to an ellipsoid (equatorial radius is not equal to the polar radius) and I want to convert the elevation values to the spheroid (equatorial radius = polar radius).
How can this be done using Gdal?
Additional info: The DTM is of a different planetary body


texts: Unable to use pgRouting Layer QGIS plugin with osm2po tool (database connection error)

texts:
By following these amazing tutorials of Anita Graser An osm2po Quick Start and pgRouting on OSM using osm2po, I was able to use osm2po tool successfully on Windows 7 (x64) with PostgreSQL version 9.5 and pgRouting version 2.1.0. I successfully generated .sql file (at_2po_4pgr.sql), imported the same into my postgres routing DB and then installed QGIS pgRouting Layer plugin through Plugins Manager. Now, I am trying to test pgr_djikstra algorithm using pgRouting Layer plugin but I am running into problem during this last step of the above mentioned tutorial.

I am already connected to my postgres DB through DB Manager of QGIS but in my pgRouting Layer dialogue box, I am never "shown as connected" to my postgres DB as shown in following figure. I tried to press refresh button and drop down list many times but no use.

The problem seems to be with database connection through Python. Because, whenever I try to run 'dijkstra' algorithm (just wanted to give it a try anyways) I am getting a Python execution error as shown below. It shows db connection problem in pgRouting python executable file (which is in pgRouting Layer plugin directory). Following are the screen shots of errors I am getting.

Although, Python2.7 is the default version and its already added to my system path (in environment variables).

Do I need to modify the indicated python file?


texts: Reading raster attribute table using gdal in python

texts:
Is there a GDAL function that can read in a raster attribute table? Something functionally similar to this function in pygeoprocessing:

pygeoprocessing.geoprocessing.get_rat_as_dictionary_uri(dataset_uri)
Returns the RAT of the first band of dataset as a dictionary.

Parameters: dataset – a GDAL dataset that has a RAT associated with the first band
Returns:    value – a 2D dictionary where the first key is the column
name and second is the row number
Return type:    dictionary


texts:
I try to reproject a layer in PosGIS using statement:

UPDATE <table> SET geom=ST_Transform(geom, 3857);

leading to error message:

ERROR: transform: couldn't project point (574415 6.22508e+006 0): latitude or longitude exceeded limits (-14)
SQL state: XX000

I thought the first number in the error message corresponds with the id (line number) in the table, but when I retrieve ST_X and ST_Y from that point then the coordinates do not correspond to the second and third number in the error message.

The projection so far is 4326 and a visual control in ArcMap doesn't show any unexpected points.

Can you tell me on how to read the error message and to overcome the issue?


texts: Calculation of TOA radiance from Sentinel-2 1C product reflectance

texts:
I was trying to convert the Sentinel-2 1C product to radiance and came across two links first once says

radiance = ((pixelValueBandX * cos(incidenceAngle) * solarIrradianceBandX) / (pi * d2)) / 10000

where d2 is 1.0/U

Where as second one says

rad = rho * cos(radians(sza)) * Es * sc / (pi * d2)
where: d2 = 1.0 / U
scale: 1 / (0.001 * 1000) = 1 (default)

So, I have few questions amongst the two which one is correct. I don't see any scaling factor in the first one. Moreover, does the incidence angle here denotes the value of Mean_Sun_Angle.ZENITH_ANGLE from file MTD_TL.xml?

Also, Let's say i do atmospheric correction using SEN2COR or DOS. I would get the BOA reflectance, right? So, then how do i convert it to a BOA radiance?


texts: I want to load the GeoTrace plugin into QGIS. I have donwloaded the ZIP file, unzippped it and created a folder C:\Users\UserName\.qgis2\python\plugins\GeoTace and copied the files in there. I still can't see the plugin if I click on plugins in the GUI. Can anybody tell me what I am doing wrong?


texts: How to save qml style including function in QGIS?

texts:
I made a style which needs function to work. In my particular case it is calculating azimuth between nodes of rectangles to display "direction".

The problem is, if I save the style, the function is not saved with it. So if somebody else wants to look at it on other computer, the function is missing and style doesn't work as expected. Is there some way to include function or workaround to not need to use function to achieve saving the style including the calculation of direction?

Poor workaround like create new column which will be populated with this function is not an option here because the key is interactive response of style on change / addition of rectangle.


texts: Programmatically alter AliasName for fields defined in layer file (not underlying source)

texts:
Is there a way to alter the AliasName on a layer file only, and not the underlying shapefile or feature class? I have a few hundred layers that define joins between hundreds of tables with census data sitting in an MS Access database, and a shapefile with census tracts. I like to keep it this way to not have to store the feature geography for each of those data files.

Each layer file maps out hundreds of variables with intuitive names such as B19325_089. All said, there are several thousand. I would like to change the NameAlias property to have meaningful aliases/variable labels based on a lookup value, but, can it be done/defined in the layer file and not in the actual geography file, as all the layer files create the join to the same shapefile?

I got off to a decent start on feature classes in a file GDB using Mark Cederholm's and Matt Wilkie's suggestions, but am, as noted, curious about whether it is possible to have the alias definition reside in the layer files.

Initially, I thought it very simple and doable with the geoprocessor object as ESRI lists the AliasName property as read/write, but that appears to not be possible, and arcobjects to be necessary.


texts:
I need help with writing a script:

The task

I have a point feature class of graffiti incidents and a polygon feature class of patrol zones with some empty attributes already created for you. I must write a script that updates the attributes of the patrol zones with:

The number of graffiti incidents falling within the patrol zone. This is an integer that goes in the INCIDENTS field.
The priority ranking for the patrol zone. This is a string that goes in the PRIORITY field. You will derive this string using some simple math that compares the number of incidents in the zone with the area of the zone.

Patrol zone priority rankings

I then will calculate a priority ranking for each zone by dividing the number of graffiti incidents in the zone by the area of the zone. My script should then examine the result and assign the appropriate priority ranking (PRIORITY). These are the priority rankings:

TOP CONCERN—15 or more incidents per square mile
HIGH CONCERN— At least 12 but less than 15 incidents per square mile
SOME CONCERN— At least 6 but less than 12 incidents per square mile
LOW CONCERN—Fewer than 6 incidents per square mile

Here's what I have so far:

import arcpy
arcpy.env.overwriteOutput= True

# layer that we are changing
patrolZone= "C:\\...\\PoliceData.gdb\\PatrolZones"
# the layer we are selecting from
graffiti= "C:\\...\\PoliceData.gdb\\GraffitiIncidents"
nameField= "NAME"
graffitiField= "OBJECTID"
incidentsField= "INCIDENTS"

# sits above the first row, starts the cursor
patrolRows= arcpy.UpdateCursor(patrolZone)

# tells it to go to the first row
patrol= patrolRows.next()
# now were in the first row doing looping

while patrol:
arcpy.MakeFeatureLayer_management(graffiti, "GraffitiLayer")

zones= patrol.getValue(nameField)
#print incidents
queryString = '"' + str(nameField) + '" = ' + "'" + str(zones) + "'"
print str(queryString)

arcpy.MakeFeatureLayer_management(patrolZone, "PatrolLayer", queryString)
arcpy.SelectLayerByLocation_management("GraffitiLayer", "CONTAINED_BY", "PatrolLayer")

numGraffiti= arcpy.GetCount_management("GraffitiLayer")
print numGraffiti

patrol= patrolRows.next()

I'm stuck on updating the INCIDENTS field


texts: Layout show label and each polygon make serrated image

texts:
I want to develop "Print lay out functionality", but in output system is giving me whole image, I want to print only selected polygon. Also it is not displaying label names on print image.

Please suggest some way.

from PyQt4.QtGui import *
from PyQt4.QtCore import *
import time
l = iface.activeLayer()

for a in l.getFeatures():
     attrs = {16: 1}
     l.dataProvider().changeAttributeValues({a.id(): attrs})
     img = QImage(QSize(1280, 1024), QImage.Format_ARGB32_Premultiplied)
     # set image's background color
     color = QColor(255, 255, 255)
     img.fill(color.rgb())
     # create painter
     p = QPainter()
     p.begin(img)
     p.setRenderHint(QPainter.Antialiasing)
     render = QgsMapRenderer()
     palyr = QgsPalLayerSettings()
     palyr.readFromLayer(l)
     lst = [l.id()]
     render.setLayerSet(lst)
     render.setLabelingEngine(QgsPalLabeling())
     # set extent
     rect = QgsRectangle(render.fullExtent())
     rect.scale(1.1)
     render.setExtent(rect)
     render.setScale(1)
     render.setOutputSize(img.size(), img.logicalDpiX())
     render.render(p)
     p.end()
     # save image
     img.save("C:/example/Ashwin/%s.jpg" % a['CHALTA_NO'],"jpg")
     attrs = {16: 0}
     l.dataProvider().changeAttributeValues({a.id(): attrs})
     time.sleep(1)


texts: Raster pixel size in QGIS automatically defaults to square

texts:
I am using QGIS 3.4 and my data consists of rasters with irregular cell sizes. In ArcGIS Pro, the cell sizes are maintained (screenshot below).

However, when I upload the same raster into QGIS, the cell size gets automatically changed to square, with one dimension of the pixel size being negative (screenshot below).

I think QGIS is automatically resampling the data and formatting to square cells, which I want to avoid. Is this a functionality that I can edit somewhere within QGIS?


texts: Strange NetCDF file with Lat / Long as 2D dimensions?

texts:
I have this NetCDF file from Germany DWD: https://opendata.dwd.de/climate_environment/CDC/grids_germany/daily/Project_TRY/air_temperature_mean/TT_201001_daymean.nc.gz
However, I don't know why this netCDF has Lat and Lon axes in 2D.
double lon(y, x) ;
double lat(y, x) ;

The result from ncdump is below:
ncdump -h TT_20100101_daymean.nc
netcdf TT_20100101_daymean {
dimensions:
    x = 720 ;
    y = 938 ;
    time = UNLIMITED ; // (31 currently)
    bnds = 2 ;
variables:
    double lon(y, x) ;
        lon:standard_name = "longitude" ;
        lon:long_name = "longitude coordinate" ;
        lon:units = "degrees_east" ;
        lon:_CoordinateAxisType = "Lon" ;
    double lat(y, x) ;
        lat:standard_name = "latitude" ;
        lat:long_name = "latitude coordinate" ;
        lat:units = "degrees_north" ;
        lat:_CoordinateAxisType = "Lat" ;
    double time(time) ;
        time:standard_name = "time" ;
        time:long_name = "time" ;
        time:bounds = "time_bnds" ;
        time:units = "hours since 2010-1-1 00:00:00" ;
        time:calendar = "proleptic_gregorian" ;
        time:axis = "T" ;
    double time_bnds(time, bnds) ;
    double datum(time) ;
        datum:long_name = "Date and time in UTC" ;
        datum:units = "YYYYMMDDHH" ;
        datum:_FillValue = 9999. ;
        datum:missing_value = 9999. ;
    short temperature(time, y, x) ;
        temperature:standard_name = "air_temperature" ;
        temperature:long_name = "mean temperature at 2 m height" ;
        temperature:units = "degC" ;
        temperature:coordinates = "lat lon" ;
        temperature:add_offset = 0.f ;
        temperature:scale_factor = 0.1f ;
        temperature:_FillValue = 9999s ;
        temperature:missing_value = 9999s ;

// global attributes:
        :CDI = "Climate Data Interface version 1.7.0 (http://mpimet.mpg.de/cdi)" ;
        :Conventions = "CF-1.4" ;
        :source = "surface and satellite observations, cosmo_090213_4.8_clm17" ;
        :institution = "Deutscher Wetterdienst" ;
        :title = "Temperature daily gridded dataset" ;
        :project_id = "TRY-advancement" ;
        :realization = "v1.0" ;
        :contact = "Stefan Kraehenmann, stefan.kraehenmann@dwd.de" ;
        :creation_date = "2016-01-23 06:07:57" ;
        :CDO = "Climate Data Operators version 1.7.0 (http://mpimet.mpg.de/cdo)" ;
        :history = "Wed Jun 22 10:24:12 2016: ncatted -a history,global,d,c, TT_201001_daymean.nc temp.nc" ;
}

Typically, I see the netCDF file has Lat and Long axes in 1D, e.g:
dimensions:
    lat = 360 ;
    lon = 720 ;
variables:

    float lat(lat) ;
        lat:units = "degrees_north" ;
        lat:long_name = "latitude" ;
    float lon(lon) ;
        lon:units = "degrees_east" ;
        lon:long_name = "longitude" ;

Anyone knows why the netCDF file Lat and Long axes in 2D?


texts:
I am used to digitizing stuff using JOSM, the Java OpenStreetMap editor. Recently, I tried some digitizing in QGIS and it felt rather clumsy.

What I was missing specifically was an equivalent to the "w"-tool in JOSM ("improve way accuracy, https://josm.openstreetmap.de/wiki/Help/Action/ImproveWayAccuracy),  which lets you add and move nodes in linestrings/polygons with single clicks. In my QGIS (3.4) to add nodes to a polygon I needed to click on the cross between two existing nodes and then on the place where I wanted the new node to be. The JOSM tool would have allowed me to ctrl-click to add new nodes. This saves one click and mouse movement per new node.

Are there plugins that allow for faster drawing and modification of lines and polygons?


texts:
I'm trying to create a web application using GeoServer.

When I try to see the file var/lib/tomcat8/webapps/geoserver/data/www/ol-demo.html, no map appears, just  the title "GeoServer OpenLayers Demo" the frame, the zoom in/out.

I have installed GeoServer 2.8.1.

When I make a selection of a layer, I can see them through GeoServer in OpenLayers, but when typing http://localcost:8080/geoserver/www/ol-demo.html or ..../mytest.html, I see the page (no error 404) but I cannot see any map.

Can you help me?


texts: Insert new rows in table using QGIS and PostGIS view and do instead trigger

texts:
I have a table with roads and use a view to populate descriptive fields from look-up tables.
I can see the view in QGIS and would like to use it to add new roads using a DO INSTEAD trigger.
Currently the trigger adds one row with the new data for every row in the view and not just for the one or two new rows ( ie if the view has 2000 rows I get 2000 new rows in the table that all look like the new row)
How do I get the trigger to do the same as the insert on the table which correctly adds a new id (serial) for each new row without manually adding the id?


texts: Extract elevation data from contour lines

texts: I am working with QGIS. I used raster-extract-contour to get contour lines from a raster image. Now I want to break up my single vector file into multiple layers based on elevation. How do I do this?


texts: Import OSM file into pgrouting, keeping certain key/values

texts:
i want to import a big (10G) OSM network into PostGIS for pgRouting. As i have to do some postprocessing with Postgres i need to load some key/values from the OSM dataset into the database as well, e.g. maxspeed=.
What's the best way for that?
Right now i tried the following ways:

Import with osm2pgrouting. Failed, didn't even import the network, most likely due to the big network file.
Import with osm2po: Works, but it seems really hard/impossible to keep certain key/values from OSM.
Import the network with osm2po AND import the .OSM file a second time with osm2pgsql, subsequently merging the tables. Takes hours as it seems impossible to use osm2psql without importing the geometries (which i don't need obviously).


texts: How to install Shapely, but got stuck into Python language

texts:
I have a MAC OSX 10.7.5 and am trying to learn how to use QGIS.
Following a couple of tutorials I came to know I have to install shapely.

I followed the instructions by the multi-cited http://tumblr.pauladamsmith.com, but I am afraid I did not understand neither how to start the installation.

I downloaded all the required files (i.e. GDAL 1.10 and its python package but did not move them in any particular directory - i.e. the complete GDAL framework seems to be installed into Library and the GDAL python package is still in the Download folder).

I did the following, thanks to the instructions of the link above:

MacBook-Pro-di-myname:Downloads myname$ cd /Users/myname/Downloads/GDAL-1.9.1
MacBook-Pro-di-myname:GDAL-1.9.1 myname$ CFLAGS=`/Library/Frameworks/GDAL.framework/Versions/1.10/Programs/gdal-config --cflags`
MacBook-Pro-di-myname:GDAL-1.9.1 myname$ LDFLAGS=`/Library/Frameworks/GDAL.framework/Versions/1.10/Programs/gdal-config --libs` python setup.py install

but then that's what I get:

running install
error: can't create or remove files in install directory

The following error occurred while trying to add or remove files in the
installation directory:

[Errno 13] Permission denied: '/Library/Python/2.7/site-packages/test-easy-install-3861.write-test'

The installation directory you specified (via --install-dir, --prefix, or the distutils default setting) was:

/Library/Python/2.7/site-packages/

Perhaps your account does not have write access to this directory?  If the
installation directory is a system-owned directory, you may need to sign in
as the administrator or "root" account.  If you do not have administrative
access to this machine, you may wish to choose a different installation
directory, preferably one that is listed in your PYTHONPATH environment
variable.

For information on other options, you may wish to consult the
documentation at:

http://peak.telecommunity.com/EasyInstall.html

Please make the appropriate changes for your system and try again.

Where is this error coming from? What I am getting wrong? Since this error impedes me   to continue I cannot follow the other instructions at the link above, so I am pretty stuck   with my work and cannot continue.

I will appreciate any type of help: I am totally ignorant on Python language and programming as a whole, indeed it was an unexpected problem during my QGIS learning.

Thanks a lot
Fra


texts: Getting the nearest point while ignoring polygon in OpenLayers 6

texts:
I have a map that contains one single layer with multiple feature. Those feature contains polygon and point.
I do not filter them because I did load every layout in the same object like the following
geoJsonObject = MyGeoJsonFile

const vectorSource: VectorSource = new VectorSource({
  features: new GeoJSON().readFeatures(geojson),
})

this.vectorLayer = new VectorLayer({
  map: this._map, // Already generated map
  source: vectorSource,
})

Now I try to find the nearest point and move to it while the user is moving
_geolocationService.onUserPositionChanged$ give me the user coordinates {lng: string, lat: string}
moveToCoordinate(coordinates: any[], zoom: number = 16) {
  this.view.animate({ center: olProj.fromLonLat(coordinates), zoom: zoom })
}

this._geolocationService.onUserPositionChanged$
  .pipe(takeUntil(this._unsubscribeAll))
  .subscribe((coordinates) => {
    if (coordinates) {
      const formattedCoordinates = [coordinates.lng, coordinates.lat]
      const fromLonLatCoordinates = olProj.fromLonLat(formattedCoordinates)

      this.positionFeature.setGeometry(coordinates ? new Point(fromLonLatCoordinates) : null)

      const nearestPoint: Feature = this.vectorLayer
        ?.getSource()
        .getClosestFeatureToCoordinate(fromLonLatCoordinates)

      if (nearestPoint) {
        this.moveToCoordinate(nearestPoint.getGeometry().getCoordinates(), 18)
      }
    }
  })

But nearestPoint will give me back the coordinate of the polyline instead of the point.
I know that getClosestFeatureToCoordinate() does accept a filter but I cannot figure out how to make it work.


texts: ggmap to get travel distance from dataframe

texts:
I am trying to use the mapdist() function from the ggmap package to get travel distances from coordinates in a dataframe.

# need an api key to use the ggmap function

register_google(key = "")

library (tidyverse)
library(ggmaps)
starting_point <- tribble(
  ~row_id, ~X, ~Y,
  1, -119.13968, 50.53766,
  2, -134.26792, 60.45440
)
destination_1 <- tribble(
  ~row_id, ~X, ~Y,
  1, -119.13968, 50.53766,
  2, -134.26792, 60.45440
)
destination_2 <- tribble(
  ~row_id, ~X, ~Y,
  1, -119.13835, 50.55011,
  2, -135.05166, 60.71958
)

I like the output produced here

# get coordinates
(wh <- as.numeric(geocode("white house")))
(lm <- as.numeric(geocode("lincoln memorial washington dc")))
# entering coordinate instead of address
(df <- mapdist(wh, lm, mode = "driving"))

I want to get two data frames back, using the starting_point as the start and the two destination data sets as the end. The row id gives to specific locations I wish to target.

I've tried entering the coordinates as a list and matrix and still seem to get an error.


texts: GRASS within QGIS - v.buffer.distance adding additional parameters

texts:
My end aim is to perform a buffer without end caps (-c) and straight corners (-s) based on a distance field. Neither of these options are available within the module options in the GUI within QGIS, is there any way to add these anywhere within the GUI?

I have also tried running GRASS within the regular processing toolbox but keep getting this error


texts:
In QGIS 3.8 I would like to export the code/structure of my data-forms.

I have created them via the "drag-and-drop"-surface.

Does anybody knows a solution to export it (as an xml- or ui-file)?


texts: Create vector grid from canvas using PyQGIS?

texts:
I've modified "Hex grid from layer bounds" script to my need. Just added cellsizex and cellsizey. Changed the hexagonal - 3 to "1" to make it as a square grid. Is it possible to include an option to select extent from canvas? Instead of getting it from layer extent, I want to manually select the extents from the canvas by drawing a square or rectangle.

I am new to scripting and though I can understand what a line of code does, I find it difficult to implement it by myself. For eg, with the help of the above said code, I was trying to get the extents from the canvas through a square or rectangle drawing.

What I did was, Create new script -> just started to type (after the initial declaration or whatever) ## get extent, ##select extent...things like that.

The problem is, i don't know which will be the right method, what is for what?

Is script editor's code different from python's? (Don't know python either).

What do I have to do to build small scripting tools by myself?


texts: Setting tool-specific parameters within the QGIS graphical modeler

texts:
Is it possible to set certain parameters of tools within set up models when opening a saved model (not edit, but double click and running)

I am talking about being able to set certain tool specific parameters such as distance in the "fixed distance buffer" tool as model parameters that then can be changed in the model window before running the model?

I treid to set up strings but the the tools withion the models don´t allow for triggering the string variables. So far only in- and output layers can be specified but not the input fields of desired tool-specific parameters


texts: Obtaining a handle for newly created map in GeoExt

texts:
Hopefully a straightforward question...

I am following the example based on this page here http://www.geoext.org/lib/GeoExt/widgets/MapPanel.html to create a map within a window, here's my code:

new Ext.Window({
    title: "GeoExt MapPanel Window",
    height: 400,
    width: 600,
    layout: "fit",
    items: [{
        xtype: "gx_mappanel",
        id: "mappanel",
        layers: [basemap],
        extent: "-5,35,15,55"
    }]
}).show();

So the above example create a new map in a floating window which is exactly what I'm after...  My question is how exactly do I get a handle to the map within this new window so I can add controls and layers after the map has been created?  Something along the lines of:

map = grab_the_map_somehow_from_Ext.Window_above();
map.addControl(my_new_control);
map.addLayers([my_new_layer]);

Thanks for your help.


texts: How to conditionnally assign a new value to pixels of a raster image?

texts:
Given a topographic GIS raster of one country crop.tif:

Given a know pixel value such as elevation (z) is +73 :

$ gdallocationinfo crop.tif 1 1 -valonly
> 73

Given an elevation threshold n = 50 (meters)

How to set the value of all/this pixels where z >= n to 50 ?

The solution should be with gdal or in console.

Related to : Gdal: How to get a pixel's value from raster image?, gdal_polygonize.py.

EDIT: I eventually used the formula --calc="50*(A>50)" --NoDataValue=0 standing for "new calc where value set as 50 when value in input A > 50, else value set as 0"

  gdal_calc.py -A input.tif --outfile=result.tif --calc="50*(A>50)" --NoDataValue=0


texts: How Can I edit the Geometry of Points when many of them lie at the same place?

texts:
I had some location data as CSV, and I have converted them to a Point Featureclass in ArcMap. Some of those records had 0,0 for lat,long, and hence the data went to Null Island.

I have now got the correct Lat-Long for some of those records, and I need to update the fields as well as geometry for those records.

When I start editing, and use the Edit tool, all of them get selected.

I can't find a way to select individual records, and edit their geometry.

How can I do this?


texts: PostGIS: problem with point count result when they are over the polygon edge

texts:
I have a PostgreSQL 12 database with PostGIS (using EPSG 4326). There are some tables with millions of points, and others with polygons.
In general , I had no problems with the queries. However, I recently noticed that in some cases, when I count points in polygons, some points are excluded or counted twice.
I use different tables with a grid of square polygons of distinct sizes, one degree, half degree, etc.
If I use:
select polygon, count(point) from point, polygon where st_within(point, polygon) group by 1

Some points located over the edge of the polygon are excluded from the count.
If I use:
select polygon, count(point) from point, polygon where st_intersects(point, polygon) group by 1

Those points get counted twice, once for each polygon.
How can I handle this situation to avoid errors in the analysis?


texts: Unable to delete folder with shutil.rmtree() with ArcPy?

texts:
Unable to delete folder with shutil.rmtree() with arcpy:

del mxd
delete_path_edit = temp_dir_extended
delete_path_forward = '/'.join(delete_path_edit.split('\\'))
shutil.rmtree(delete_path_forward)

ERROR

Traceback (most recent call last):
  File "D:\PythonScripts\exportToPDFFolderV1.py", line 166, in <module>
    shutil.rmtree(delete_path_forward)
  File "C:\Python27\ArcGIS10.3\Lib\shutil.py", line 252, in rmtree
    onerror(os.remove, fullname, sys.exc_info())
  File "C:\Python27\ArcGIS10.3\Lib\shutil.py", line 250, in rmtree
    os.remove(fullname)
WindowsError: [Error 32] The process cannot access the file because it is being used by another process: u'D:/MaharashtraIDSI/TempFolder/temp_custom\\Clip_IDSI_SA_2006_001_masked.img'

Failed to execute (ExportToPDFFolder).


texts: What language (or languages) are the next 'big thing'?

texts:
Having (mostly) learned how to navigate Python programming in GIS, I've gotten to thinking about where GIS programming is going from here.  I know VB used to be the go-to language for ArcGIS, but that has been supplanted by Python.  I love Python's utility, ease of use, and support in a variety of spatial and non-spatial programs,  but I know these are no guarantee of it sticking around.

I'm curious about which languages I should keep an eye on and consider dabbling in to be ahead of the curve.  I've seen a few options out there.  For example Geoscript has libraries for Scala and Groovy in addition to js and Python.

What up-and-coming programming languages are trending towards increased use, increasing geospatial functionality, and cross-application scripting?


texts: Spatial analysis and clustering near features

texts:
Im working on marking behaviour in group living animals and im interested in how marking behaviour is affected by certain characteristics of the neighbouring groups. I have plotted the territories of each group from 95% density isopleths created using a-LoCoH. The marks of each individual from a group are plotted as a separate points layer. Marking behaviour does not appear to be random, and they do seem to cluster around boundaries with neighbouring groups. What I would like to investigate whether individuals mark more near some groups than near others, for example, do individuals mark less near groups that are closely related? is it possible to show that points are clustering near some features but not near others?

Any help that can be provided would be greatly appreciated,

Dave Humphries

##### update

For demonstration, below is an image of four territories (A,B,C, and D) and the marking locations of two individuals from group A. Groups A and B are closely related. Groups C and D are unrelated to group A.


texts: Remove details from loaded OSM-layer? - Openlayers-3

texts:
I have created a tile server with help from this guide.
Sadly I still have no software to actually create the tiles and because of that I load in a online map (HOT) from OSM via Openlayers:

var map = new ol.Map({
      target: 'map',
      controls: [],
      layers: [
        new ol.layer.Tile({
          source: new ol.source.OSM({layer:'hot'})
        }), ... etc

Is it possible to hide details from this source?

I will eventually start using the data from my tile-server but I'm on a short deadline, so I'd like to evaluate my options...


texts: Open street map overlay white background

texts:
I have been working through trying to get the NOAA ridge radar image to act as an overlay. Looking at other entries on the site  I have worked with the image and have it projecting correctly but the png has a semi opaque white background I can't seems to get rid of.

I have tried setting transparent = true

it appears that the warp is adding the white background

gdalwarp -t_srs '+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs' radartemp1.png radartemp2.png

i tried also adding -dstalpha but then the layer comes out black...

anyone know how to solve this problem


texts: QGIS 3.10: forms_2 parents table/1 child table

texts:
I have 2 parents table (shape): one for trees (pk = id_tree); another for parks (pk = id_park).
I have 1 child table to list our works (fk = fk_works related to 'id_tree' and 'id_park').
In both cases, it's a one-to-many relations.
I declared both the relations in the properties of my QGIS project.
I build my forms.
But I can't add new work because, in the child table, the values of fk_works disappear.
I don't understand why and what to do to add new work for a special tree or park?


texts: How to adjust coordinates of duplicated layers on QGIS

texts:
I am having an issue with QGIS.

There are two layers (opaque green and brighter green) which are in different positions but are showing the same coordinates (tables below).

How do I update layer coordinates?

OBS.: those layers were drawn based on an orthomosaic and one layer is a moved duplication of the other.


texts: Opening KML file from ArcGIS server with GDAL?

texts:
I am able to open the google KML sample by running:

ds = gdal.OpenEx( fn, 0 )

But when using a specific file that I get from my country's database, it doesn't work. I've solved the problem by converting the file into a GPX using GPSBabel but I am sure there must be a better way of doing it. By the way when I run on my .gpx:

from osgeo import gdal, ogr

# File converted from KML to GPX by GPSBabel.
fn = r"D:\Downloads\POT_DONDE_VIVE_LA_GENTE.gpx"

ds = gdal.OpenEx( fn, 0 )
cnt = ds.GetLayerCount()

for i in range( 0, cnt ):

    print( ds.GetLayer(i).GetExtent() )

I get:

(0.0, 0.0, 0.0, 0.0)
(0.0, 0.0, 0.0, 0.0)
(0.0, 0.0, 0.0, 0.0)
(0.0, 0.0, 0.0, 0.0)
(0.0, 0.0, 0.0, 0.0)

That tells me that there must something wrong. The file works perfectly when using the ArcGISWebViewer.

So my questions are, how can I read KML files that output:

IOError: Could not open 'fileName'

with GDAL/OSGEO?, and if the only way is to convert them first to a different file format(how to it within python?) how can I be sure that the output of the conversion matches the original file?


texts: Can a layer's style be modified dynamically using an SLD in OpenLayers?

texts:
I am a complete novice with GeoServer and OpenLayers. I'm trying to ascertain at what point a style is applied to a GeoServer layer which is drawn in OpenLayers, when using an SLD.

The Introduction to SLD topic says:

  When adding a layer and a style to GeoServer at the same time, the
  style should be added first, so that the new layer can be associated
  with the style immediately

which implies that the style is stored with the layer in GeoServer. Does this mean that I can't change the layer's style dynamically in OpenLayers?

The Attribute-Based Polygon example shows how to symbolize a polygon based on its attributes. Could a user change this symbology dynamically, within an OpenLayers map which is running on their browser?

The OpenLayers Styling page says:

  When a feature is added to a map, its style information can come from
  one of three sources:

Does this mean that it would be necessary to remove a layer, and add it to the map with a new SLD, in order to change its style?

Thanks


texts: Using QGIS Vector->Geoprocessing Tools->Difference?

texts:
I'm trying to use the Vector Geoprocessing Difference Tool in QGIS to "take a bite" out of one polygon using another in another layer. The tool menu Difference icon is exactly what I'm trying to do - it shows a circle that it "taking a bite" out of the overlapping corner of a rectangle beneath it.

I've tried every combination of input in the dialog (swaping input layer and difference layer, using selected polygons, etc.) I get three results: The union of the two layers (!), all of layer one, or all of layer two. I only get results by selecting "use only selected features". When I don't "use only selected features" it simply churns away and never produces anything at all.

The dialog seems simple and straight-forward, the operation is pretty basic, so I wonder what I am doing wrong.

Is there some twist or trick to this that I'm missing?

I've looked for detailed tutorials and/or instructions, but have not found any.


texts: Problems using GRASS tools for topology validation in QGIS

texts:
I'm working with a shapefile of Enumeration Districts (Setores Censitários) of the 2010 Brazilian Census. Enumeration Districts (ED) are smallest areas for which Brazilian Census data is available. The 2010 Census had 318 thousand EDs. Created to facilitate Census logistics, EDs defined as areas containing from 250 to 350 households in urban areas and 150 to 250 in rural areas, so that they can be surveyed by one of Census surveyor. To avoid double counting or ignoring households, EDs must divide the space in a topologically consistent manner, with no overlaps or gaps.

To fix eventual overlaps and gaps, i'm using GRASS tools in QGIS, but the results are changing significantly the original topology.
Here is an example of topological errors in the shape: .

To import the shape vector to GRASS mapset i'm using the v.in.ogr.qgis tool  with 0.0001 snap threshold and 0.1 mininum size of area.
Just a few cases of overlaps and gaps are fixed and increasing those values didn't solve the problem, so i tried to use de v.clean tools to make these corrections in a controlled manner.

The v.clean.bpol didn't make much difference.
The v.clean.snap with 0.0008 threshold solve the gaps problems, but it changed the original topology in other areas, as you can see in the image:

Lower thresholds keeps the topology, but don't solve the gap problems.

The v.clean.rmarea with 10000 threshold solve some of the cases of overlap, but it erases some of the original areas that area very small as well. Lower threshold won't dissolve the overlap areas and greater threshold will dissolving original small areas.

There is anyway to fix these topology problems using GRASS or another tool but without changing the original topology in the process? A more controlled manner to do this validation?

Here is available a small portion of the shapefile that i've been working:
https://www.dropbox.com/sh/daodgcww4afr2w1/AABCUiHMnscNJPE3zRmjKfVJa?dl=0


texts: Adding labels with symbols to legend in ArcMap?

texts:
I'm creating a map with ArcGIS Desktop v10.2.2, showing museums as points. The museums layer has no symbol and the labels are made of a circle filled with color with the ID inside, taken from the [ID] field.

I want to add the entire list of museums in the layout map legend, having a list of my symbols with their IDs and the description (their name) taken from the [LABEL] field, without creating it manually.

I know that a possible manual solution would be to set the layer symbology to "Unique values", then add all values (more than 100!) and set once for each point the same label symbology. But I hope there's a better solution!

Take a look at my map below.

The layer is "Musei" (Museums) and as you can see on the legend on the left, it has no symbol. On the right you can see the descripted legend symbology, a Dark gray circle with Red text inside taken from the [ID] field. I want to add to my layout a column containing all the circles with IDs and the corresponding label.

-- UPDATE --

This is the Field Calculator settings I used following the solution proposed by @FelixIP.

But the points are not positioned in vertical align as I expect. Instead they're moved a bit here and there. I also tried using different step values: 10, 100, 1000, 10000.


texts: How to match Google, OpenLayers and shapefile data?

texts:
I've installed OpenLayers plugin and opened Google satellite layer in QGIS.

For a point in Beijing china, it's shown as 12961957,4852605 coordinate in blue status bar, and it should be 116.33710,39.83818 in longitude and latitude.

I also have some shapefile which are in lat & long coordinate for Beijing area. So how can I overlay the shapefile data on top of Google layer in correct position?

Anyone can help me or provide some hint?


texts: Moving symbols in Print layout, not in Main Window

texts:
I am working on production of Fire Evacuation maps. There are three assets, Fire Indicator Panel, Fire Mimic Panel and Emergency Exit that are located one very close to the other on a wall. They are all styled with SVG.

When I create a print layout the symbols overlap with each other (see picture below):

Is it possible to move the symbols in Print Layout so they look neatly one next to the other in a row (or one below the other in a column), without affecting the actual location of these assets?

I am running QGIS Noosa on Windows 10 Enterprise.


texts: Edit SQL Server table (not registered with geodatabase) with arcpy

texts:
I am trying to update some rows in a SQL Server database table (not registered with geodatabase) with arcpy.da.UpdateCursor. However, when running the deleteRow(), I get the error return without exception set error. Updating geodatabase tables works fine.

I guess it might be a limitation of arcpy which can't edit a non-registered with geodatabase RDBMS table (in this case, SQL Server), but I want to make sure I don't miss anything.

Can ArcSDESQLExecute class help me to update SQL Server table that are not registered in geodatabase (have found no info on that in the Help)? Any other workarounds?

EDIT:

I have ended up using ArcSDESQLExecute which works fine for executing SQL commands (including deleting rows) for tables that are not registered with geodatabase.

sde_conn = arcpy.ArcSDESQLExecute(r"C:\GIS\Temp\nongeo.sde")
sql_statement = """delete from added_table WHERE AddID = 2 """


texts:
I have a PostGIS database where I would like to store multipolygon data of certain cities.

The type of the table column is geometry(MultiPolygon,4326).

I'm now looking for a source where I can get the multipolygon data for cities like London or Berlin or a tool that supports me in creating this data.


texts: Why use st_union after st_combine with sf?

texts:
The documentation of the st_union function of sf is not clear (for me!) on one point: what are the action and result of the function when applied to an sfc with a single feature, with the parameter by_feature=FALSE (default)? It seems that, according to the text of the documentation, nothing should happen and that the value returned should be the same as the one-feature sfc provided.
I ask this question because it is rather common to find proposals, advices for chaining st_union to st_combine, especially in the context of the use of st_difference (see the sf manual, at geos_binary_ops, or ?st_intersection):
the use of "st_difference(SF1,st_union(st_combine(SF2)))" is proposed.
Two questions in order to understand:
. Because st_combine returns a single-feature sfc, how could st_union have any effect on such an object?
. Is there any simple case, any example, where st_union, with by_feature=FALSE, applied to a single sfc has any effect?
Extract of the documentation of st_union: "If st_union is called with a single argument, x, and by_feature is FALSE all geometries are unioned together and an sfg or single-geometry sfc object is returned. If by_feature is TRUE each feature geometry is unioned. This can for instance be used to resolve internal boundaries after polygons were combined using st_combine."


texts: Adding Rasters as Columns in Shapefile using ArcGIS for Desktop?

texts:
I am currently working on a research project in a lab, and I am fairly familiar with ArcGIS (I'm only in high school and have done some previous work using ArcGIS). I was asked to add 3 columns onto a data set (a csv file with XY coordinates). These three columns are supposed to come from rasters, and so far I have had no luck. So far I attempted to:

Extract data by points from the rasters
Join the files
Convert the csv file to a shape file and then attempted to join the three rasters (before I projected the rasters and after)
Adding XY coordinates to the shape file and then joining
Projecting the rasters and then join based on the FID/Object ID

I feel like this is impossible.


texts:
I've clipped a polyline shapefile with polygons, which produced a new feature class with obviously less polylines. However, the attribute table for the new feature class calculates the length of the polyline from the previous attribute table, which is clearly wrong.
For example, the new attribute table gives me the length of the previous polyline, even though only some of is highlighted in the new feature class. How can I find the length of only the highlighted polylines?


texts: number of pixels inside a polygon in Google Earth Engine

texts:
I'm trying to know how many pixels are inside my polygon when I calculate the mean (and other statistics) by using the "count" reducer. I guess many things go into this: what GEE does with pixels that are split, the right choice of scale, etc. I read somewhere that a pixel will be considered in a calculation when more than 50% of it is included in the polygon. I'm trying to understand what is going on with the following code:

https://code.earthengine.google.com/6144ac3715d195827378795c98994139

In have two polygons, polygon 1 and 2 that are inside the corresponding pixel. Polygon 1 has an area of 316 m2 and Polygon 2 an area of 360 m2, yet when I use the count reducer I get a value of 0 for the first one and a value of 1 for the second one. What is the rule here? Moreover, I still get a value for the mean even though the count is zero. What am I doing wrong?

I was trying to read the scale documentation about this, but I cannot understand what is going on.


texts: Why is an original ArcMap-project twice as big as its copies?

texts: I'm currently working with ArcMap 10.7. Therein, I regularly do File > Save a Copy of an existing ArcMap-project to create a copy for another ArcMap version. However, the file size of the original ArcMap-project (e.g. 3692 KB) is roughly double of the size of its copy (i.e. 1842 KB). The inconsistency of file size occurs in any direction, no matter if I copy ArcMap 10.3 > ArcMap 10.7,  ArcMap 10.5 > ArcMap 10.3 etc. What's the reason that the original file is of double size of its copies?


texts: Visual edit geographic data with geotools

texts:
With JMapFrame I have options for to visualize my geographic datas, but I don´t know how to add tools for to edit this data, equals conventional GIS. It´s possible to do this?

att,

Augusto


texts: Why is the vector style in the PDF export different than in print composer?

texts:
I'm exporting a map to PDF in the print composer.
See the map in the print composer

As you can see, the top map is different from the bottom; the top has 10 tracts shaded in, the bottom has 6.

Now see how it is when I export it to  PDF

Both the top and the bottom have the same tracts filled in!

Is it an export problem? Is it a layer problem? Is it another problem?


texts: If you create a series of circles/buffers in QGIS, what type of geometry do they go into PostGIS as?


texts: Tool, API of or code samples for calculating spatial similarity between Lines

texts:
I am migrating a database with a road network to an other external data provider.
We have joined our own data about roadworks, temporary speed-limits, etc., using ID and linear referencing on this data. So I need to find the same segment in the newer dataset and its ID to migrate the data.

There are small spatial differences and no matching attributes between these to datasets, so I can't use a regular spatial join. Parallel streets close to each other make using buffers a problem.

I looking for a way to find the similarity between lines, similar to fuzzy string searching with text. This way I can find the closest matching line, add the new ID, and recalculate the lengths in our tables.

Is there a tool, API, or codesamples (I prefer Python, JavaScript or PostGIS-SQL) that shows how to do this?


texts:
Does ArcGIS 10 have a Geoprocessor Programming Model by using Arcpy?
is it possible to use ArcGIS 9.3 Geoprocessor Programming Model within ArcGIS 10?


texts: Carto: How to set pop-up (hover) to display differently according to layers and zoom level

texts:
I'm working in Carto.com but I'm still new with this, therefore my question might be somewhat banal.
I have tried to generate a map of US with two layers: one layer reffering to the state and the other one for ZIP Codes. The layers change depending on the zoom level - if the zoom level is higher than 5 then you see the zip codes; otherwise you see the states.

I want the pop-up (hover) to be active for both layers, but only when each of them activate. But now the pop-up for the ZIP Code level is also displayed in the state level and vice versa.
This means that when I look at the US map by state and hover over it,I get displayed data about the ZIP Codes.

This is the CSS code I've added for the layer change (I worked with the visibility of the layer):

#idToHide[zoom<6] {

  polygon-opacity: 0;
  raster-opacity: 0;
  line-opacity: 0;

}

What should I do to set the pop-ups visibility accordingly to the layer?


texts: Check point in line segment not accurate in Shapely

texts:
I use Shapely and I've run into an issue where I can find the intersection of two lines. However once I have the intersection point, I don't seem to be able and cross-reference that it belongs in both lines.

Example

from shapely.geometry import Point, LineString

# Test two lines
l1 = LineString([(0, 300), (999, 300)])
l2 = LineString([(379, 250), (450, 320)])
point = l1.intersection(l2)

assert point.intersects(l1)  # True
assert point.intersects(l2)  # False

I wonder if this is a rounding issue? And in such a case how can I get around it? I.e. if I have a point, how can I guarantee that the point is in the line with an error of say 2 decimals?


texts: Map Refresh in Infrastructure Map Server 2015

texts:
I have a MapGuide application build in MapGuide 2011 where I refreshed my map using the following script

function MapRefresh() {
       debugger;
       ExecuteMapAction(20);
        parent.parent.mapFrame.Refresh();

    }

Now I am upgrading my application to Autodesk Infrastructure Map Server 2015 where the above script does not work,
Please guide me as to how I can refresh my main map from a panel?


texts: Apply Canny edge filter on raster failes with scikit-image

texts:
I have raster with one band and I would like to do edge detection with Canny edge filter.
This is my raster:

I have opened it with rasterio and covert it into array :
img=srs.open("raster.tif")
array=img.read()

Then i'm trying to apply the Canny filter as following (based on this: https://scikit-image.org/docs/dev/auto_examples/edges/plot_canny.html):
# Compute the Canny filter for two values of sigma
edges1 = feature.canny(array)
edges2 = feature.canny(array,sigma=3)

but it fails with this error:

ValueError: The parameter image must be a 2-dimensional array

I  read that this can happpen if my image is not in grey scale, but I have shown it as grey scale:
show((img,1),cmap='Greys_r')

so I tried to checl what is my array type:
type(array)
>>>numpy.ndarray

so I have no idea why I get this error message saying the image is not 2 dimensoinal as it is 1 band with X  and Y axis.
My end goal - to be able to apply Canny edge on my image and to understand why I get this error


texts: Updating ol.layer.Tile parameters in OpenLayers 3.9.0

texts:
I have a TileWMS:

var mapWMSTile = new ol.source.TileWMS(
            {
                url: '.....',
                params: {'LAYERS': 'l1,l2', VERSION: '1.3.0', 'TILED': true},
                serverType: 'mapserver'
            });

        var fieldTile = new ol.layer.Tile(
            {
                source: mapWMSTile,
                name: 'fieldTile'
            });

Now I want to add a layers to my WMS parameters (e.g. add l3 in 'LAYERS'). Is this possible?

I have tried the below codes without success:

First try:

mapWMSTile.updateParams(
                {
                    params: {'LAYERS': 'l3', VERSION: '1.3.0', 'TILED': true}
                });

Second try:

fieldTile.getSource().updateParams({'LAYERS': 'l3', VERSION: '1.3.0', 'TILED': true});

I am using OpenLayers 3.9.0.