texts: Intramaps Roam GPS_LOG not honoring Tracking settings

texts:
I have Roam 2.5.2 and enabled the gps_log layer as discussed here: https://roam-docs.readthedocs.io/en/latest/features/Features-2.1/#gps-logging-and-tracking

For most part it works well, however it seems to be recording points every second, even though the setting 'GPS Feature Tracking' is set to Every 10 second.

After an hour or two we have almost a Gig of point data, which is unusable in that it takes forever to load while panning the map.

Is there another way to set this value?


texts: Blend differing line thicknesses in QGIS?

texts: I'm mapping traffic flow along a series of roads (links), styling the line thickness to reflect traffic flow. Along the same road (adjoining links), traffic flow can differ, as other roads feed traffic in or out. I've currently got something that looks like the top part of the image below, but I want something that looks like the middle or bottom bit (where the two line thicknesses merge). Is this possible in QGIS (2.18.6), and if so, how?


texts: Creating an intersecting linestring layer from two polygon layers using some tolerance

texts:
I have two postgis layers:

The first layer contains the sidewalks of my town (polygons in blue).
The second one contains the pavements (roads) of my town (polygons in brown).

The map looks like this in QGIS with both layers activated:

Ok, now I have to obtain a layer with the curbs (lines) in the town starting from these two layers. The new layer must be composed of linestrings. Each element in this new layer should be a linestring where every of the greometries of the two layers intersect.

In the following snapshot I manually painted the geometries that should be created in this particual area:

Sorry about my bad pulse. I hope it helps to understand me better.

Ok, if I just add the layers to QGIS, for example, and intersect one layer with another I get the following layer:

The new green layer is the obtained curbs layer. I have two snags with this layer:

It is composed of polygons
Not all the curvs are drawn.

With this last evidence I discover that the sidewalks and roads layers were not precisely delineated one next to the other, so I need some kind of tolerance in the operation to get the layer I am looking for. As I told before it should be a linestring layer.

I don't mind if the query is directly done in postgis or using QGIS or other software's tools, I just can't figure out right now how to solve this.


texts:
Edited:

I want to ilustrate my question. Suppose that you're in "Point A" and want to go to "Point B". This points wouldn't be in "at_2po_4pgr" table cause aren't source/targets nodes.
Then, I would search for the nearest node for points A and B (green points). After this, I might perform a shortest_path call using green points ids and I would obtain "orange" path.
But to obtain the real path cost (distance) in first case I would have to substract "offsetA" and in second case add de "offset B". To calculate distance between red-points and green-points, I run the following query:

SELECT * FROM st_distance(
ST_GeomFromText('POINT(-3.6963314 42.3498066)',4326),
ST_GeomFromText('POINT(-3.6954276 42.3479634)',4326)).

How I would know when add or substract the offset?

Sorry for my english!


texts:
I am creating a model where the first step is creating a geodatabase with location and name set as parameters. After that I created another step (Table to Table) to import a single table to the created GDB. After that, this table is going through many steps and tools where I have many outputs afterwards.

How I can assign the location of the many outputs produced to be in the created GDB at the beginning?

I tried the "%" sign but it didn't work. If the "%" is the solution, please give me an example on how to use it correctly.

I am using ArcGIS 10.2.2 for Desktop.


texts: What is the meaning of distance between a point and polygon while using ST_Distance? I'm aware that it means the shortest distance between them but what is here meant by shortest distance?


texts: What is syntax of bool in arcpy.sa.ExtractByAttributes?

texts:
I want to extract specific values of raster layer. Use codes as follows:

attExtract=arcpy.sa.ExtractByAttributes(outraster,'VALUE==8 | \
                                            VALUE==9 | VALUE==10 | VALUE==16')

system hint the SQL is invalid:

ExecuteError: ERROR 999999: Error executing function.
An invalid SQL statement was used.
An invalid SQL statement was used.
Failed to execute (ExtractByAttributes).

How to write the operation bool code?


texts: Trimble GPS data is not displaying in ArcMap attribute table

texts:
I am currently having difficulty trouble shooting a problem. While using a TRIMBLE GEOXH (GeoExplorer 2008 Series) I have the problem of missing data after exporting it to ArcMap from my GPS receiver. It seems my GPS receiver is NOT collecting; Elevation, Vert Prec, Horz Prec, Latitude and Longitude for some GPS points.

How does one get this data to return??


texts: Preserving Vector Attribute Types When Processing Vectors

texts:
Qgis seems to change vector attribute types for certain processes. Is there a way to force preservation of the attribute types?

My situation:
I have five vector layers I want to merge, all shapefiles, four of which have identical geometry and attribute types, and one whose only difference is that the geometry is Multipoint rather than Point. (The attributes are all Strings, except one which is Integer64.) I used "Vector/Geometry Tools/Multipart to singlepart" to convert the one layer, then tried using "Vector/Data Management Tools/Merge vector layers" on the four original layers and one new one.

What happened:
The merge failed because the attribute data types no longer matched: the new layer had a String type for what should have been the Integer64 attribute. I fixed my attributes, ran the merge again, and this time it worked fine, except that the merge also changed the Integer64 attribute to String.

The layers I made with these processes were all temporary layers. I then tried saving to a Spatialite table, to a PostGIS table, and to a permanent (named) shapefile; the only one which preserved the attribute type was the save to a file.

Is there something wrong here, or am I missing something when I'm doing this?  Other than saving to a file, how do I preserve my attribute data types?

(BTW I am using QGIS 2.18.8 on Linux Mint 17.3.)

EDIT
Sample shapefiles can be found zipped here. "Test1" is the multipoint shape, with processing result in "multi2single_t1," while "test2" and "test3" are merged in "merged_t2_t3."


texts: How to make a polygon and point shapefile line up to do a spatial join?

texts:
I have a polygon with a North_America_Lambert_Conformal_Conic coordinate system and a point shapefile which had a WGS 1984(DD) coordinate system.

I defined projection and projected both of them to North_America_Lambert_Conformal_Conic coordinate system. The two shapefiles do not overlay. If checked in the screen shot of ArcGIS 10.2, only one of the layer appears (whichever I have zoomed to).

Can someone please guide me on how I can get this to shapes to line up or overlay?

My intention is to do a spatial join on them so that I can count some point features within the polygons.


texts: Creating point file along line (trail) at regular distances (0.20 miles) using ArcGIS Desktop?

texts:
I'm using ArcGIS 10.1. I have a trail network that I'm trying to establish signs along the trails every 0.20 miles.

Is there a way to create the point file from the line feature?


texts: How to convert int to string in CQL filter GeoServer?

texts:
I want to filter columns in GeoServer WMS layer using 'ilike'.  I tried following cql query.

"column_name" ilike 'string'

It works for columns of datatype string, not working for integer columns.

How to convert integer column into string using CQL ?

UPDATE

Usecase Example:
 Lets say I have two columns in a WMS layer such as one is address and another is zipcode. Both address and zip code has numbers in their column, but zipcode is of datatype integer and address is of datatype string. Below is sample column.

Column

address

338/4 -yyy,zzz.

444/3 - zxy

323 - abc.

zipcode

33431

385493

30293

Now I want to run CQL query to pick rows even if one column has a search string. My query for search string '38' is like follows.

Query

address ilike '38' or zipcode ilike '38'

For the above query I want the result which contains first and second rows.

Result

address

338/4 -yyy,zzz.

444/3 - zxy

zipcode

33431

385493

Like this many layers, and so many columns. So changing query for each column is not easy.


texts: Expanding polygons and clipping at midpoint?

texts:
I have several separate shapefiles with polygon features. The end goal is to maintain separate shapefiles for each polygon but form a contiguous, non-overlapping coverage for my area. Basically, I guess I am trying to fill in gaps using a consistent, repeatable process.

The polygons are irregular:

Initially I thought I could buffer each polygon and do something like this for the overlaps: Cut off overlapping polygon areas at their cutting edge but was unable to perform using v.clean. Is there another way?  Most methods I see seem to favor one polygon over the other and I'm trying for a standard or middle distance.

If a buffer is not used, I want to expand each polygon in a standardized way and have them meet at a 'midpoint' between each polygon. Is this possible? Maybe something with node snapping tolerances? I am currently using QGIS 2.16 and don't have any experience using python.


texts: How do I identify and gracefully kill hung Gp service processes?

texts:
I've published a gp tool as an asynchronous non-pooled gp service.

When the gp tool hangs I would like to allow the user to kill it, without impacting other services.  If I put a webservice on the SOM that kills a process based on an ID, will it harm other services?

Related questions:

Is there a way for a webservice running on the SOM machine to kill a process on a (different) SOC machine?

The first thing the gp tool does is write out the hostID and processID to a log file, but I'd like a system to work for any gp service, not just ones that log this info.

Is there a request I can make to the gp server that returns a hostID & processID for an arbitrary jobID?

Update:  I've decided to first implement a method to Cancel GP jobs, and then only pursue the process killing strategy if canceling proves ineffective.  Since the REST API doesn't currently have a Cancel method, I will be doing this via SOAP.  I've also been shown how to log the process IDs, so that way I won't need to rely on my own logging.


texts:
I am trying to interpolate distribution value of rain using IDW (inverse distance weighting) method using quantum gis. What steps do I need to undertake?

The input layer data which I used there are 2 as same as the shapefile data.

"02StaHujan" as point input
"Balai" as the boundaries which I'll interpolate.

I have some problems when I'm trying.

Why I can't set the Z field? the Z field isn't active
How values must be filled in number, row, cell X and Y columns? I want the results of interpolation cell size is 30m
why the results of my experiment like that? Looks pict. "Capture4"


texts:
I want to export a GeoTIFF with PyGRASS but i can't figure out how I can do that.
my code so far:

# import grass python libraries
from grass.pygrass.modules.shortcuts import general as g
from grass.pygrass.modules.shortcuts import raster as r
from grass.pygrass.modules.shortcuts import vector as v
from grass_session import Session

inputfile = './data/gtifs/13.tif'

with Session(gisdb="./tmp", location="location", create_opts="EPSG:4326"):
    r.external(input=inputfile, output="layer_1", overwrite = True)
    print(r.info(map="layer_1").get_dict())

    r.mapcalc(expression='if(layer_1>0, 1, null())',overwrite=True)

    print(g.list(type='raster', flags="m").get_dict())
    # neither of the approaches below work
    r.external.out(directory='./', format="GTiff")
    r.out.gdal(input='layer_1', output='./test.tif')

r.out.gdal() produces an FileNotFoundError: [Errno 2] No such file or directory: 'r.out': 'r.out' and an GrassError: Error running: 'r.out --interface-description'. error

r.external.out() does not seem to exist as a method

how do I save my grass raster as an GeoTIFF? I'm using the Grass GIS 7.8 binaries and python 3.6.9


texts: How to move line segment to certain direction using ArcGIS for Desktop?

texts:
How can I move a Line Segment to a certain direction?
In my case i want to move the Line up north, so i need to change the y-coordinate somehow, I'm just not sure how to do this in ArcMap. Im using ArcGIS 10.1


texts: Setting up Concatenation using Python expression in field calculator?

texts:
I am having trouble creating a Python expression in Field Calculator in ArcMap 10.3.1. I’m attempting to concatenate two field values (double and short integer) to a third (double) field, (CON_IN). I need to allow only two digits to be added to the system number, so if the ET_NO value is greater than or equal to 10, do not add “0”. If the ET_NO is less than or equal to 9, add “0” and the value.

The final output should look like this:     2161400187509  or  2161400187528

I’ve tried variations of this script but get basically the same results.
The resulting error shows a syntax error and that SYS_NO is not defined.

Does anyone have any advice?

Pre-logic Script Code:

def   Concat(x,y):

  x= SYS_NO

  y= IN_ETNO

  if y >= 10:

    return x +  y

  elif y<= 9:

    return x +  "0" +  y

Expression:

CON_IN =  Concat( !SYS_NO! , !IN_ETNO! )

I see I didn't need to define the x,y. I see by the other posts that I am trying to concatenate a number and a string value, but when I remove the quotation marks on the string it just sums the SYS_NO and ET_NO values which is not what I need. I am not adding the values together, I am trying to add a two digit code to the end of a system number based on unique values for each record.


texts: Finding shortest route for cyclist that can switch between road layer and cycle path anywhere they intersect?

texts:
I have two data layers that I want to use in a network analysis for the shortest route analysis for a cyclist:

road network
cycle paths

I want to find a shortest route for a cyclist that can switch between road layer and cycle path anytime where they intersect.

Problem is that sometimes the cycle network lies on the top of the road network but does not align properly. I am not sure how to use the data in network analysis.

Should I merge these two data sets first and use the feature to line tool to create intersections of the street and cycle network, or is this a case of a multimodal network?

I work with Arc GIS 10.2


texts: Solving Error " destination provider" when trying to save WMS as GeoTIFF

texts:
I'm new and I use QGIS 3.10

I would like to save the WMS output as GeoTIFF to be able to load and clip the raster by a mask layer (polygon). When I try to save the WMS output as GeoTIFF (Export--> Save as) I get the message: "cannot write raster. Error: destination provider".

This is the window when I set the parameters:

In "name file" it makes me choose a folder in which to save the raster.
When I click "ok", I get the error: "cannot write raster. Error: destination provider"

Can you help me?

I tried to save in other folders but the error remains.


texts:
I have a layer with many points on a country, for example, Peru, and I have done a grid grouping all the country. Each hole of the grid has an ID so it's easy to differentiate them. My idea is to set that ID for each point from the other layer that way I can work with all the point on the same hole. In SQL it would be to make a new column called 'ID_Hole' and write the grid ID, but I'm not getting it on QGIS.

Any idea about what should I do?

I attach a more explicative image

The point is obtained from a.csv where I can see the name of the settlement(city) and its located by latitude and longitude on the map.


texts: Only Integers Transferred when Joining by Location

texts:
When joining attributes from a points layer to a polygon layer of Postcodes, only the integer fields from the points layer are transferred. See screenshots below.

Points Layer

Points Layer Continued

Postcode Layer

Setting up the algorithm...


texts: Methods for hot spot or cluster analysis for multiple measurements at the same geographic unit?

texts:
I would like to find geographic clusters of zip codes where each zip code has a measurement of a person's age (integer), sex (categorical/binary), and if they recovered from a certain medical condition (categorical/binary).  Each zip code could have multiple people in it.

Will having multiple records associated with each geographic unit pose a problem for a cluster or hotspot analysis?  Are there methods for dealing with this situation?  Aggregating the data for each zip code (average or counts) doesn't seem appropriate.

I'm not tied to any specific software.


texts: polylines between two locations over sea and not land

texts:
I have several thousand locations in a database, Location A travelling to Location B, New York to London, Paris to Buenos Aires.

I'm trying to plot a path between two places only over sea. I'm not looking for strict shipping lanes or that the start and end points are in a port, if the start/end point is close to New York then that's good enough.

I've been looking for an API service to do this work, Google maps provides a path between two locations but you can't restrict the path to ocean only.

Are there any services where you can plot a very simple polyline between two locations over sea programatically?


texts: Extracting MultiValuesToPoints yields identical value for each band?

texts:
I'm attempting to extract hyperspectral remote sensing data to a point featureclass using the ExtractMultiValuesToPoints ArcGIS tool.

The value of each band (228 bands) should be appended as a field to each point in the input feature class. The problem is: each appended field contain identical values for each point.

So it seems that only the value of the first band is extracted and assigned to each next band.

I've scripted the tool as follows.

import arcpy
import os
inRasts= "L:/MyData/StudySite.IMG"
fcDir= "L:/MyData/Geo.gdb"
fcIn= "Plots"
inPoint= os.path.join(fcDir, fcIn)
arcpy.sa.ExtractMultiValuesToPoints(in_point_features= inPoint, in_rasters= inRasts, bilinear_interpolate_values= "NONE")

The Remote sensing image properties are:

data type: File System Raster
Format: IMAGINE Image
pixel type: floating point
Admittedly, the raster data may be a bit weird. The NoData pixels have a pixel value of -3.40282346639e+038 I don't quite know how I ended up with that, must be a rounding error somewhere.

I'm using ArcGIS 10.1 on W7 computer. I cannot find any service packs installed.


texts: Converting GeoPDF to ArcGIS with annotations (arcmap original labels)?

texts:
I've got a GeoPDF file that contains vector and annotation data (arcmap labels). I need to convert it into a feature class.

I've tried to follow the former questions that been asked here:

Adding PDF (especially GeoPDF) to ArcMap as data?
Converting GeoPDF to vector format like CAD or shapefile?

I've succeeded to convert the file into .dwg/.dxf files by using the "Any PDF to DWG" software. The output was an ungeoreferenced and annotations-less .dwg.

How to convert the GeoPDF file to .dwg/.dxf without loosing the annotations layer?
How to keep GeoPDF's source georeferencing?


texts: What is the the color gradient at borders called?

texts:
Many paper maps have an effect around the borders where there is a gradient from the darker color of the borderline to the lighter color of the area fill. What is this effect called?

See https://www.wwu.edu/huxley/spatial/tut/cart/adobe_illustrator.htm#inset for an example.


texts: Converting mapsforge .map file to OSM or shp

texts:
I have found solution for converting OSM file into .map file but is there a way for converting .map file back to osm or shp or any other well-known vector file type?

by .map file i mean mapsforge .map file, not mapinfo .map file


texts: How large is the entire ASTER Global DEM on disk?

texts:
I see from various online sources (including GIS.SE's own How to download the entire ASTER GDEMv2 dataset?) that the ASTER GDEM has 22'000 tiles, each 1°×1° and 3601×3601 pixels (each pixel being 16-bit integer), with both elevation and quality rasters.

A handy napkin says the uncompressed data would be 22000 * 3601 * 3601 * 2 * 2 / 1024**3, or ~1063 GB. But how big is the compressed data that you download from the official providers?


texts: draw lines with an angle inside a rectangle in python

texts:
I am trying to make a plugin for Qgis with python. I have the four corner coordinates for a rectangle. I want to give the user the inputs, line space and angle (0 -180 degrees).

With these inputs, I want to fill this rectangle with lines in the direction of the specified angle and with the specified line spacing between the lines.

I guess what I need is the same type of code as for making a grid. But I can only find information about grids which are vertical or horizontal. I will try to code it myself, but if anyone have a webpage with information, I would gladly to read it.


texts: Debug In-dialog provided python code in QGIS attribute form 'Python Init Code Configuration'

texts:
In a QGIS layers attribute form configuration I have the following Python Init Code Configuration:

Code in plain text:

# -*- coding: utf-8 -*-

from qgis.PyQt.QtWidgets import QWidget, QLineEdit
from qgis.PyQt.QtCore import Qt

def atlas_punkte_init_values(dialog, layer, feature):
    #geom = feature.geometry()
    control = dialog.findChild(QLineEdit, 'str')
    print(control)
    l = dialog.findChildren(QLabel)
    print(l)

As one can see at first sight, this code is erroneous, because QLabel is not imported. By means of this, I would expect an exception somewhat like ...NameError: name 'QLabel' is not defined, but the only thing I get is this:

Because far trickier error constellations are imaginable than a forgotten import, debugging of more extensive scripts is nearly impossible that way.

First idea is to look at python errors in the log or the python console, but thats a fail.

Do I miss a location where qualified exceptions are thrown by these scripts or is there really no chance to properly debug these?


texts:
I have two questions related to calculate building density.

I need to calculate the building density in a map which Data is in shape format. The buildings and the area are in two different layers.
Same question as above but the map/area is raster format.

Thank you.


texts: Reprojecting coordinates using MakePoint() with SpatiaLite?

texts:
In my SpataLite-DB, I've got a table holding coordinates in EPSG:2056, splitted in two columns "x_coord" and "y_coord".

I want to copy them into another table, using MakePoint to convert the two coordinate columns into a geometry (EPSG:21781) column.

Is it possible to convert the coordinates during the copy process, or do I have to first convert and then copy  them?


texts: How to add vertices to a polygon on mouse click in OpenLayers

texts:
I am pretty new to OpenLayers.  I have a polygon drawn on a vector layer, along with ModifyFeature and SelectFeature, I need to implement a functionality when user clicks on the polygon, a new vertex should be added to the polygon at the position where use clicked.  With the new vertex the user should be able to  reshape the polygon. Is there any existing feature to do this? if not how can we implement this?

function init() {
    map = new OpenLayers.Map('map');

    var wms = new OpenLayers.Layer.WMS("OpenLayers WMS",
        "http://vmap0.tiles.osgeo.org/wms/vmap0?", { layers: 'basic' });

    OpenLayers.Feature.Vector.style['default']['strokeWidth'] = '2';

    var renderer = OpenLayers.Util.getParameters(window.location.href).renderer;
    renderer = (renderer) ? [renderer] : OpenLayers.Layer.Vector.prototype.renderers;

    vectors = new OpenLayers.Layer.Vector("Vector Layer", {
        renderers: renderer
    });

    map.events.register("mousemove", map, function (e) {
        var position = this.events.getMousePosition(e);
        OpenLayers.Util.getElement("coords").innerHTML = position;
    });

    map.addLayers([wms, vectors]);
    map.addControl(new OpenLayers.Control.LayerSwitcher());
    map.addControl(new OpenLayers.Control.MousePosition());

    if (console && console.log) {
        function report(event) {
            if (event.type === "beforefeaturemodified")
                controls.polygon.deactivate();
            console.log(event.type, event.feature ? event.feature.id : event.components);
        }
        vectors.events.on({
            "beforefeaturemodified": report,
            "featuremodified": report,
            "afterfeaturemodified": report,
            "vertexmodified": report,
            "sketchmodified": report,
            "sketchstarted": report,
            "sketchcomplete": report
        });
    }

    controls = {
        polygon: new OpenLayers.Control.DrawFeature(vectors,
                    OpenLayers.Handler.Polygon),
        modify: new OpenLayers.Control.ModifyFeature(vectors),
    };

    var featureSelect = new OpenLayers.Control.SelectFeature(vectors,
    {
        hover: false,
        clickout: false,
        multiple: false,
        toggle: true,
        //box: true,
        eventListeners: {
            featureunhighlighted: function () {
                //Do whatever you want here.
            }
        }
        //multipleKey: "shiftKey"
    });
    map.addControl(featureSelect);
    featureSelect.activate();

    for (var key in controls) {
        map.addControl(controls[key]);
        controls[key].activate();
    }

    map.setCenter(new OpenLayers.LonLat(0, 0), 3);

    //vectors.events.register('click', vectors, function (e) {
    //    console.log(e);
    //});

    vectors.events.register('sketchcomplete', vectors, function (e) {
        controls.polygon.deactivate();
    });
}

Thanks in advance


texts: MapBox Text Layer is adding multiple text instead of single centered text

texts:
I am trying to use data from my tileset (.zipped shape file) to show Text Label over my tile layer. However, instead of adding a centered label, it's adding multiple-texts sometimes.

For example 38 or 37.

If I zoom, it adds more

How can I force it to have only single label in the center of the tile instead?


texts: Multi-User offline edit synchronization using GPKG

texts:
I have been trying to find an article like http://postgis.net/workshops/postgis-intro/history_tracking.html that explains how I can do the following with a GPKG (ideally in QGIS)

Sync edits made on and offline by multiple users
See which user edited which records without having an auto update column in the table with this info or using http://geogig.org/

I understand that some of this is part of QField Sync but when I try it using QGIS 3.11, it just opens all the tables additional times and I don't see the modified data. I also want to understand how syncing and having an audit trail in a geopackage works as explained in this talk https://media.ccc.de/v/bucharest-262-field-data-collection-strategies-using-qfield-and-qgis#t=743.

I have reviewed http://www.geopackage.org/spec121/ and it mentions syncing but there is no info on how it's practically achieved.

How do I sync edits made offline by multiple users and have an audit trail?


texts: How to get raw values from OGC Web Coverage Service (WCS) via ArcGIS for Desktop 10.1?

texts:
I posted this a few days ago to the ArcGIS Discussion Forum but have had no response so thought I would try here as well.

I'm new to Web Coverage Services (WCS) but have found it easy to add one using ArcGIS for Desktop 10.1 and this service URL:

http://nsidc.org/cgi-bin/atlas_south?service=WCS&request=GetCapabilities&version=1.1.1

I add the service by Liu et al. and get a nice display.

However, when I use Identify to look at the values for any cell on the grid/coverage all Stretched Values and Pixel Values are in the range 0-255. I would expect to see 0-255 used for the Stretched Values but thought the Pixel Value would be holding raw data with values well in excess of 255.

When I export data from the WCS to a file geodatabase raster the values there are 0-255 too.

Is it possible to access the raw values held in a WCS?

(Hopefully I am just misunderstanding how to get my desired result)

NSIDC have said that their server is configured to provide data in GeoTIFF format. Their docs make it appear that GeoTIFFInt16 and GeoTIFFFloat32 are both options but they suspect that the GetCoverage request generated by ArcGIS for Desktop 10.1 is using a format that maps the data into an 8-bit format instead.

I do not know of an option within ArcGIS for Desktop to select the GeoTIFF format used by GetCoverage or a way to trap the GetCoverage URL that it is passing to the WCS server.


texts: Cannot create geodataframe from feature class due to IndexError: list index out of range

texts:
I have a function that helps me read in large feature classes and them create a geodataframe:

line_feat = "S_Studies_Ln"
cnms_cols = ['REACH_ID', 'STUDY_ID', 'CO_FIPS', 'CID', 'FLD_ZONE',
         'VALIDATION_STATUS', 'MILES', 'STATUS_TYP','STATUS_DATE',
         'STUDY_TYPE', 'TIER','LINE_TYPE', 'BS_ZONE', 'BS_STDYTYP',
         'PRELIM_DATE']

def records(filename, cols):
    """DESCRIPTION HERE"""
    with fiona.open(filename, layer=line_feat) as source:
        for feature in source:
            f = {k: feature[k] for k in ['id', 'geometry']}
            f['properties'] = {k: feature['properties'][k] for k in cols}
            yield f

cnms_df = gpd.GeoDataFrame.from_features(records(input_gdb, cnms_cols))

This works fine on all my feature except one. For a feature class called "FL", it throws the follwong error: IndexError: list index out of range. Any ideas as to what may be causing this in a feature class or how to debug it?
Here is the full output:

> IndexError                                Traceback (most recent call
> last) <ipython-input-12-8059b0bcbbb8> in <module>
> ----> 1 cnms_df = gpd.GeoDataFrame.from_features(records(input_gdb, cnms_cols))
>       2 #cnms_df = pd.DataFrame(cnms_df.drop(columns='geometry'))
>       3 cnms_df.tail()
>
> ~\AppData\Local\conda\conda\envs\geo_env\lib\site-packages\geopandas\geodataframe.py
> in from_features(cls, features, crs, columns)
>     231                 f = f
>     232
> --> 233             d = {'geometry': shape(f['geometry']) if f['geometry'] else None}
>     234             d.update(f['properties'])
>     235             rows.append(d)
>
> ~\AppData\Local\conda\conda\envs\geo_env\lib\site-packages\shapely\geometry\geo.py
> in shape(context)
>      37         return MultiPoint(ob["coordinates"])
>      38     elif geom_type == "multilinestring":
> ---> 39         return MultiLineString(ob["coordinates"])
>      40     elif geom_type == "multipolygon":
>      41         return MultiPolygon(ob["coordinates"], context_type='geojson')
>
> ~\AppData\Local\conda\conda\envs\geo_env\lib\site-packages\shapely\geometry\multilinestring.py
> in __init__(self, lines)
>      50             pass
>      51         else:
> ---> 52             self._geom, self._ndim = geos_multilinestring_from_py(lines)
>      53
>      54     def shape_factory(self, *args):
>
> ~\AppData\Local\conda\conda\envs\geo_env\lib\site-packages\shapely\geometry\multilinestring.py
> in geos_multilinestring_from_py(ob)
>     121     exemplar = obs[0]
>     122     try:
> --> 123         N = len(exemplar[0])
>     124     except TypeError:
>     125         N = exemplar._ndim
>
> IndexError: list index out of range

I have tried checking the attribute table for the feature class in ArcMap but can't see any difference between it and the other feature classes.


texts: US Census Bureau Shapefiles with district IDs

texts:
I am loading senate district polygon shapefiles from https://www.census.gov/cgi-bin/geo/shapefiles/index.php?year=2019&layergroup=State+Legislative+Districts into my Leaflet map. Well, to be accurate, I am first converting these shapefiles to GeoJSON with mapshaper.org. Anyway, we would like to color individual districts based on statistics such as income level. However, the shapefiles do not contain district IDs.

Are the polygons in these shapefiles organized in a consistent sequence that can be mapped to district IDs?

If not, is there another way to know to which district each polygon corresponds?


texts: How to correct sentinel data with sen2cor in R using a json with parameters?

texts:
Reading the sen2r documentation I saw the next example:

library(sen2r)
library(sf)

# Download an L1C SAFE product
example_s2_list <- s2_list(
  spatial_extent = st_sfc(st_point(c(12.0, 44.8)), crs=st_crs(4326)),
  tile = "32TQQ",
  time_interval = as.Date(c("2017-05-01","2017-07-30"))
)
s2_download(example_s2_list, outdir = tempdir())

# Correct it applying a topographic correction
sen2cor(
  names(example_s2_list)[1],
  l1c_dir = tempdir(),
  outdir = tempdir(),
  use_dem = TRUE
)

In the example the code first checks the products available, then it downloads them and finally with sen2cor it applies some corrections.
My question is if there is a way to put the parameters to correct that from a json previously exported from the GUI with the desired parameters


texts:
I have a parcel layer (blue polygons) and a road centerline layer (black lines).

I want to generate polygons that represent the portion of the road right-of-way that is adjacent to each parcel (orange polygons).

Question:

Is it possible to generate road-adjacency polygons via ArcGIS Desktop 10.7.1?

Idea:

Convert the parcels polygons to lines.
Eliminate the part of the parcel lines that are adjacent to the road right-of-way.
Use the remaining parcel lines to project-out to the road centerline.
Use those projected lines to cut out the parcel-adjacency polygons.

Would something like that work?

Note: I also have a road-right-of-way polygon layer that might be useful.


texts: How to travel x (e.g. 10) miles and get new co-ordinate having direction towards another location?

texts:
I'm trying to create a large but optimal circle that covers two intersecting circles. For that, I need to calculate new radius and new center point. I'm able to get the optimum radius by how can I get the new center point coordinate? I think we need to travel some X miles from one point to another. I have the calculation for X miles, but how to get the coordinate?
See below illustration. I have points A and B as geo coordinates. I wanted to travel x miles (let's say 10 miles) from Point A to Point B and get the co-ordinate of point P.
We can easily get the midpoint coordinates by (x1 + x2) / 2 formula which is 50% traveled from A to B but what if I want to travel 10% from point A towards point B?


texts: Splitting polygon in half with lines in QGIS 3.12

texts:
I have a polygon and a line, I want this line to split the polygon in half, i.e the same thing I could manually do with "Split Parts".

I have tried every Geoprocessing tool; most answers I found suggest using a clip or intersection tool, but to me, that only results in a new line, which is the opposite of what I'm trying to achieve.


texts: PostGIS, no use of gist for geometry column in a query

texts:
My query is:

  SELECT count(datasetid)
  FROM scenes
  WHERE acquisitiondate between '2014-06-01' and '2015-06-01'
  AND ST_Intersects(ST_GeographyFromText('SRID=4326; POLYGON((' || '10.66437' || ' ' || '49.981172' || ', ' || '15.46437' || ' ' || '49.981172' || ', ' || '15.46437' || ' ' || '54.781172' || ', ' || '10.66437' || ' ' || '54.781172' || ', ' || '10.66437' || ' ' || '49.981172' || '))'), bounds)
  GROUP BY datasetid;

datasetid has the data type smallint and acquisitiondate is a timestamp (without time zone).

There is a gist-index on the column bounds. The index won't be used if the column has geometry data type. But if the column has geography data type than the gist index is used instead of a sequential scan.

What is the reason for this?

Maybe this can be helpful:
PostGIS Version is 2.1.8 and PostgreSQL 9.4.5.


texts: arcpy.MakeFeatureLayer_management: query throws error 000230/000358

texts:
Im trying to create a feature layer on which i want to perform a select by location, however:

arcpy.MakeFeatureLayer_management (perp, "selected","LINK_ID = 23002109")

throws:

ExecuteError: ERROR 000230: Failed selecting with LINK_ID = 23002109

and

arcpy.MakeFeatureLayer_management (perp, "selected")
arcpy.SelectLayerByAttribute_management("selected", "NEW_SELECTION", "LINK_ID = 23002109")

throws:

ExecuteError: ERROR 000358: Invalid expression

The same query works fine when done manually in ArcMap. The field exists and its datatype is DOUBLE.

My problem seems to be exactly the same as in ArcPy inconsistently accepts SQL queries for MakeFeatureLayer or SelectLayerByAttribute The difference is I'm working with a local GDB and not with a SDE connection. What is wrong here?


texts: Inserting only one band in Google Earth Engine and modify that band

texts:
I am trying to see the heat inside Laos watershed and I would like to only use Band 10 in Landsat image. But the watershed is located within the intersection of four different Landsat images. So if I want to process "Band 10" for the four images, how to process that in GEE?

This is my snippet code.

 //Display laos admin from shapefile
var laosadmin;
Map.addLayer(laosadmin, {color: 'green'}, 'laosadmin');

//import image colllection
var lcCol = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR');

//Filter image by date and region of interest
var filterdate = lcCol.filterDate('2019-01-01','2019-03-01').filterBounds(laoswshed);

//Display landsat into the screen
Map.addLayer(filterdate);
//Display watershed area
Map.addLayer(laoswshed);
//Center the view to Laos Admin for macrosight
Map.centerObject(laosadmin);

//see the list of the date at the console
print(filterdate);

//Filter image collection by the IMAGE_QUALITY property
var filtered = lcCol
.filterMetadata('IMAGE_QUALITY','equals',2);
print(filtered)

// I just randomly pick one landsat image ID from the image collection, loading one image
//var image = ee.Image('LANDSAT/LC08/C01/T1_SR/LC08_128046_20190209');

//Convert the brightness temperature (Kelvin 0.1) to Celsius
var tempC = image.select("B10")
.multiply(0.1)
.subtract(273.5);

//Load and See the different of pixel value in raw landsat and tempC
Map.addLayer(tempC);

And here is my snippet script:
https://code.earthengine.google.com/c28efdb63bd2733d7fa27bf2e48e58e6

Specifically, the var image is only for one Landsat image. But I want to process for four Landsat images simultaneously. Is there any suggestion?

//Convert the brightness temperature (Kelvin 0.1) to Celsius
var tempC = image.select("B10")
.multiply(0.1)
.subtract(273.5);


texts: QGIS save layer as converts Layer Field Data Type

texts:
In QGIS I open a geodatabase and save one layer to a SHP file to work on.  I use the "Save Layer As" option with the default settings.  However I discovered that the SHP Field data type changes, in particular a "Date" type is saved to a "Text" type. I only discovered the issue when I was trying to use the field calculator and there were issues when using text in calculations.

Is there a setting I am missing that maintains the original data types.

Original Data Type from ArcCatalogue

QGIS data type


texts:
I created a shapefile of Sundarban area Bangladesh. Here it is:

Here is the projection information about it:

I want this shapefile to project upon an image which is here [I downloaded it from earthexplorer USGS where the path and row is 138,45 respectively]
while I was trying to project it using ERDAS software the shapefile does not appear.

Here is the image:

The image meta data is here:

and the shapefile metadata is here:

 .


texts: Creating feature layer in InMemoryWorkspace using ArcObjects?

texts:
I'm writing a WPF application, in C#, using ArcObjects.

I have an ESRI.ArcGIS.Controls.AxMapControl on my form, and I'm trying to draw some point features on top of it, in code.

It's been suggested that I create a point feature layer, in an InMemoryWorkspace, and add that to the map.

I've not been able to find any decent examples of how to do this, online.

Can anyone give me an example of how to create a point feature layer, in an InMemoryWorkspace, or point me to an example elsewhere?


texts:
I have a geojson file, but it's way too heavy for the app I need to develop

I'd like to know what alternative do I have to make it lighter.

I was thinking about simplifying it's shapes, and also turning to topojson.

Could you recommend me some tool to achieve this?

--

Edit: I saw the other thread stating how to convert to topojson, but I'd also like to simplify the shapes, they are much too detailed for the visualization I'm working on.


texts: How to define coordinate system and geographic transformation for a CAD drawing?

texts:
I am a GIS guy and I occasionally work with CAD files as well. Now I need to know how I can define a coordinate system and a geographic transformation for a CAD Drawing. According to ESRI:

"You can define a coordinate system and a geographic transformation and create and edit metadata for CAD datasets. The information is stored in external support files and moves with the CAD dataset when it is managed in the Catalog window."

But don't know how to do this! Can you please let me know or point me to a tutorial to understand this?


texts: simple raster format to use in Web applications?

texts:
I'm building a Web application in which users create and edit very coarse raster data. They're assigning values to cells in a georeferenced grid. Right now I'm storing the values in a multidimensional JavaScript array; for example:

[
  ['yes', 'no', 'maybe'],
  ['no', 'no', 'yes']
]

and I'm storing latitudes and longitudes for the bounding box, along with the number of rows and columns, as separate values on the overall object.

That works, but I'm just curious to know if there are any standard GIS data formats amenable to JavaScript I can use. In effect, my question is: shapefile is to GeoJSON as Esri grid is to what?


texts: Add MGRS custom Coordinate Reference System to QGIS

texts:
I am looking for a way to add MGRS custom Coordinate Reference System to QGIS.

How do I get the Proj.4 cartographic parameters for Military Grid Reference System (MGRS) to add via Settings >> Custom CRS?


texts: Calculate area of interest based on X events within Y kilometers within Z days

texts:
I'm currently in Afghanistan and need some geospatial support. I have a customer that needs to know where Areas of Interest (AOI) exist. AOIs are defined as 3 types of events within 1KM of each other within 30 days. I can figure out how to find X number of events within Y kilometers, but can't figure out the time part.  Any ideas?

Also, the 30 days does not mean within the last 30 days. An AOI may have existed back in 15May-15June, for example.

It is a vast amount of data, and I cannot share it, but I have an MGRS, lat/long in DD, and a date field.

I really appreciate any advice on how to tackle this. On the network I operate on, I only have the base install of ArcGIS 10, and the data is stored in a file geodatabase. I won't be able to install any other software such as PostGIS, but we do have python installed.

Thanks!


texts: Why can't a MapView be instantiated with a specific spatialReference? (ArcGIS Javascript API 4.4)

texts:
I'm attempting to set the spatialReference property of the view to 102667. The console logs the error:

  init.js:321 [esri.views.MapView] The view could not be initialized
  with the spatialReference 102667. Try specifying an extent or a center
  and scale

I'm using the sandbox environment as an example, because I'm not understanding why this wouldn't work. The data I am working with is all projected with this spatial reference, so I thought it would be prudent to set the View to match. The excerpt comes from here:

https://developers.arcgis.com/javascript/latest/sample-code/sandbox/index.html?sample=get-started-mapview

<!DOCTYPE html>

  <meta charset="utf-8">
  <meta name="viewport" content="initial-scale=1,maximum-scale=1,user-scalable=no">
  <title>Get started with MapView - Create a 2D map - 4.4</title>

  <link rel="stylesheet" href="https://js.arcgis.com/4.4/esri/css/main.css">


texts: Combining polygon features and retaining both geometries and attributes using `sf`

texts:
I use the sf R package regularly, but I'm having trouble finding the correct process to use in sf to get a parallel result to what I would get performing a union in ArcGIS. Given the image below, where I have a set of polygon features in purple, and another polygon feature shown by the heavy black outline, what I want is what I would get from using union in ArcGIS - that is, the purple polygon features sliced on the border of the black polygon, and for the purple features within the black polygon, the attributes would now also contain attributes of my black polygon.  For purple features outside the black polygon, they would have the attribute fields of the black polygon but they would have no values populated.

No matter how I parameterize st_union in sf, my interior purple polygons are dissolved with the bounding black polygon, and ALL features get all the attritues of the black polygon populated, even those outside the black polygon.  What would be the analogous process in sf to get the result one would have using union in ArcGIS? It's simply a combination of both polygon features I need, where purple polygons are sliced by the black polygon, and features inside get attributes of black polygon feature, while those outside get the attribute fields of the black polygon but un-populated  values.


texts: Print axMapControl contents with preview

texts:
I have a map with many feature layers on top of it. I would like to be able to print the map and feature layers either scaling it to fit a certain page size or printing it to actual scale thus creating a mosaic of printed pages.
Ideally I would like to preview the results before sending them to printer and I am having two major issues doing so.

What I did so far is get the map onto a axPageLayoutControl so that I can add artifacts about which part of the mosaic each page is, scalebars, titles etc. Although getting the axMapControl.ActiveView.FocusMap which appears as landscape, in the PageLayout (orientation set to landscape) the map appears in portait (actually showing a slightly different extend). Changing the pageLayout.page.orientation has no effect on the map. Also I tried drawing all the features from the axMapControl onto the pageLayout but the spread in a much bigger area than the map itself, but this time in landscape orientation.

Another problem is that I don't know how to preview the mosaic pages after setting the esriPageToPrinterMapping.esriPageMappingTile value to the axPageLayoutControl.Page.PageToPrinterMapping property as this does change the number of pages axPageLayoutControl.get_PrinterPageCount() to more than 1, but has no effect on the pageLayout or the printout itself.


texts:
I'm trying to automate the actual geo-referencing in PyQgis/PyQt5.
I went manually over the built-in GDAL Georeferencer in QGIS and saved a GCP-file. Also, the actual commands used for the conversion are available:

gdal_translate -of GTiff <gcpString> <inputFile> <tempFile>
gdalwarp -r cubic -order 2 -co COMPRESS=NONE <tempFile> <outputFile>

I tried to 'copy' this behaviour in PyQt5 using subprocess.Popen() but it seems from PyQgis I do not have permissions or access to the gdal cmd-tools.

Then, I stumbled upon the Python GDAL library which is (I think) built-in with QGIS, and seems a good option. However, I cannot get this to work with the given options as above and documentation is not very easy.
What would be the best approach for referencing the image with a given GCP-file and options?


texts: How to calculate area from only boundaries length without point coordinates

texts:
How can I calculate polygon area if I have only boundaries lengths?
Let check this simple example:

I have all boundaries length from this picture in the following format:

K6-K5: 24.15m
K5-K10: 10.54m
K10-K7: 23.63m
K7-K6: 12.29m

What is suitable approach for calculating polygon area?
I think we have only 2 ways:

Calculate current area of above plot, then scaling up
Draw a "real polygon" by boundaries lengths, then calculate area on "real polygon"

Is there another approach for getting more suitable and accurate way?


texts:
I'm trying to use shapely as simple cad kernel and have a problem with parallel_offset's robustness: Problem-GIST

When offsetting the geometry, it depends on the offset distance whether the transformation fails and I can not recognize any pattern.


texts: Is it possible to get the points used for the convex hull in osm2po?

texts:
I have a non-convex hull algorithm that I want to use for an app. With OSM2PO, it's possible to get a convexHull, but is there a way to get the points before the hull is constructed?

More, it would be awesome for me if I can get only the points that are "path end", I mean points that are part of segments that are connected to no more than one other segment.

(If you're intererrested by the non-convex hull algorithm, aka Alpha-shapes, I can share...)

Thanks for your quick answer, but I don't understand well your solution. I've seen this table, and it might be helpfull but I don't see how to query PGSQL (pass it millions indices to check ?) and my aim is to avoid PGSQL.

I'dont know how you dit it but when you traverse the graph, there's probably something that say "Hey stop traversing this branch, you've reached the max cost or this branch has no more points". It should be possible to grab these "last" points, don't you think ?

In order to be sure we understand well, here's a little example.

Imagine this "graph" :

                       7
                        /
                     6
                    /
1---2---3---5---8----9---10
            |
            4

router.getVisited() will give me 1,2,3,4,5,6,7,8,9,10.

My dream : only get 1,4,7,10

Hope you can help me. Many thanks again.

(sorry for answering my own question but too many chars to be posted as a comment...)


texts: rasterio rasterize does not take geometry

texts:
rasterize(allshapes)  return;

"ValueError: No valid geometry objects found for rasterize"

Testung a single shape:

print(shape.GetGeomFieldRef(0))
POINT (56.91288 27.33519)

print(is_valid_geom(shape.geometry()))
False

in rasterios github the function checks for following :

geom_type = geom['type']
    if geom_type not in geom_types.union({'GeometryCollection'}):
        return False

cheking for this by hand :

print(shape.geometry()['type'])
TypeError: 'Geometry' object is not subscriptable

and for :

if 'type' not in geom:
    return False`

`

cheking by Hand also returnes False

Changing to Polygon wont help

print(shape.geometry())
POLYGON ((....))

print(is_valid_geom(shape.geometry()))
False

I dont get this, pls help


texts: How to show a toolip over a feature with OpenLayers?

texts:
This is a continuation (part 2) of this previously asked GIS question.

I have an array of objects. Each object has a name, some score and some Well Known Text.

For each object, I'm displaying the WTK over an OpenLayers map.

I'm trying to have a tooltip over each shape, but only when the mouse hovers over that shape/feature.

if I have this code here..

var data = [
    { 'aaaa', 100, <some WKT>' },
    { 'bbbb', 32, <some WKT>' },
    { 'cccc', 93, <some WKT>' },
    { 'dddd', 46, <some WKT>' } ]

and for each object in that array I need to generate the feature.. which is more or less like this..

var feature = format.read(data.WellKnownText);
shapeLayer.addFeatures(feature);

var highlightControl = new OpenLayers.Control.SelectFeature(shapeLayer, {
    hover: true,
    highlightOnly: true,
    renderIntent: "temporary"
});

map.addControl(highlightControl);
highlightControl.activate();

I'm guessing I need to define something after the renderIntent? Like onhover call this function?

Don't forget, I also need to pass to that function data.Name and data.Score ...

Any suggestions, folks?


texts: Doing COUNTIFS equivalent in ArcGIS Field Calculator?

texts:
I have an attribute table with location data for a series of individual animals. They are often found at the same locations, so I have created a column in my attribute table where I list the rock ID. The Rock Id counts are repeated for each animal in order to keep data collection simple in the field, so rock 1 for animal 1 is not the same as rock 1 for animal 2. (In any instance where two animals used the same rock - I combined the rock id numbering system for their rocks, but this was rare.)  I would like to get a count for the number of times an animal is found at each rock, so I can use symbology on my map that varies in size based on this attribute. (More times an animal uses a rock, larger the symbol).

I have a data set that looks something like this and I want to populate rock count to look something like what I have below

I know I can use CountIFS in excel - and then potentially import this in, but Ideally I could calculate this for the data already in ArcGIS using code and the field calculator. The summarize tool works well for telling me the unique totals with two case fields- but I would like to replicate this value in each cell. The output of that table is smaller than my original, since obviously there are some repeats.

I don't think I can do a join from this summary table, because each possible joining attribute (Rock_ID, Animal_ID) has more than one row in both tables.


texts:
My task is to reproduce a map of a specific place. The problem are the borders which I can´t reproduce in QGIS.

The attached picture shows the line. It looks like a polygon filled with slashes.

Any ideas how I could reproduce this?


texts: Keep .gdb structure to postgresql + PostGIS (Feature dataset)

texts:
I have a .gdb file with a specific structure. I will extract all shp files for uploading to a postgtresql db. But I want to keep the structure of gdb. How can I keep this structure to postGIS+Postgresql? Is it impossible to keep the Feature dataset?

I use Qgis, ArcGIS 10.1, Postgresql and PostGIS 1.5.x


texts:
I'm hoping to analyze google timeline data in R for a final project in a GIS class,

Google can output timeline data as geojson and KML but it looks like the KML file includes a lot more info (type of transportation distance traveled etc) than the JSON file. Additionally, JSON is an option for the entire time lines but to download a single day/week/month  it looks like I have to use KML.

I understand from
How to efficiently read a kml file into R
and
Importing KML files to R
that I need to specify the layer info as well as the kml file name to readOGR(), what I'm a little confused about is exactly how the layer names are included in a kml file.

Looks like the <name> tag is associated with the layer name, but there are 122 name tags in the file so its clearly not exclusive to layer. Fine.

using layers <- ogrinfo(data_source) gets me

[1] "INFO: Open of `C:\\Users\\Documents\\GIS_Coursework_3\\history-2018-09-28.kml'"
[2] "      using driver `LIBKML' successful."
[3] "1:  Location history from 2018-09-28 to 2018-09-28 "

then using Location_History <- readOGR(data_source,  "Location history-2018-09-28 to 2018-09-28 ") gives:

Error in ogrInfo(dsn = dsn, layer = layer, encoding = encoding, use_iconv = use_iconv, :
Multiple incompatible geometries: wkbPoint: 12; wkbLineString: 12

The problem is that there are multiple "sub" layers,

in QGIS i can see that when I open the layer, there is a points layer and a lines layer

I dont' see either of these text strings in the KML files anywhere when I open it as a text file.

I could probably just copy those, but it's not particularly useful to me as a one off. I need to get the layer info programmatically, rather than opening QGIS every time.

Is it practical for me to start exploring xml parsing in R?

Is there a package I haven't been able to find that handles this stuff successfully?

Am I missing something obvious about how to read KML layer info?

This is the only feature I've found lacking in R compared to QGIS or ArcGIS. They've been pretty comparable so far which I've found impressive.


texts: UpdateCursor by unique values from another Field

texts:
I am trying to update a field in an attribute table in a feature Layer with arcpy.da.UpdateCursor.
I am running the Python script as part of a model in ArcGIS 10.6.
But I am lacking some of the basic concept of how arcpy.da.UpdateCursor edits the attribute table of my feature layer.
My Python script is based on https://community.esri.com/thread/237881-how-to-calculate-a-field-using-other-field:
import arcpy
import os

# function returns list of unique values in a given field of a table
def unique_values(table , field):
    with arcpy.da.SearchCursor(table, [field]) as cursor:
        return sorted({row[0] for row in cursor})

fc = os.path.abspath(arcpy.GetParameterAsText(0))
lst = unique_values(fc, 'NetworkID')
with arcpy.da.UpdateCursor(fc, ['NetworkID', 'outlet']) as cursor:
    for row in cursor:
        row[1] = lst.index(row[0])
        cursor.updateRow(row)

The field I want to edit is called 'outlet' and mostly consists of NULLs and it should be edited based on the value in 'NetworkID'. NetworkID has no NULLS and categorizes the outlets into groups (lets call them networkgroups).
How do I get all values in outlet of the same networkgroup to have a value based on the only not null value in outlet in the same networkgroup?
For example:
outlet NetworkID
null     1
null     1
null     1
S_21     1  #<- as this is the only non NULL in the networkgroup '1' I want to use it to define the other outlet in the same networkgroup '1'
null     2
S_45     2
...      ...

should look like:
outlet NetworkID
S_21     1
S_21     1
S_21     1
S_21     1
S_45     2
S_45     2
...      ...


texts: Where to find US National Park Service basemap as WMTS?

texts:
I'd like to locate the US National Park Service basemap, ideally as a WMS or WMTS to use in a leaflet map.
While browsing the NPS docs: https://www.nps.gov/maps/tools/npmap.js/api/#wms
I am repeatedly met with broken links of working examples. An on top of that, the docs keep using a current weather WMS in their documentation
(example: http://nowcoast.noaa.gov/wms/com.esri.wms.Esrimap/obs
However, it seems that this basemap does exist because the NPS is using it as their actual basemap in all their examples:

So where can I find the basemap itself?


texts:
I'm trying to teach myself OpenLayers and took an introductory course in which we received the whole code. For my first own project I want to do something similar and recycle some of the code. But there is some part of the code I don't understand and hope that someone can explain it to me. The idea is that the user can click on the vector element on the map (schools in this case) and recive information about its attributes. Here is the code:

map.on('singleclick', function(e) {
  var feature = map.forEachFeatureAtPixel(e.pixel, function(feature) {
    return feature;
  });

  var schoolName = feature ? feature.get('name') : '';
  alert(schoolName);
});

All right, when the user clicks on the map, it fires the function e. In this function we create a variable with the name feature and this variable recieves the function forEachFeatureAtPixel, which returns the variable feature. Thus, the way I understand it the value of the variable feature is itself. Thus, it shouldn't have a value. Where is my logic mistake?

Furthermore I don't understand the next line either. We create a variable schoolName and this variable gets the value "feature ? feature.get('name') : ' ';

But what does the "?" do and why do we need : ' ';

Can someone explain these lines of code?


texts:
I downloaded some Wikimapia data using the wikimapia.api. Everything is A-OK since I've opened this in Google Earth and everything is in place. However, I need to turn it into a shapefile, or get the data onto QGIS, since there are impt polygons.

Using GDAL Complete for OS X from kyngchaos, I try to run the ogr2ogr command.

$ ogr2ogr -f "ESRI Shapefile" mountains.shp api.wikimapia.org-4.kml
ERROR 4: Cannot open api.wikimapia.org-4.kml
FAILURE:
Unable to open datasource `api.wikimapia.org-4.kml' with the following drivers.

etcetc. These are named mountains which I wish to import onto QGIS.


texts: Visualize large amount of geopoints with Geoserver and PostGIS

texts:
I have a postgis database with over 50 million rows. Every row consists of a POINT corresponding to lon,lat and a value (number 1-100)

I'm trying to create tiles for zoom levels 5-20 that I want to represent on google maps with the help of GeoServer.

Depending on the value of the point I want to paint it in a different color. As many of points are close to each other I want to paint based on their combined average value. Obviously which points are combined is depending on the current zoom level. Does Geoserver/PostGIS provide a good solution to this problem?

My current solution is to use ST_SnapToGrid and use different tables for every zoom level with different grid sizes. A solution I'm not very happy with.


texts: Convert or add hand drawn geometries to FeatureCollection as they are drawn on GEE

texts:
Is there any way to convert hand-drawn geometries into a ee.FeatureCollection, ee.Geometry or something similar; in order to be used at the first instance and without having to re-run the script after drawing?

For example, the following code, slightly modified from one of the GEE examples (this is just a demonstration of the issue):

// Load an image.
var image = ee.Image('LANDSAT/LC08/C01/T1_TOA/LC08_044034_20140318');

// Define the visualization parameters.
var vizParams = {
  bands: ['B5', 'B4', 'B3'],
  min: 0,
  max: 0.5,
  gamma: [0.95, 1.1, 1]
};

// Center the map and display the image.
Map.setCenter(-122.1899, 37.5010, 10); // San Francisco Bay
Map.addLayer(image, vizParams, 'false color composite');

var clipButton = ui.Button({
    label: 'Clip image',
    onClick: function() {
      image = image.clip(geometry);
      Map.addLayer(image, vizParams, 'clipped');
    },
    style: {
        color: 'blue',
    },
});

print(clipButton)

If I press Run, later draw a polygon and then press the button the "geometry is not defined" error raises (Despite var geometry is added to imports section). In my case geometries cannot be defined before running the code because areas of interest are discovered after some calculations, so they are random and should be established by the user.

Summarizing, could hand-drawn geometries be recognized after running the code just once as some way and later be used through a button or other?


texts: Executing intersect algorithm with pyqgis?

texts:
I want to run the process "qgis:intersection" in a pyton-script and try it like this:

import processing
a1 = QgsMapLayerRegistry.instance().mapLayersByName('Layer1')[0]
a2 = QgsMapLayerRegistry.instance().mapLayersByName('Layer2')[0]
processing.runandload("qgis:intersection", a1, a2, "memory:myCut")

But there comes the "Problem executing algorithm" message: "GEOS geoprocessing error: One or more input features have invalid geometry."

If i do the same by the qgis-gui (Vector/Geoprocessing Tools/Intersect...) in QGIS 2.14.21 like this:

then everything is ok and qgis do the intersection.

Doing the intersectin with QGIS 2.18.13 with the gui, I become also an error like above.

"Layer2" is a large layer with some hundredthousands polygons an if i use the "Topology Checker" i get about 600 invalid geometries.

Why does the intersection work only with the qgis-gui in 2.14.21 and not somewhere else?


texts: Decrypt strange encoded characters in the attribute table of a shapefile

texts:
In QGIS, I can't read characters of the shapefile in my geological layer attribute table : I don't understand or there are misunderstandable characters of shapefile.

Below a screenshot of the attribute table.. Does anyone know how I could get the original values back ?


texts: How to compile GDAL with external libraries? (Rasterlite2)

texts:
I want to compile GDAL with Rasterlite2 Support, but I can't find any info on how to do this.

I already downloaded librasterlite2, libspatialite and libsqlite3.

When configuring with --with-rasterlite2=path, I only get an error:

checking for rl2_load_raw_tiles_into_dbms in -lrasterlite2... no
configure: error: did not find rl2_load_raw_tiles_into_dbms

Is the downloaded source from https://www.gaia-gis.it/fossil/librasterlite2/index wrong?

I'm on Windows and want to compile with mSys/MinGW


texts: Calculate attributes to table based on other attributes

texts:
I am going through thousands and thousands of road lines and giving them a value based on QC with imagery.  So for instance for the add_data field I give it one of these five numbers:

0 = not examined
1 = include
2 = consider excluding
3 = exclude
4 = DC

But for each of these I also need to fill out comments and type of road data.  So I was wondering if there is a way to calculate the comment and other fields based on the add_value field? So maybe it might look like this for the add data:

"1" = "include" comment = looks good  RD_Typ = Major Road.
  "2" = "include" comment = looks good  RD_Typ = Dirt Road.
  "3" = "include" comment = Joined data to other data  RD_Typ = Major Road.
  "4" = "include" comment =  Joined data to other data = Dirt Road.

..... and so on you get the point.  So the question is can I calculate the other fields based on the number I calculate in the add data field.  I know I could do this with a join and then calc field but that just does not seem to be as much fun.


texts: Adding delta symbol in field cell for shapefile in QGIS

texts: I'm just trying to add a delta symbol to label a feature for a shapefile in QGIS. I do this by typing the label in the appropriate cell during an editing session. I've tried copying and pasting Δ and it turns into a question mark for my label after saving the edits. I've also tried the ALT+30 method that produces the filled-in variety ▲. The same thing happens; it turns into a question mark after saving the edit. Thanks for your help.


texts: mapbox - disabling viewing of unrendered white tiles outside of the map

texts:
I'm looking to disable being able to move the map in such a way that one can move outside of the map area:

I don't want to be able to see the surrounding white areas when the map is dragged. I have turned off wrapping as I don't want any horizontal duplication.

Code:

    // Create a map in the div #map
    L.mapbox.map('map', 'example', {
        minZoom: 2,//need to change it to 1 maybe 0 for mobile devices
        // These options apply to the tile layer in the map.
        tileLayer: {
            // This map option disables world wrapping. by default, it is false.
            continuousWorld: false,
            // This option disables loading tiles outside of the world bounds.
            noWrap: true
        }
    });

I have a suspicion I might be able to utilise maxBounds or fitBounds but I'm no quite sure how to go about it.


texts: Can I use R to do a buffer inside polygons (shrink polygons / "negative" buffer)?

texts:
I want to use R (working with SpatialPolygons) to create a buffer inside my polygons - or actually shrink them by a value; ideally only if something is left after shrinking.
The buffer function I know (e.g. raster::buffer) allows only positive values (creating a buffer around features).
I thought of a workaround: Converting polygon to polylines (outline), creating buffer, subtracting buffer from original polygon. But that seems kind of dirty. I am wondering whether there is a package / function I can use.


texts: Driving directions that favor a particular highway

texts:
Given a set of geo coordinates and a geo line (KML file or GeoJSON data) is there an API or algorithm that returns a line that goes along that highway as much as possible while stopping at each coordinate on the way?

With most APIs the driving directions favor the shortest route between points.

What I'm looking for is a route that is not necessarily the shortest but that goes along a particular highway as much as possible.


texts:
I add WMS connection to openLayers, it has 127 different layers. I add one Layer which occurred to be small icons symbolizing Pictures on map. But the styling is ugly, according to the specification  of layers (GetCapabilities) "Styles" is set to "Default". I want to style each icon on this WMS layer, but i cant understand how.

My sample code below:

new TileLayer({
    id             : 5,
    name           : 'test picture_1',
    visible : false,
    source  : new TileWMS({
        url     : WMSUrl,
        params  : {
            'TILED'     : true,
            'LAYERS'    : 118
        },
        transparent: 'true',
        projection : "EPSG:3057"
    }),
}),

How could i style all features (which are small icons on the map) in this layer ?


texts: Adding, removing and calculating fields of existing vector layers using QGIS Modeler

texts:
I got a feeling I'm missing something very simple here, but how can I use the modeler in QGIS to do the following tasks:

Adding/Removing a field
Looping through a field
Calculate field values based on other fields

I can only find how to do new layers and work with them - but most of the time, I already have the layer loaded in QGIS - I just need to edit the fields. Processing modeler to update table suggests to use the Network plugin, but even then I can't get it to work.

My understanding was, that the layer has to be connected with the field calculator in order to work, like in QGIS 3 Graphical Modeler field calculator with attribute's parameters but when I choose a vector layer as Input in the model builder, I can't choose my existing layer.
How am I supposed to connect them then?
Until now, I used PyQGis to do these basic tasks, but my shapefile is too big and I was hoping the modeler loops through a big attribute table a bit faster.


texts: How to find multiple centers of gravity?

texts:
Say I have a number of customers in a geographical area, marked "c":

Each customer has a Lat/Long. I would like to find x number of points, in the same space, that best serve the customers, ie. they are closest to the customers.

Its easy to do when x = 1, because that is the geographical center of gravity, and there are a few ways to find that out. However I am struggling to think of how I would approach this problem if x > 1

Say x = 2, then I would like to know, where would I place the two points so the distances to customers are minimised, how would I even approach this problem?

I would eventually try and solve this using a programming language, so I am not after a purely mathematical formula, even though that would help.


texts: Produce a list of polygons that a line crosses

texts: I have a polyline feature class and a polygon feature class. What I want to do is add a field to the lines that lists all polygons that touch the line. This is effectively performing a select-by-location finding the polygons that touch each individual line then add the names of each of these to an attribute. Does anyone know how I can do this in ArcGIS Desktop or FME? My Python is not great but I can give it a go if that's the only way.


texts: Disabling line in Select by expression dialog of QGIS?

texts:
I would like to disable some lines in the 'Select by expression' dialog (QGIS 2.6) but I can't find whether it is possible or not, and if yes, how to achieve that.

I have tried several characters such as /, //, /*, #, REM.

Note that /* and REM transform the line into green text, which seems to be disabled, but the query result is empty and I get in invalid expression error below the query field.


texts: Importing grib, grib2, netcdf to PostGIS

texts:
According to grib2 to raster? (grib2 -> raster2pgsql -> PostGIS 2.0 -> GeoServer), it seems like that grib2 can be treated as raster directly. Is it the so called "standard" way?

How about grib and netcdf?


texts: Converting Feature Collection of one element (Multi Point Geometry) into Feature Collection of different elements

texts:
CODE:

var l8 = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
       .filterBounds(geometry);

// TEST IMAGE

 var first = ee.Image(l8.first()).clip(geometry)

// get image projection
var proj = first.select([0]).projection()

// get coordinates image
var latlon = ee.Image.pixelLonLat().reproject(proj)

// put each lon lat in a list
var coords = latlon.select(['longitude', 'latitude'])
             .reduceRegion({
 reducer: ee.Reducer.toList(),
 geometry: geometry,
 scale: 30
})

// get lat & lon
var lat = ee.List(coords.get('latitude'))
var lon = ee.List(coords.get('longitude'))

// zip them. Example: zip([1, 3],[2, 4]) --> [[1, 2], [3,4]]
var point_list = lon.zip(lat)
print(point_list);
// Create points
var mp = ee.FeatureCollection(ee.Geometry.MultiPoint(point_list));
print(mp);
Map.addLayer(mp,{}, 'All Points')
Map.centerObject(geometry);

Here when I print mp it is showing it as a Feature collection of a single element of MultiPoint geometry.

Example for feature collection of different elements:

var ran = ee.FeatureCollection.randomPoints(geometry,4516,10);
print(ran);
Map.addLayer(ran,{}, 'Random Points')

What I want it is the same as the above one which feature collection of different elements.

What changes should I do to convert my feature collection same as this one?

IMPORTS:

var geometry = /* color: #d63000 */ee.Geometry.Polygon(
    [[[-67.76435852050781, -22.213792811479802],
      [-67.76641845703125, -22.220785143658368],
      [-67.75646209716797, -22.228412744835992],
      [-67.74410247802734, -22.230319580297703],
      [-67.7420425415039, -22.221420792936158],
      [-67.74341583251953, -22.212203596469955],
      [-67.74787902832031, -22.207753698686023],
      [-67.75543212890625, -22.213474969917876]]]);


texts: Error attempting to Run ArcPy Script in Anaconda's Spyder

texts:
I'm attempting to switch over from running ArcPy scripts in Idle to running ArcPy scripts in Anaconda's Spyder GUI.

When I copy over the code and attempt to execute it, I get the following error:

File "C:\Program Files     (x86)\ArcGIS\Desktop10.3\arcpy\arcpy\geoprocessing\__init__.py", line 2, in <module>
from _base import *
ImportError: No module named '_base'

Here is the attached code:

import os
import sys
import string
import arcpy
from arcpy.sa import *
arcpy.CheckOutExtension("spatial")

workpath = "C:/workspace/gates/THRG/tmin_monthly_summaries/"
arcpy.env.workspace = workpath
tiffs = arcpy.ListRasters("*")
print (tiffs)


texts:
I'm using Mapfish Print Tool in Windows and everyhing works correctly. But when I print the same template on Linux server, Mapfish Print doesn't print all elements when printing on formats like png or jpg. It does work only with pdf.

Do you have any ideas, what's the problem?


texts: Geoserver - problem with MySQL support (impossible to create MySQL store at all)

texts:
I need to access spatial data in MySQL database with Geoserver. Official user manual states: just download package, extract it into WEB-INF/lib directory and all will be OK. Instead, just nothing happens - MySQL option doesn't appear when creating new data source, I can't find any relevant messages in logs and Google doesn't provide a hint.

What do I do:

Download MySQL extension, following recommended link. I have Geoserver version 2.4.3 (uses jre7, downloaded recently). Downloaded archive declares to be of version 2.4.3 (and manual specifically states that versions should match).
Unzip archive and copy it contents to [Geoserver installation directory]\webapps\geoserver\WEB-INF\lib (also tried other paths mentioned in the manual and extension readme, though they don't make any sense - result the same).
Restart Geoserver service (also tried rebooting system), than go to Data -> Stores -> Add new store in Geoserver web interface (logged in as admin). Manual states that MySQL option should be there, but it isn't.

As mentioned, I use Geoserver 2.4.3 (latest at current time), OS is Windows 7 x32. Geoserver is installed as a service and starts automatically.

Logging profile of server is set to VERBOSE_LOGGING, no traces of "mysql" are found. System logs don't contain anything but messages from Service Control Manager about service being started/stopped.
Trying to mess with other JARs in the same directory causes Geoserver to fail. There was installed Geoserver 2.1-beta3 some time ago, but current installation is made into different directory.

I'm interested in any advice how to get Geoserver-MySQL combination working, or how to diagnose why it is problematic/impossible.


texts: is it possible to place a label automatically outside its polygon?


texts: Impossible to Add geometry attributes on a .shp file with QGIS

texts:
I am using QGIS version 3.8.3-Zanzibar
I have a shapefile and I want to compute the area by using Vector/Geometry tool/Add Geometry attributes.

Runnig this tool I get :

Feature (11) has invalid geometry. Please fix the geometry or change the Processing setting to the "Ignore invalid input features" option.
Execution failed after 0.06 seconds

Any suggestion how to solve this?


texts: Select only raster overlapping shape files using Rasterio

texts:
I'm cutting out several polygons from different rasters using the IDs of a shapefile and I would like the images to be just those polygons that overlap the raster without saving the empty images. Rasterio already stops as soon as it encounters an ID outside the image, returning the message:

  ValueError: Input shapes do not overlap raster.

How can I solve this?

My code:

for polygon in polygons:
    with fiona.open(polygon,'r')as shp:
        shape=[feature["geometry"] for feature in shp]

        for raster_path in rasterfiles:
            output_path = outputPath + 'ID_{}'.format(shp[0]) + '.tif'

            with rasterio.open(raster_path) as tile:
                out_image, out_transform = rasterio.mask.mask(tile, shape, crop=True)
                out_meta = tile.meta

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

                with rasterio.open(output_path, "w", **out_meta) as dest:
                    dest.write(out_image)


texts: Copying dbf data into shapefile using ArcPy?

texts:
How do I perform a one-to-many join in a python script, from a dbf to a shapefile?

I saw (Greg Krakow's method) using QGIS, but I can't use it (I'm working on it, but QGIS is not yet recognize as an "official and usable" option within our organisation).

My dbf contains many records with the same ID and my source shapefile contains the geometry for each ID.

My approach is to create from scratch a shapefile ("outFile" in the script) by copying the structure of my dbf (adding field by field), using arcpy.FieldList()

dbfFieldList = arcpy.ListFields(dbf)

for field in dbfFieldList:
    fieldTypeShp = validFieldType(field.type)
    arcpy.AddField_management (outFile, field.name, fieldTypeShp, "", "", field.length)

validFieldType(field.type) is a function "converting" the field type of my dbf to a valid field type for arcpy. For example, I have a "smallint" which came back as "SHORT" so I could use AddField_management...

I figured that if I could read my dbf line by line, I could, field by field (since both my dbf and shp have the exact same structure), udpate the value of my shapefile.

#-- Seach cursor for my .dbf
searchCur = arcpy.SearchCursor(dbf)

for sCur in searchCur:

    #-- Insert cursor for my shapefile
    cur = arcpy.InsertCursor(outFile)
    row = cur.newRow()

    #-- Reading for the current row, the fields value, one by one.
    for field in dbfFieldList:
        fieldName = field.name
        val = sCur.getValue(fieldName)
        row.setValue(fieldName,val)

    #-- Inserting the new row into my shapefile
    cur.InsertRow(row)

But I have an error executing the function "setValue". If I remove that line and the cur.InsertRow (or put a "print val" it works perfectly. I have created a shapefile with the same number of row as my dbf, but I can't add the field values.

After that part, I would use the ID field to retrieve the geometry from my source shapefile with the same ID, recreating my 1-to-many join.

If only I could make that part work.


texts: Reclass raster with range input in Model Builder (QGIS)

texts:
I would like to use the QGIS model builder to preform reclassifications on several rasters. In QGIS I use r.reclass to do the reclassification by a given range (f.e.: 100 thru 150 = 1)

In Model Builder I'm not able to set the range-input (Sunshine duration range in the screenshot) to the r.reclass rule. As you can see in the screenshot. Of course, r.reclass is expecting a rule like "100 thru 150 = 1".

How could I do a reclassification by a range input on a raster using the model builder?


texts: Wrong Coloramp function in ArcGIS 10.2.1

texts:
I want make elevation map use new method Learn a new method

But i have problem with Coloramp function whatever the choice of color ramps.

What's wrong?


texts: Looking for a web service for change detection

texts: I'm looking for service that provides me a change detection service based on a database of remotely sensed data. Has anybody happen to know such a service? The service I am looking for is a geoportal that alarms the user of a specific change (like oil spill or major coastal change or even fire in forest) in a specific area or globally. It can be either freely available service or on-demand service.


texts: ModelBuilder Copy Features tool not updating shape when run as geoprocessing service

texts:
I have a model that uses the "copy features" tool and writes down a shapefile. When I run it as a tool, the shapefile is created. But when I export the tool as a service, if invoke the service the shapefile is not updated.

If I delete the shapefile created when I ran the tool, the service invocation does not create the shapefile either. No errors, no logs, no nothing. Maybe I'm missing some conceptual stuff here.

It seems that the image service is not available :(

f25_solucion = point geometry
departamentos  = polygons

The whole idea is to get all points that fall within certain polygons. The model runs fine, but runs only when ran as a tool.


texts: JSON To Features tool does not obey symbology?

texts:
I am using the JSON to features tool in the ArcMap conversion toolbox and i cannot seem to get the symbology right.

This JSON i am using is a GPFeatureRecordSetLayer which i have written to a .json file

{
    "displayFieldName": "pointFeature",
    "hasM": false,
    "hasZ": false,
    "geometryType": "esriGeometryPoint",
    "spatialReference": {
        "wkid": 102100,
        "latestWkid": 3857
    },
    "fields": [
        {
            "name": "OBJECTID",
            "type": "esriFieldTypeOID",
            "alias": "OBJECTID"
        },
        ...
    ],
    "features": [
        {
            "geometry": {
                "type": "point",
                "x": 19865753.933378372,
                "y": -2051854.6857393598,
                "spatialReference": {
                    "wkid": 102100,
                    "latestWkid": 3857
                }
            },
            "symbol": {
                "color": {
                    "r": 0,
                    "g": 0,
                    "b": 0,
                    "a": 0.001176470588235294
                },
                "size": 13.333333333333332,
                "type": "simplemarkersymbol",
                "style": "circle",
                "outline": {
                    "color": {
                        "r": 0,
                        "g": 255,
                        "b": 255,
                        "a": 1
                    },
                    "width": 2.6666666666666665,
                    "type": "simplelinesymbol",
                    "style": "solid"
                },
                "xoffset": 0,
                "yoffset": 0
            },
            "attributes": {
                "OBJECTID": "4131",
                ...

            }
        }
    ],
    "exceededTransferLimit": true
}

When i run the Features to JSON tool on the generated features i get this JSON string which totally skips any symbology information.

{
    "displayFieldName": "pointFeatures",
    "fieldAliases": {
        "OBJECTID": "OBJECTID"
    },
    "geometryType": "esriGeometryPoint",
    "spatialReference": {
        "wkid": 102100,
        "latestWkid": 3587
    },
    "fields": [
        {
            "name": "OBJECTID",
            "type": "esriFieldTypeOID",
            "alias": "OBJECTID"
        },
        ...
    ],
    "features": [
        {
            "attributes": {
                "OBJECTID": "4131",
                ...
            },
            "geometry": {
                "x": 19865753.933378372,
                "y": -2051854.6857393598
            }
        }
    ]
}

Everything else works fine. Is this a known limitation?


texts:
I am using the GeoExt example and this is the js code and trying to get rid of the multiple lat and lon that get added into the PopUp on multiple clicks.

Not sure which is the code that I need to comment out? I thought that removing,

this.handlerOptions = OpenLayers.Util.extend(
            options && options.handlerOptions || {},
            this.defaultHandlerOptions
        );

Would do the job but it seems that is not right as I am still getting the multiple popups and commenting out that code made no difference. Not sure what that code was for?


texts: I am somewhat new to programming and the SQLite world. I am looking for a way to take data from and SQLite DB and populate feature classes. The feature classes are all polygons and the intent is to get their shape from the Well Known Text field within the SQLite DB. I also need to populate two other fields in the feature class where the data is coming from the SQLite DB.


texts: Creating empty DBF table and populate it with PyQGIS

texts:
I need to create a standalone DBF file with some features. I know how to add the fields but I'm not able to add features and to create only a standalone DBF, without any other files such as .shp or .prj. I've tried creating a LineString in memory and then writing it with QgsVectorFileWriter, using ESRI Shapefile as driver, but it creates a simple shapefile with no features. The GDAL docs say that it can manage standalone DBF but I'm struggling with it and I can't find the way to do it.
This is the code I tried to create it and the export it.
self.work_line_table = QgsVectorLayer('LineString', 'Line_table', 'memory')

[Code for adding fields]

# Add features using the attribute from another layer. It's not working
self.work_line_table.startEditing()
for line in self.work_line_layer.getFeatures():
    line_id = line['IdLinia']
    feature = QgsFeature()
    feature.setAttributes([line_id, self.municipi_codi_ine])
    self.work_line_table.dataProvider().addFeatures([feature])

self.work_line_table.commitChanges()

# Export the table
QgsVectorFileWriter.writeAsVectorFormat(self.work_line_table,
                                        os.path.join(GENERADOR_WORK_DIR, 'MM_LiniesTaula.dbf'),
                                        'utf-8', self.crs, 'ESRI Shapefile')

I've thought in creating it as shapefile and then removing all the other files associated to it except the DBF file, but anyway I can't add any feature to the DBF.


texts:
What's a provider in PyQGIS and how many types of providers exist?
I am really not sure what a provider does. There is an example below from PyQGIS cookbook where a provider is used to load a QgsVectorLayer.
#get the path to the shapefile e.g. /home/project/data/ports.shp
path_to_airports_layer = "testdata/airports.shp"

#The format is:
vlayer = QgsVectorLayer(data_source, layer_name, provider_name)

vlayer = QgsVectorLayer(path_to_airports_layer, "Airports layer", "ogr")


texts: Can't compile plugin builder plugin in QGIS 2.4 (Windows)

texts:
I'm trying to build a basic plugin using plugin builder but I can't even compile it. I tried following a few tutorials and just trying to compile whatever but that always results in the same problem.

Error messages for rcc4.exe and uic4.bin:

"C:\Program Files\QGIS Chugiak\bin\pyrcc4" -o my_plugin.py my_plugin_dialog_base.ui
No resources in resource description.

"C:\Program Files\QGIS Chugiak\bin\pyuic4" -o my_plugin.py my_plugin_dialog_base.ui
ImportError: No module named site

relevant files in the plugin folder:

my_plugin.py
my_plugin_dialog.py
my_plugin_dialog_base.ui
resources.qrc

What's in the /bin folder:

pyrcc4.exe
pyuic4.bat
pyuic4.bat.tmpl

Any idea what's wrong? I'm limited to using whatever's installed with QGIS to compile this, by the way.


texts: Using ArcPy to find all points that lie within certain radius of other points?

texts:
I am trying to find a way to count all the points that are only within a specified distance of 500 meters of another set of points.

My first shapefile contains a series of points of various types of businesses in an area. And the other shapefile contains only points of churches in the same area. I want to overlay the two shapefiles programmatically using python to determine how many of my original set of points lie within a 500 radius of the church points.

I've attempted to solve this is but I am running into a few problems. First is I am not sure how to extract a particular type of business from my shapefile, say barbershops, I am not sure how I would specify to python that I only want to find the number of barbershops that lie within a radius 500 Meters of all Churches in the area.

Here is what I am attempting:

import arcpy

# Unsure what exactly Buffer_analysis does, so not sure if this would work
arypy.Buffer_analysis("LocalBusinesses.shp", "C:\Users\new\pointDataBuffer.shp", "500 Meters", "#", "#", "NONE")

target_features = "C:\Users\new\churches.shp"

join_features = "C:\Users\new\LocalBusinesses.shp"

out_features = "C:\Users\new\pointDataBuffer.shp"

arcpy.SpatialJoin_analysis(target_features, join_features, out_features,"#","#","#","#","500 Meters","#")

Problems I am facing are to specific that I only want barbershops in the first shapefile. I am not sure if 500 Meters there will take care of that.

Is there a more efficient way of performing this?

Wasn't sure how to use this info, but it seems relevant:

Finding points that are within a set distance of other points
Finding points a certain distance from each other

Also wondering if this might work:

http://resources.arcgis.com/en/help/main/10.2/#/Point_Distance/00080000001r000000/


texts:
GEOSVoronoiDiagram is not published in the Shapely API. However, it is available in the newer GEOS libraries. I'm trying to use it with

func = getattr(geos._lgeos, 'GEOSVoronoiDiagram_r')
attr = ftools.partial(func, geos.lgeos.geos_handle)
attr.__name__ = func.__name__
setattr(geos.lgeos, 'GEOSVoronoiDiagram', attr)
geos.lgeos.GEOSVoronoiDiagram.restype = c_void_p
geos.lgeos.GEOSVoronoiDiagram.argtypes = [c_void_p, c_void_p, c_double, c_int]

gc = geom_factory(
    geos.lgeos.GEOSVoronoiDiagram(
        geom._geom,
        env._geom,
        tolerance,
        int(edges)
    )
)
g for g in gc.geoms]

Which is copied from relevant code in Shapely, but that gives the error

ctypes.ArgumentError: argument 4: <type 'exceptions.TypeError'>: Don't know how to convert parameter 4

What am I missing?


texts: Error reported when trying to import prj for data points in New Brunswick Canada into QGIS

texts:
I have data points recorded with a Trimble GIS. The projection is New Brunswick Sterographic Double Projection, which as near as I can find via a Google search is EPSG:2953. However, when I try to import the proj (see below) I get the error "not a valid or recognised data source". The same files import just fine into ARCMAP 10. How do I get these files to read correctly into QGIS 1.7.3 (Windows 7) or does this mean there is a bug in the projection data base of QGIS? I guess I could drop the projection and specify it in the program but I am concerned that the projection may be wrong if the files don't import correctly in the first place. Thanks

Alyre

Projection file:

PROJCS["NAD_1983_CSRS_New_Brunswick_Stereographic",GEOGCS["GCS_North_American_1983_CSRS",DATUM["D_North_American_1983_CSRS",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Double_Stereographic"],PARAMETER["False_Easting",2500000.0],PARAMETER["False_Northing",7500000.0],PARAMETER["Central_Meridian",-66.5],PARAMETER["Scale_Factor",0.999912],PARAMETER["Latitude_Of_Origin",46.5],UNIT["Meter",1.0]]


texts: Is PostGis suitable for Global data routing?

texts:
I'm searching quite some time now...

Is it possible to load the global routing data into PostGis?
What can I expect from a performance point of view?
It's for routing purposes only.

I've downloaded the global data, extracted it, saw the 300Gig size and started wondering...   the .osm is 300Gig, but still has to be converted to pgRouting.

Searching on the internet ends up mostly in examples whereas pgRouting is used for boundaries like cities...  hence my question over here.

Thanks for answering, any type of help is useful!

Johan


texts: Using findAndReplaceWorkspacePaths to change datasource for CAD files- ArcPy

texts:
I try to use "findAndReplaceWorkspacePaths" in order to change datasource for CAD files. This code works fine but when there no cad layers to work with, despite this the code will print 'findAndReplaceWorkspacePaths'. How can i get this print message only when the code will find cad layer?

import arcpy,os,sys,string
from arcpy import env

env.workspace = r"G:\PROJECTS\daras\gis"
for mxdname in arcpy.ListFiles("*.mxd"):
    print mxdname # print list of mxd's in the folder
    mxd = arcpy.mapping.MapDocument(r"G:\PROJECTS\daras\gis\\" + mxdname)
    mxd.findAndReplaceWorkspacePaths(r"E:\daras\gis\layers", r"G:\PROJECTS\daras\gis\layers")
    print 'findAndReplaceWorkspacePaths'
    mxd.save()
del mxd

this is the result:

airplane noise2 a3.mxd
findAndReplaceWorkspacePaths
environment.mxd
findAndReplaceWorkspacePaths
environment011.mxd
findAndReplaceWorkspacePaths
environment012.mxd
findAndReplaceWorkspacePaths


texts: How to use an IF THEN statement in an sql Select expression?

texts:
I need to select the top 3 records from a table, but depending on an inline variable from an iteration, occasionally one of those top 3 records should not be selected and should be replaced by the top 4th.

I'm using sql expression in a personal geodatabase with Arc.

So far I have this:

SELECT * FROM MyTable WHERE
[Field1] in(SELECT TOP 3 ([Field1])
FROM MyTable
ORDER BY [Field1] DESC)
AND NOT [Field2] = '%Variable%'

This just selects the top three or the top 2 when [Field2] equals an inline variable. I'm a syntax sql NOOB, but I have a basic understanding of what can be done and that it is probably a simple solution.

Maybe:

If
Count([Field1] in(SELECT TOP 3 ([Field1])
FROM MyTable
ORDER BY [Field1] DESC)
AND NOT [Field2] = '%Variable%') = 3
SELECT * FROM MyTable WHERE
[Field1] in(SELECT TOP 3 ([Field1])
FROM MyTable
ORDER BY [Field1] DESC)
AND NOT [Field2] = '%Variable%'
END
Else
[Field1] in(SELECT TOP 4 ([Field1])
FROM MyTable
ORDER BY [Field1] DESC)
AND NOT [Field2] = '%Variable%'
END

Something like that?


texts: Reusing leaflet map elements in custom functions in Shiny/R

texts:
I have an application that reuses several leaflet elements, so I'm trying to build functions to streamline my code.

All of my logical attempts failed, so I utilized an example I was familiar with here. I copied the format and created functions for markers and polygons, seen in the example below:

library(leaflet)
library(shiny)

markerFunction <- function(map, x, y){
  addCircleMarkers(
    map, x, y,
    color = "red"
  )
}

polygonFunction <- function(map, polygon){
  addPolygons(
    map, polygon,
    color = "red"
  )
}

# shinyapp
shinyApp(
  # ui
  ui = fixedPage(
    # leaflet output
    leafletOutput("map", height = 600)
  ),  # END UI

  # server
  server = function(input, output, session){

    # map output
    output$map <- renderLeaflet({
      leaflet() %>%
        addTiles
    })  # END RENDER LEAFLET

    # map proxy
    proxy <- leafletProxy("map")

    # observe
    observe({
      proxy %>%
        markerFunction(40, 50) %>%
        polygonFunction(gadmCHE)
    })
  }  # END SERVER
)

The marker function works fine, but the polygon function spits the error:

  Warning: Error in derivePolygons: addPolygons must be called with both
  lng and lat, or with neither.

which is the error I get when I don't include data = explicitly. When I add data = polygon to the function definition and data = gadmCHE to the polygonFunction() call, I get the error:

  Warning: Error in polygonFunction: unused argument (data = gadmCHE)

Running the markerFunction() on its own yields the expected result, but I'm perplexed how it works without including the map argument. When I remove it from the function, it no longer works.

I've tried many configurations and can't get any of them to work except for the markerFunction() as it exists in the above example app. I have tried including and excluding the map argument, explicitly defining other arguments (like x, y, and data) in both the function definition and function calls, excluding the explicit arguments with no luck. I want to create similar functions for polygons as well, but haven't found a solution.


texts:
Has anyone had luck adding text to a button in Leaflet Draw? I would like to add text next to the icons on the buttons I have, but there doesn't seem a way to specify that information in the documentation.

I have the following code to customize and add the drawing button:

drawControlFull = new L.Control.Draw({
    position: 'topright',
    draw: {
        polygon: false,
        polyline: false,
        circle: false,
        marker: false,
        circlemarker: false,
        rectangle: {
            shapeOptions: {
                color: '#1976D2'
            },
            repeatMode: false,
        },
    },
    edit: false,
});

L.drawLocal.draw.toolbar.buttons.rectangle = 'Draw Bounding Box';

map.addControl(drawControlFull);

I notice that this L.drawLocal.draw.toolbar.buttons.rectangle = 'Draw Bounding Box'; generates a span that has a class of sr-only and therefore the text doesn't appear.

I have tried modifying the .innerHTML of the a tag for .leaflet-draw-draw-rectangle using
$(".leaflet-draw-draw-rectangle").html("Draw");

but that also doesn't work. Any ideas? Trying to get something like this: