texts: Dividing a trapezium into quarters with QGIS?

texts:
Using QGIS 2.18.3, I need to divide a trapezium (Math is Fun says that this is called a trapezoid in the UK - in any event, a 4-sided polygon with no sides parallel) into 4 polygons.  Please refer to my crude MS Paint approximation to visualize what I need.  I suspect that this would be accomplished by finding the midpoint along each input side and connecting it with the midpoint of the opposite side, thereby creating quarters.

Purpose: I have polygonal vector data of surveyed USA land at the section (square mile) level where each polygon is usually square, but often containing trapezium-shaped sections, that I need to divide into approximate quarter-sections.  Since square polygons would exist in the data, the ability to subdivide a square would thus also be a requirement.

The ability to accomplish this subdivision with multiple sections at one time would be a bonus.

A search through the plugins didn't turn up anything, nor did a Google search.  Sadly, I have no programming skills to accomplish this.  Suggestions?


texts: Adding feature class name automatically from overlay layer into attribute of clipped layer while iterating clip

texts:
I am trying to clip the road layer to district boundary layer in QGIS using clip and enabling 'iterate' option and I get the clipped road layers of individual districts. Later, I need to add the district names into the road layer attributes under field 'dist_name' manually. In order to add the dist_name I am opening the road layer and entering the district name using field calculator one by one.

Is there any way how I can add it automatically doing the clipping process?


texts: Leaflet - Hiding labels when any popup window is shown?

texts:
I am new to Leaflet and have found the following solution to add static labels to features in my map:

L.CircleMarker.include({
    bindLabel: function (content, options) {
        if (!this._label || this._label.options !== options) {
            this._label = new L.Label(options, this);
        }

        this._label.setContent(content);
        this._labelNoHide = options && options.noHide;

        if (!this._showLabelAdded) {
            if (this._labelNoHide) {
                this
                    .on('remove', this.hideLabel, this)
                    .on('move', this._moveLabel, this);
                this._showLabel({latlng: this.getLatLng()});
            } else {
                this
                    .on('mouseover', this._showLabel, this)
                    .on('mousemove', this._moveLabel, this)
                    .on('mouseout remove', this._hideLabel, this);
                if (L.Browser.touch) {
                    this.on('click', this._showLabel, this);
                }
            }
            this._showLabelAdded = true;
        }

        return this;
    },

    unbindLabel: function () {
        if (this._label) {
            this._hideLabel();
            this._label = null;
            this._showLabelAdded = false;
            if (this._labelNoHide) {
                this
                    .off('remove', this._hideLabel, this)
                    .off('move', this._moveLabel, this);
            } else {
                this
                    .off('mouseover', this._showLabel, this)
                    .off('mousemove', this._moveLabel, this)
                    .off('mouseout remove', this._hideLabel, this);
            }
        }
        return this;
    }
});

L.marker([-37.7772, 175.2606]).bindLabel('Look revealing label!').addTo(map);

The thing is that if the label is too big the labels are displayed above the popup window, not behind.

I don't know how to hide those labels when the popup appears?


texts: Publishing a web tile layer on ArcGIS Online with python

texts:
I have a tiff layer in an ArcGIS Pro project. I would like to publish it on my ArcGIS Online portal, with Arcpy. I try to publish it as a Web Tile Layer, but could be another layer type, if anyone has another suggestion.

Previously I published it manually from ArGIS Pro, as TPK. Then I updated the TPK and the Tile Layer from my python script.
However now, I would like to directly create the Web Tile Layer with my script - that is the layer does not yet exist on the portal.

First I created an ArcGIS Pro project. I add the tiff layer to the project. This tiff layer will change each year and the script will publish this new layer as a new web layer on the portal.

from arcgis import features
from arcgis import mapping
import arcpy
import sys

# Connection to AGOL
gis = agol_connect()

dir = 'K:/WEB-SIG/PROD/URBAINE/'
agglo = 'LILLE'
year = '2018'
pol='PM25'
metrique='moy'

project_name ="TEST_%s_%s" % (agglo,year)
project_path = dir + r"ARCGIS_projects\%s\%s.aprx" % (project_name,project_name)
map_name = agglo+'_'+ pol +'_'+ metrique

# Parameters for Arcgis Pro
arcpy.env.overwriteOutput = True
arcpy.env.workspace = dir + "ARCGIS_projects"

# Connection to the project in ArcGIS pro
project = arcpy.mp.ArcGISProject(project_path)
map = project.listMaps(map_name)[0]
layer = map.listLayers()[0] #this is the tiff layer

Now I tried to publish it but I didn’t succeed ..

First attempt :
Using the WebMap function of mapping module. But it does not work with tile Layer.

wm = mapping.WebMap()
wm.add_layer(layer)

Second attempt : I tried to create a TileSharingDraft and to publish it. I got the following error.

# Create TileSharingDraft and set service properties
service = "LILLE_2018_PM25_moy"
sharing_draft = map.getWebLayerSharingDraft("HOSTING_SERVER", "TILE", "TileDraftSharing", layer)
sharing_draft.summary = service
sharing_draft.tags = "modele"
sharing_draft.description = "Modelisation urbaine"
sharing_draft.credits = "Atmo Hauts-de-France"

# Create Service Definition Draft file
draft_output = dir + "WebLayer/" + service + ".ssdraft"
sharing_draft.exportToSDDraft(draft_output)

# Stage Service
sd_filename = service + ".sd"
sd_output =  dir + "WebLayer/" + sd_filename
arcpy.StageService_server(draft_output, sd_output)

# Share to portal
print("Uploading Service Definition...")
arcpy.UploadServiceDefinition_server(sd_output, "My Hosted Services")

   However this way seems a bit complicated for the simple publication I want to do. I supposed that there is a more direct way to do it, as for TPK, a bit like this :

#with the arcGIS Pro Map
arcpy.CreateMapTileLayer_management(map,layer_name_online)
#or with the layer
arcpy.CreateMapTileLayer_management(layer, layer_name_online)

Does anyone have an idea?


texts:
I have two shapefiles. One is the DMA and another of states.

What I'm trying to do is "expand" a DMA boundary to match the boundary of the states, since the DMA map doesn't extend all the way. To illustrate:

This is the map of the Mid Atlantic and Midwest region states:

And this is the DMA for the same region:

And this is the overlay between the two shapefiles:

Is there any way I can create a query that will artifically increase the size of the DMA to match the state boundaries? The reason why I need this is because I plan on going a large spatial join with coordinates and I know alot of people will not be counted for if we use the current DMA shapefile.


texts: Converting annotation into a line feature?

texts:
This might seem like an odd question, but the end-result will be a GPX file that is to be used in delineating areas in a Mapathon, in which the editor only uses GPX files.

My question is:

Does anyone know how to convert my annotation layer representing grid names inside polygons to a feature class, representing the text as lines?
Meaning; the text has to be "written" with lines that contains actual spatial information.
My hope is to create a feature class that contains the lines representing the text, as it would look in my annotation layer, but made of lines with coordinates. The image below only shows the annotations within my gridcells.

I am using ArcMAP, but Qgis and other potential softwares could be included. Please feel free to ask, if you want me to elaborate the problem statement any further.


texts:
Is it possible to input a zipped .kml file (.kmz) int Leaflet?

I was trying the code:

fetch('Cov.kmz')
      .then( res => res.text() )
      .then( kmltext => {

            // Create new kml overlay
            parser = new DOMParser();
            kml = parser.parseFromString(kmltext,"text/xml");

            console.log(kml)

            const track = new L.KML(kml)
            map.addLayer(track)

            // Adjust map to show the kml
            const bounds = track.getBounds()
            map.fitBounds( bounds )

      })

But it doesn't appear on the map at all.

There is the plugin for that:

https://github.com/engrabid411/kmzol3

which doesn't work for the time being.

I also found an advice here:

https://stackoverflow.com/questions/57147510/kml-layer-fetch-error

But to be honest, extracting .kmz file into .kml and another ones, comprising of the .kmz file not help me enough.


texts: Label size changes when exporting or printing WMTS layers

texts:
I am using QGIS2.4.0 on a Windows computer. When I export a map containing OSM landscape, the labels in the output map are smaller then in the view in the composer. (even if the view in the composer is at zoom 100%)
This is not the only difference, here an example (click to magnify images):
The view in the composer

The exported file

You can see that the labels are now impossible to read.
I want to be able to print and export the view as I see it in the composer.
I tried the solutions for the question How to improve OSM layer resolution in OpenLayers plugin while exporting? but that didn't help.
So the problem is not the openlayers plugin.


texts: Computing the list of POI "visited" by a given track

texts:
The problem: finding "visited POIs"
(1) I have on the one hand a set X of let's say Points of Interest (POI) defined by (label, latitude, longitude), for instance
Pall = {("Punab", 45.18287, 5.71691), ("RailwayStation", 45.19145, 5.71507) ... }.

I've on the other hand a track T1 generated by a mobile app. It is represented as a GPX file.
I want to compute (in a java server) the sequence PVisited of points that were "visited" by the track T1. That is I would like to implement a java function "Visited" like that (100m is the precision):
Visited(Pall,T1,100)=PVisited

(2) As a matter of clarification here is an example of map with (I'm not using google map in the project, this is just to provide an example) :

one track T1 (orange),
a set of POI (red),
some of them being "visited" (in yellow).

Let's say that "visited" means being 100m from the track. Just because it may matter for the sake of precision and approximation I know that the area of all POIs is 20km x 20km.
In practice, I'm not using Google Map. I just want to upload the gpx file in a java server, match the track against the set of POIs (in a database if possible) and add a flag for each visited point, with ideally, the approximate time of the visit. Something like :
PVisited = [
   ("Punab", "Saturday, March 27, 2021 5:10 PM CET"),
   ("Caserne', "Saturday, March 27, 2021 5:15 PM CET"),
   ("Hoche", "Saturday, March 27, 2021 5:25 PM CET")}

I want to build a java server that allow uploading GPX tracks Tx and for each track Tx compute the set PVistedByTx.
What I've tried
I'm a newbie. In the last weeks I tried many different tracking applications and just realize  that none of them was providing this feature. So this may not be so simple. Anyway I want to integrate this "computed visited POI" features on my java server.
I tried to look at GPX libraries but so far I didn't find something useful. But I may have missed something.
global or local ? I'm wondering if the global function Visited is available somehow (taking all POIs at once, and the full track at once), or if I need to iterate with a local function like this (iterating on each POI and on each track segment) :
DistancePS(POI, Segment) : Integer

I've seen various formulae for Distance(P1,P2) but I gues that what I need is the distance between a Point and a Segment.
converting to (x,y) plan and regular geometry ? my plan was to convert everything in (x,y) coordinates and then do some geometry computing. I've tried to convert my POIs to (x,y) coordinates using different formulas. Some are quite complicated. I guess that with having a 20x20km square and something like 100m could help having something rather simple. But solutions I've tried (for instance this) does not seem to work. I don't know if the "convert to (x,y) + geometry" plan is the good one.
Conclusion
I would like to know if a (simply) library provide the function
Visited(ps:Set(Point),Track) : List(Point)

If not what is there a (simple) library that provide the function
DistancePS(p:Point,s:Segment) : Real

Or What is the formulae to compute this ? Since the POIs are in a 20x20km square and a precision  of 100m, may be something simple could work ?
I do not want to embed a GIS into the server. I need something as light as possible since this is the only geographical-related feature I need.
I thought computing this was a piece of cake, then it is note. I've searched and now I'm lost. I need someone to show me where north is.


texts: Vertical datum conversion in ArcGIS Pro not recognizing vertical coordinate system

texts:
I have a GeoTiff surface model and the z values reference EGM96.  I would like to use the vertical datum conversion capabilities in ArcGIS Pro 2.7 to convert the GeoTiff to reference the NAVD88 vertical datum.
I see that the Project Raster tool includes a checkbox for applying vertical transformations. This checkbox disappears when I select my input raster.  After reading the tool documentation.
I understand the vertical option is only available if and only if the data have vertical datum information.
Thinking that the data must lack vertical information read that I can add vertical datum information using the Define Projection tool.  Trouble is, when I add my surface model to the Define Projection tool the Define Projection tool states that the data have vertical information.  See the images below.  Thinking that there must be something about the GeoTiff that Pro is not recognizing I ran the Define Projection despite the warning on the tool.
After running "Define Projection" the Project Raster tool still removed the vertical option once I added the GeoTiff suggesting the tool is again not seeing the vertical information.
I see this from the tool help:

To perform a vertical transformation, check the optional Vertical parameter on the dialog box. By default, the Vertical parameter is unavailable and is only available when the input and output coordinate systems have a vertical coordinate system (VCS), and the input feature class coordinates have z-values. Also, additional data (coordinate systems data) setup must be installed on the system.

What is this additional coordinate system data come from and where is it installed?  In short, how do I apply a vertical transformation to these data?


texts: ogr2ogr TopologyException invalid input geometry; ring self-intersection

texts:
I'm facing a problem when using osgeo/gdal:ubuntu-small-3.2.1 docker image and running these commands on the Natural Earth vector dataset:
# Get and unzip the data
# (inspired from https://github.com/mapbox/natural-earth-tm2/blob/master/util/setup_db.sh)
wget -q -N http://naciscdn.org/naturalearth/packages/natural_earth_vector.sqlite.zip
unzip -o -d natural_earth_vector natural_earth_vector.sqlite.zip

# Then, insert data into PG using ogr2ogr:
ogr2ogr \
    -progress \
    -f Postgresql \
    -s_srs EPSG:4326 \
    -t_srs EPSG:3857 \
    -clipsrc -180 -85.05 180.1 85.05 \
    PG:"host=localhost port=5432 dbname=mydb user=postgres password=****" \
    -lco GEOMETRY_NAME=geometry \
    -lco OVERWRITE=YES \
    -lco DIM=2 \
    -nlt GEOMETRY \
    -overwrite \
    "natural_earth_vector/packages/natural_earth_vector.sqlite"

The PostgreSQL/PostGIS database I'm trying to inject this data into is based on postgis/postgis:11-3.1  docker image.
This is the error I have:
0ERROR 1: TopologyException: Input geom 0 is invalid: Ring Self-intersection at or near point 30.796820538582722 61.704133379158407 at 30.796820538582722 61.704133379158407

This seems to be the culprit watershed in that particular case.
But this error actually appears more than 1'000 thousands times on other geometries    .
I'm not sure where to search for a fix. Any ideas?
If this helps, here is the GDAL version number from inside the container:
# ogr2ogr --version
GDAL 3.2.1, released 2020/12/29

This thread sounds related but it's from R: Getting TopologyException: Input geom 1 is invalid which is due to self-intersection in R? and I don't know, if it's the same error, how could I transpose the solution here.


texts:
I'm working on a project where we would like to add QGIS features to an existing Qt application.

The Qt application has been built with Qt5.2.1 using MinGW 4.8 and thus I think I need to build QGIS with MinGW 4.8. If I'm wrong on this I'd be glad to hear what can be done.

If I'm right, my problem is that the QGIS building-from-source tutorial (https://raw.github.com/qgis/QGIS/master/INSTALL) is indeed outdated as far as MinGW is concerned and the given msys environnement is not working with the current QGIS cmake files. I tried to edit those but some libraries are not in msys (spatialite and spatialindex) and it seems linking to Qt5 instead of Qt4 isn't working (problem with QT4_WRAP_CPP).

I thought I'd just try to build spatialite myself first but the procedure seems to take ages (http://www.gaia-gis.it/gaia-sins/mingw_how_to.html) and I'm not even sure it would work in the end.

If there's anything else I can do (a QGIS plugin for Qt for instance) I'd forget about building without an inch of regret so do not hesitate to share. The only thing to consider is that I'd rather keep Qt5.2.1 because it includes OpenGL and Android support which are both probably going to be needed for further developments.

In case nothing can be done for QGIS to work with my project, do you know any other GIS API that has comparable features ?

Thanks for your help !


texts: Converting Fiona Collection to Geopandas GeoDataFrame?

texts:
I have a feature class in a file geodatabase that I attempted to open with Geopandas, but I got the "list index out of range" error.  To bypass this error, I'm reading the feature class directly with Fiona.  This worked without a problem, and now I have a Fiona collection.  What I can't figure out, is how to continue processing this file in Geopandas.  I have been unable to find anything in the respective documentations about how to process an existing object.

I also submitted a bug issue at github/geopandas.


texts: ArcGIS Pro- no option to add graphic layer

texts:
According to the documentation, it should be easy to draw graphics onto a map....
"On the Map tab, in the Layer group, click Add Graphics Layer.
A new graphics layer appears in the Contents pane."
In my Map tab- Layer Group, I only have 'Basemap, Add Data, Add Preset"
I have search for commands, ect but simply can't find a solution.
Any suggestions?


texts: Getting information from one layer in attribute table of another layer

texts:
In the picture you can see polygon number 39a in the center. This polygon is drawn on one layer. The layer underneath holds information with soil conditions (green and red). You can see that polygon 39a is laying in 'vochtig zand'(red) and a little bit in 'Nat zandleem'(green).

In the 39a polygon attribute table, I want to see the soil condition which is most present within the surface of the polygon. Not by making a new 'clipping' shapefile but by a formula or expression. Is this possible?
Can anybody help me with the expression formula?


texts: Esri-Leaflet getFeature returns 'undefined'

texts:
I'm using Esri-Leaflet and whenever I try to use getFeature() the console.log returns "undefined".
There are no other errors to accompany the undefined message.
var service = L.esri.featureLayer({url: 'https://services1.arcgis.com/example/arcgis/rest/services/LayerName/FeatureServer/0'});
var objID = parseInt("4");
var featurePoly = service.getFeature(objID);
console.log(service.featurePoly);

A feature of ObjectID of 4 does exist. The feature service is reachable.
Any thoughts on where I might be going wrong?


texts: Switching between background 64 bit processing and 32 (so I can use 3rd party ie XTools)

texts:
I often run scripts now with 64 bit background processing enabled but when I want to use certain tools like Xtools table to text I have to turn it off as these tools are 32 bit. Is there a way to switch this within a script say... (Geoprocessing_Enable is fictional)

#start at 64 then change to 32 for Xtools
arcpy.Geoprocessing_Enable(32)

arcpy.XToolsPro_Table2Text(in, out, fields)

#go back to 64 after 32 bit tool ran successfully
arcpy.Geoprocessing_Enable(64)

Is there anything equivalent to doing this?


texts: Merge Vectorfile Outputs from Model Batchprocess in QGIS

texts:
I have written a model in QGIS 3.6 that allows to vectorize the outlines of raster layers containing measurement data. I run this model as a batch process to vectorize several field outlines in one go.

In the end, I want to have a single file containing my field outlines. The Batchprocess, however, creates one output for each input. Until now I just store all these output files (output_1.gpkg, output_2.gpkg,..., output_X.gpkg) in a temporary folder on my desktop and use the SAGA Merge Vector Layers to merge all the outputs into one file.

I am now trying to put this last step directly into the script, but since the outputs are the result of the script itself and only created in the batchmode, there is no file I can hand to the merge tool as an input. Is there a way to use a path to the temporary folder as an input in the Modeler?

I use QGIS 3.6 on a Windows 7 PC.


texts: Make imagecollection from list of images through a for loop in Earth Engine Python API

texts:
I try to code Mosaicking a Image Collection by Date (day) in Google Earth Engine this issue with Python.

#Load Sentinel 1 and filter data
def load_dataset(ImageCollection_ID,begin,end,aoi):
    ic = ee.ImageCollection(ImageCollection_ID).filterDate(begin,end).filterBounds(aoi)
    return ic

def filter_sentinel1(ImageCollection,polarisation,instrumentMode,resolution):
    ic = ImageCollection.filter(ee.Filter.listContains('transmitterReceiverPolarisation',polarisation)).filter(ee.Filter.eq('instrumentMode',instrumentMode)).filterMetadata('resolution_meters','equals', resolution)
    return ic

def seperate_look_angels(ImageCollection,polarisation):
    Ascending = ImageCollection.filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING')).select(polarisation)
    Descending = ImageCollection.filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING')).select(polarisation)
    return Ascending,Descending

def get_properties(ImageCollection):
    features = ImageCollection.getInfo()['features']
    dict_list = []
    for f in features:
        prop = f['properties']
        dict_list.append(prop)
    df = pd.DataFrame.from_records(dict_list).drop(['system:footprint','transmitterReceiverPolarisation'],axis=1)
    #Pandas Series of unique distinc values in df
    unique = df.nunique()
    im_id_list = [item.get('id') for item in ImageCollection.getInfo().get('features')]
    date_list = [datetime.datetime.strptime(x[35:43],'%Y%m%d') for x in im_id_list]
    #property_names = list(df.columns.values)
    return unique, im_id_list, date_list

def make_mosaic(date):
    date = ee.Date(date['value'])
    filterCollection = VV_Ascending.filterDate(date, date.advance(1,'day'))
    image = ee.Image(filterCollection.mosaic()).copyProperties(filterCollection.first(),["system:time_start"])
    return image

#Time of interest
begin = ee.Date.fromYMD(2016,1,1)
end = ee.Date.fromYMD(2016,3,1)
date_range = end.difference(begin, 'day')

#Source dataset
ried_225_222 = ee.FeatureCollection('users/tillmueller1990/ried_225_222')
sentinel1 = load_dataset('COPERNICUS/S1_GRD',begin,end,ried_225_222)
#Filter dataset for High resolution and Vertical transmitt vertical receive
sentinel1_VV = filter_sentinel1(sentinel1,'VV','IW',10)
#Filter for different look angles
VV_Ascending,VV_Descending = seperate_look_angels(sentinel1_VV,'VV')

#Get list of ids,dates and unique count of prop
unique, im_id_list, date_list = get_properties(VV_Ascending)
date_list = ee.List([ee.Date(x) for x in date_list])
newList = ee.List([])

for date in date_list.getInfo():
    mosaic = ee.Image(make_mosaic(date))
    print(mosaic.getInfo())
    newList.add(mosaic)
newcol = ee.ImageCollection(newList)
print(newcol.getInfo())

the function make_mosaic(date) return for every date a image. The output of print(mosaic.getInfo()) is

{'type': 'Image', 'bands': [{'id': 'VV', 'data_type': {'type': 'PixelType', 'precision': 'double'}, 'crs': 'EPSG:4326', 'crs_transform': [1.0, 0.0, 0.0, 0.0, 1.0, 0.0]}], 'properties': {'system:time_start': 1452618462260.0}}

I try to add this image to an empty ee.List([]) named newList to call later ee.ImageCollection(newList) but the newList remains empty. Why i cant add the images to the list with newList.add(mosaic)?


texts: Saving ArcGIS Pro symbology to PostGIS DB for each layer

texts:
I want to have a default style/symbology assigned to each layer in a PostGIS database so that whenever it is loaded in ArcGIS Pro, it will come up with the default style.
There is a way in QGIS to save styles to PostGIS DB so that whenever the layer is dropped to the canvas it will come with the default style, but the style saved in DB will be a qml file and hence will not work with ArcGIS Pro.
I am looking for a way to add style to DB such that they load in ArcGIS Pro as most of the users for the tool I am developing will be using ArcGIS Pro to connect to the DB.
Is this possible?


texts: Tracking progress of running geoprocessing tools in ArcPy?

texts:
I am looking for a way to access a GP tool running in a Python script to be able to see how much of the job has already been done.

I have found SetProgressor function, but this one is available only when running script tools via ArcGIS Desktop GUI. I am in contrary interested in obtaining the information on the progress via the code. This is to be able to track on the execution of the tool (this will be reported to a user which has triggered running the code). I don't run the code inside the ArcGIS Desktop application session.

So, for instance, when calling the Buffer GP tool on a polygon feature class (via arcpy.analysis.Buffer), I want to get an estimate from the tool after several seconds how much (in %) of the work was already done and this information should be possible access via the code.


texts:
I have only one shapefile with multiple polygons. Every week I make updates to this shapefile, including new polygons, changing the shape of others and deleting other polygons, that is, this shapefile is always changing.
When preparing a layout for this shapefile, I define a geographic scope that fits all polygons and also define a scale.

As this shapefile changes frequently, I have to go back layout and again adjust the geographic scope and scale again, with each change I have to revise the layout.

How could I do it in QGIS so that the scope and scale is updated automatically, as the shapefile has changes?


texts: Show both label and icon for amenities using Carto CSS

texts:
I'm trying to style, e.g. museums, so that both their icon and their name shows up. However, there seems to be some conflict. Here's my - related - code:

.points {
    ::tourism {
        [tourism = 'museum'][zoom >= 17] {
            marker-placement: interior; marker-fill: @icon_fill;
            marker-file: url('icons/museum.svg');
            marker-clip: false;
        }
    }
}

#area_label[zoom>=10] {
    [type = 'museum'][zoom >= 17] {
        text-name: "[name]";
        text-fill: #111;
        text-size: 8;
        text-wrap-width: 30;
        text-line-spacing: -1.5;
        text-face-name: @sans;
        text-placement: interior;
    }
}

Now, the thing is:

If I leave only the icon, it looks fine
If I leave only the label, it still looks fine
If I leave both, most likely only the label shows (though at times, it's only the icon that shows)
If I add marker-allow-overlap, it becomes a mess

I'm looking at the "official" OpenStreetMap theme, and it seems to be handling it great: At lower zoom levels, the label may not show up, but if you zoom enough (my case), both show. (To be precise, the icon is on top, and below you can find the label).

How do they achieve that? Any hints?


texts: Control scale dependent visibility of symbols created using geometry generator in QGIS

texts:
I am using QGIS v3.2 "Bonn".
I currently have symbols that have labels and call-out lines created using the geometry generator. I can make the labels follow scale dependent visibility rules, and would like to have the call-out lines follow the same rules but haven't found a way to do so.
Is there a way to set scale dependent visibility for such objects? Below is the geometry generator interface.


texts:
i'm newbie in GIS.
i'm trying to make webgis using leaflet.
i already make layer with icon in leaflet.

/*layer*/
    var layer = new L.GeoJSON.AJAX("GIS/layer/request_layer.php", {
        middleware:function(data){
        return L.geoJson(data, {
            onEachFeature: function(feature, layer){
                layer.bindPopup("<center>" + feature.properties.id + "<center>"),
                layer.setIcon(icon);
            }
        }).addTo(map);
        }
    });

/*icon layer*/
    var icon = L.icon({
            iconUrl : 'GIS/leaflet/images/icon.png',
            iconSize : [15,15],
            popupAnchor : [0, -28]
        });

how to add buffer in this layer ?


texts: Getting data type of each column from PostGIS table?

texts:
I need to get the column data types of all columns in a table, including the geometry types. What I wanted to know is if there is a function or SQL that gives something like this:

column_name | data_type
------------+--------------
gid         | integer
descr       | character varying(32)
class       | character varying(10)
area        | double precision
geom        | geometry(Polygon,3763)

From a few answers on stackexchange and gis.stackexchange I know that I can get some of the information with the following query:

SELECT
    g.column_name,
    g.data_type,
    g.character_maximum_length,
    g.udt_name,
    f.type,
    f.srid
FROM
     information_schema.columns as g JOIN
     geometry_columns AS f
         ON (g.table_schema = f.f_table_schema and g.table_name = f.f_table_name )
WHERE
    table_schema = 'my_schema_name' and
    table_name = 'my_table_name'

Result:

column_name | data_type         | character_maximum_length | udt_name | type    | srid
------------+-------------------+--------------------------+----------+---------+------
gid         | integer           |                          |          |         |
descr       | character varying | 32                       |          |         |
class       | character varying | 10                       |          |         |
area        | double precision  |                          |
geom        | USER-DEFINED      |                          | geometry | Polygon | 3763

But, is there a proper more practical way to retrieve the information in the format I need? Or must I enter the "world" of CASE WHEN structures and string concatenation to gather all column attributes in one single column in that format?

My fear is if a non expected data type surprises me by needing other attribute from the information_schema.columns table. I.e, in the example table before, I did not used any numeric (15,2) data type, that would need to use another attributes (numeric_precision and numeric_scale) to be parsed by a CASE WHEN.


texts: Add Item and Publish Service with ArcREST or ArcPy

texts:
I would like to have an automated process which does the following:

Adds a CSV to our portaL/AGOL organization account.
Publishes the item as a hosted feature service.
Consumes hosted feature back unto EGIS environment.

I have checked out ArcREST as a possible solution, however, one hang up that I can imagine is the recurring ID number that will change as the added CSV is updated on a daily basis via task scheduler.

Is this solution possible with ArcREST?

I will post my code sample when I return to my desk.

I have used a piece of your code Kevin in addition to the addItem method from ArcREST.  My item adds, and my token and item ID print, however, when I try to pass the itemID to the function below it does not publish and no errors are returned.

import arcrest
import json
import requests
import urllib
import urllib2
if __name__ == "__main__":
    username = "myusername"
    password = "mypassword"
    portalId = "myportalnumber"
    url = "http://myorg.maps.arcgis.com/"
    thumbnail_url = ''
    file = "myflie.csv"

    securityHandler = arcrest.AGOLTokenSecurityHandler(username,
                                                       password)
    #   Create the administration connection
    #
    admin = arcrest.manageorg.Administration(url, securityHandler)
    #   Connect to the portal
    #
    portal = admin.portals(portalId)
    #   Access the content properties to add the item
    #
    content = admin.content
    #   Provide the item parameters
    #

    itemParams = arcrest.manageorg.ItemParameter()
    itemParams.thumbnailurl = thumbnail_url
    itemParams.title = "CSVTest"
    itemParams.type = "CSV"
    itemParams.tags = "test"
    #   Enter in the username you wish to load the item to
    #
    usercontent = content.usercontent(username=username)
    # print usercontent.addItem(filePath=file,
    #                           itemParameters=itemParams, overwrite=True)
    agolID= usercontent.addItem(itemParameters=itemParams)['id']

    print agolID

    itemId = agolID
    tags = "Demo, Publishing"
    securityHandler = arcrest.AGOLTokenSecurityHandler(username,
                                                       password)

### Generate Token ###
gtUrl = 'https://www.arcgis.com/sharing/rest/generateToken'
gtValues = {'username' : 'myuser',
'password' : 'mypass',
'referer' : 'http://www.arcgis.com',
'f' : 'json' }
gtData = urllib.urlencode(gtValues)
gtRequest = urllib2.Request(gtUrl, gtData)
gtResponse = urllib2.urlopen(gtRequest)
gtJson = json.load(gtResponse)
token = gtJson['token']

print token

def PublishService(Short, agolID, username, token, itemName, publishParams, XField, YField):
    ''' Publishes the input itemID (uploaded CSV) this is using
    a post request from urllib, the JSON input can be generated by
    the analyze function. Item name, X and Y fields are specified'''

    publishURL = 'http://myOrg.maps.arcgis.com/sharing/rest/content/users/myUser/publish'.format(Short, username)
    publishParams['CSVTest'] = itemName
    publishParams['locationType'] = 'coordinates'
    publishParams['Y_COR'] = YField
    publishParams['X_COR'] = XField
    query_dict = {
        'itemID': agolID,
        'filetype': 'csv',
        'f': 'json',
        'token': token,
        'publishParameters':publishParams}

    query_dict['publishParameters'] = json.dumps(query_dict['publishParameters'], sort_keys=False)
    request = requests.post(publishURL, data=query_dict)
    print request.json()`


texts: Importing 3D DXF file into QGIS with z information as attribute value?

texts:
I have a dxf file of 3d contours. after import it does have only  layer,x,y attribute and doesn't have Z attribute.

How can I provide Z attribute to create a TIN file?


texts:
adding markers works fine using OpenLayers 2 but now I need to add markers using OpenLayers 4, how can I do this? my code as follows:

var markers = <?php echo json_encode($data)?>;

    var points = new OpenLayers.Layer.Markers( "Markers" );

    for (var i = 0; i < markers.length; i++) {

    var lon = markers[i][0];
    var lat = markers[i][1];

        var lonLat = new OpenLayers.LonLat(lon,lat)
        .transform(
        new OpenLayers.Projection("EPSG:4326"), // transform from WGS 1984
        map.getProjectionObject() // to Spherical Mercator Projection

        );

    points.addMarker(new OpenLayers.Marker(lonLat));

    };

    map.addLayer(points);

    var poiBounds = points.getExtent();
    var lonLat = poiBounds.getCenterLonLat()
    map.setCenter (lonLat, 8)


texts: Plotting shapefile into Google Maps and then plot some of the points

texts:
I am currently doing spatial analysis project where I have a shapefile (in Delhi region). That shapefile also consists of ward level data. Sample data set is given below:

A    B    C   Ward_No
1    2    3     6
4    5    6    12
43   4    82   26
..    ...  ....

I want to plot that file into my Google map (based of Ward_No) and then plot some points on it.

So far I have done the following parts in R:

# Reading a shapefile with rgdal
library(rgdal)
library(ggplot2)
library(tmap)
s1 <- readOGR(".", "Delhi")
summary(s1)

# Convert it into data frame
y = as.data.frame(s1)

# Plotting that shapefile... Not in top of Google map
z = tm_shape(s1) +
  tm_borders("black") +
  tm_layout(frame = FALSE,
            inner.margins = c(0, 0.25, 0, 0.01))

# Create a data frame which having lat/lon data

# No     lon         lat
# 1   77.138015   28.6516444
# 2   77.0738028  28.6387532
# 3   77.1337863  28.7146212
# 4   77.0714699  28.5897667
# 5   77.0919553  28.626676
# 6   77.2860315  28.6320635

# Plot coordinates above the shapefile.

pt <- SpatialPointsDataFrame(coords = geodata[,c("lon","lat")], data = geodata,proj4string = CRS("+proj=longlat +datum=WGS84"))
z = tm_shape(s1) +
  tm_borders("black") +
  tm_layout(frame = FALSE,
            inner.margins = c(0, 0.25, 0, 0.01)) +
  tm_shape(pt) +
  tm_dots(size=1, col="yellow", border.col="black")

It's working fine but when I try to plot this in my Google map it's not working. The code is:

sgmap <- get_map(location="New Delhi", zoom=11,
                 source="google")

p <- ggmap(sgmap) + z
p


texts: Is it possible to connect to an Access db via ModelBuilder?

texts:
I created an OLE DB Connection in ArcCatalog Database Connections. If I right click on one of my db tables I can "Create a Feature Class > From XY Table" and continue with my data processing.

However, if I try to incorporate this step into a geoprocessing model (ModelBuilder) I can't seem to access the tables through my .odc connection. See image.

Is there another way to those Access tables within ModelBuilder, or am I doing this wrong way?

MS Access version is 2010.ArcGIS version is 10.1. I connected to .accdb.


texts: How to simplify borders between polygons in shapely?

texts:
I have some fairly complex, and not completely accurate shapes stored as shapely Polygon objects. I'm trying to find the common borders, which should be as easy as finding the LineString returned from poly1.intersection(poly2).

However, due the the slight inaccuracies I'm finding that intersection is returning Polygon and Point objects, or sometimes a GeometryCollection object. Never a simple LineString except in my simple test code where the Polygon shapes are accurate.

Is there any way of simplifying the borders, like a "snap-to" method, or a tolerance parameter that might ease my pain?


texts: Configurable Map Viewer (CMV) - Configure dynamic layer to not display all sub-layers?

texts:
Using the example below, from the CMV, what would be the best way to not show the Fire Stations layer in the Louisville Public Safety map service assuming I wanted to see everything else and I wanted sublayertoggle true?

            type: 'dynamic',
            url: 'http://sampleserver1.arcgisonline.com/ArcGIS/rest/services/PublicSafety/PublicSafetyOperationalLayers/MapServer',
            title: 'Louisville Public Safety',
            slider: true,
            noLegend: false,
            collapsed: false,
            sublayerToggle: true, //true to automatically turn on sublayers
            options: {
                id: 'louisvillePubSafety',
                opacity: 1.0,
                visible: true,
                imageParameters: imageParameters
            },
            identifyLayerInfos: {
                layerIds: [2, 4, 5, 8, 12, 21]
            }

I've added these 3 lines after the declaring and setting of imageParameters:

var imageParameters2 = new ImageParameters();
    imageParameters2.layerIds = [1];
    imageParameters2.layerOption = ImageParameters.LAYER_OPTION_EXCLUDE;

and changed imageParameters: imageParameters to imageParameters: imageParameters2 in the layer config but the Fire Station layer still shows up.


texts:
I'm trying to create a BalloonPopup that has 2 different buttons or interactive elements, but I don't seem to manage a way to do this. The guide I'm following is here https://carto.com/developers/mobile-sdk/guides/vector-objects-on-map/#balloonpopup in the "Add a BalloonPopup" section.

I did manage to detect onClick events on the BalloonPopup setting a setVectorElementEventListener in the associated Layer. But I don't find a way of detecting which element inside of the BalloonPopup was clicked.

Basically what I want is to have 2 different buttons inside the Balloon, and then detect which of them has been clicked, so I can perform 2 different actions.

I could open an Android floating dialog, but it won't look as great as this would.

Any help will be more than welcome.


texts: Import GPX file to QGIS and correct apparent CRS problems

texts: I am using Easy Trails software on an iPhone to output a gpx file of a walk. When I import the file (import vector layer) to QGIS, into a project with OSGB information and the correct CRS for OSGB data, the imported track appears well to the West and South of where it should be. I have tried changing the CRS in both layers, tried changing the project CRS, enabled and disabled on the fly projection but no combination of what I try will get the track into the right place. I feel I must be doing something fundamentally wrong but cannot work out what it might be. Am I missing some vital step out?


texts: Correlating number values with text values from different fields in ArcMap?

texts:
I need to create a vector layer for representing the geology of an area. I have the geological map, a raster file.

I made some polygons and I named them in the attribute table, but I  don`t want to spend time writing the name of the rocks associated with each polygon. I made another field named code and I want to correlate those codes with the polygon names, so when I introduce that code, the name associated appears.

I think that I can do this in Field Calculator, but I don`t know how.

Is there any possibility to do this ?


texts: New Shapefile (Polygon) Disappearing in Edit Session ArcGIS 10.6.1

texts:
I created a shapefile in ArcCatalog and assigned it its coordinate system as I always do. I selected "polygon" as the type.

I added a generic topographic basemap to ArcMap.

Upon editing the new shapefile, if I hit save edits, the entire thing disappears and is not in the attribute table either. This has happened 4x now with new shapefiles, using a different basemap, and having nothing on but a basemap and new file.

Has anyone experienced this before or can pinpoint what I am doing wrong? I made a new polygon shapefile last month and did not have this problem! I have rebooted my computer as well to see if that was a glitch.


texts: Error in supervised classification output (Earth Engine)

texts:
I have a problem in the supervised classification output. The classification itself is giving any problem, but when I try to print or display the image I obtain this error:

  Image (Error)
  Unable to use a collection in an algorithm that requires a feature or image. This may happen when trying to use a collection of collections where a collection of features is expected; use flatten, or map a function to convert inner collections to features. Use clipToCollection (instead of clip) to clip an image to a collection.

Here is the code:
https://code.earthengine.google.com/9dae0afc574251bd102e81e10850dd82

I'm pretty sure that the problem is from line 103 to 128. Here I created a polygon from a raster, and set it as a FeatureCollection in order to set a class WATER with value = 1.
If I use polygons created by giving the coordinates of nodes, I don't have this problem in the classification. I'll attach another link where I use polygons and don't have any problem:
https://code.earthengine.google.com/6b03174ebd1c75a9281944c4a17c78b3

The problem is that I need to use the polygon from the raster and I have no clue how to solve this problem.


texts: Raster Clip becomes pixelated in ArcGIS for Desktop?

texts:
I'm running ArcGIS 10.1 (Basic License). I have a large aerial imagery for a region that I need to clip to smaller polygon. I have 28 different polygons that need to be cut out from the one aerial imagery. All but 5 area worked perfectly. I'm using the clip tool found in Data management Tools > Raster > Raster Processing > Clip. I have a definition query on my polygon shapefile (but I've clipped the others with the definition query on).

Below is a picture of the area that I want to clip from my aerial (in red).

This is what my resultant clip looks like:

Does anyone have an idea of where I might have gone wrong with my clip?


texts:
I am looking for a general approach to "refreshing" a NETCDF store using the REST API.

I have a NETCDF file that is generated every 24 hours.

I need to update the file, but would prefer not having to recreate each layer, as there are many.

I thought perhaps I could rename the coverage store and replace the backing file, then name it back to the original coverage store name, but I am unsure if that will trigger geoserver to reprocess the contained information.


texts: Loading raster results in wrong spatial location

texts:
I have read all sorts of entries about CRS inconsistencies of different layers but it just didnt help to solve my problem.

So far I have always used ArcGIS and recently changed to QGIS. My data consists of a DEM (find attached: 'rangsdorf.ovr') and several sounding points (find one example attached: 'Layer1_Verlust_mit_Maechtigkeit.shp'). All those layers have the same coordinate system (UTM) and looking at the data in ArcGIS everything is visualised the way I want it to be - the points are located within the DEM.

However, loading it into QGIS, the raster is located somewhere completly different - I really dont get it. As I said, it should all have the same coordinate reference system, and anyway I am using OTF so I´m absolutly clueless.

Link to the data: https://www.dropbox.com/sh/31xyddr78qklsnf/AAA68tY06D-B4TISvaMX602Aa?dl=0


texts: Unfolding polygons and other feature classes in Antarctica-like fashion using ArcGIS Desktop?

texts:
Is there a way in ArcGIS Desktop 10.4 to 'unfold' (reshape/modify) feature classes, e.g. a polygon of a circular city, to make them look similar to a Mercator or plate carrée projection of Antarctica, i.e. opened out into a strip at the bottom of a map? Distortions away from the feature of interest are fine.

This is my specific example:
I have a polygon of an almost circular city and other feature classes showing transport links into the city's periphery. Rather than producing a map with the empty polygon of the city in the centre of the map, I would like to unfold the city's polygon and turn it into a strip at the bottom of the map with all the other feature classes extending upwards from this strip. I want the city and its periphery to look similar to how Antarctica is displayed on maps with a Mercator or plate carrée projection.
This question is not about changing the projection of a data file.


texts:
I published around 1000 layers in geoserver without enabling tile caching for layer.

How to check a cache is enabled for a layer or not?

How to enable tile caching (check box) for a layer under Tile Caching tab using python or Curl?


texts:
I have this real Big Bandet Image with 25000 Bands.
I want to Sum over the Bands, Pixelwise.

is there a Way to either :

-convert this ee.image, bandwise to an ee.ImageCollection

-Collapse all Bands of means into 1 Summed Band ?


texts: Counting number of pixels with negative NDVI value within polygon using Google Earth Engine?

texts:
I have made a map of vegetation loss and gain for a polygon in Google Earth Engine using NDVI and would like to calculate the area where loss occurred. I attempted to use the reducer to get the sum of pixels but I somehow need to specify the range of pixels to count only negative values where loss occurred and then total in square feet. My code is posted below, bottom part is where I'm having difficulty.

Can I use the stats.get function to only return negative values somehow?

var collection = ee.ImageCollection('LANDSAT/LE07/C01/T1_8DAY_NDVI');

var filtered2000 = collection.filterDate('2000-01-01', '2000-12-31');

var ndvi = filtered2000.median();

print(ndvi)

var geometry = ee.FeatureCollection('ft:1peQqvrT-WKfsXa_c63gI5h0-Wo_LmBuMnrEg3C').filterMetadata('geometry_vertex_count', 'equals', 10642)

Map.addLayer(ndvi, {palette: '000000, 00FF00', min: 0, max:.8});

//compare the ndvi in 2000 and 2010

var collection = ee.ImageCollection('LANDSAT/LE07/C01/T1_8DAY_NDVI');

//filter for the year 2000

var filtered2000 = collection.filterDate('2000-01-01', '2000-12-31')

//filter for the year 2010

var filtered2010 = collection.filterDate('2010-01-01', '2010-12-31');

//identify the median pixel value per year for 2000 and 2010

var ndvi2000 = filtered2000.median();

var ndvi2010 = filtered2010.median();

//subtract the 2000 ndvi values from the 2010 ndvi values

var difference = ndvi2010.subtract(ndvi2000).clip(geometry);

//add layer with loss in red, gain in green

Map.addLayer(difference,
 {palette: 'FF0000, 000000, 00FF00', min: -0.3, max: 0.3});

//count pixels where loss occurred

var stats = difference.reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: geometry,
  scale: 30,
  maxPixels: 1e9
});

print (stats);
print('Area where vegetation loss ocurred: ', stats.get(''), 'square meters')


texts:
I'm trying to complete the following steps. What tools should I use?

Convert 2km buffer polygons into polylines
Select random point along these lines for each feature using Random Point tool
Generate new buffer of radius 4km around these points


texts: Creating ArcMap toolbar button from model tool?

texts:
Using Arc 10.3.1 ModelBuilder, I've just successfully created my first tool, which is stored in a toolbox within "My Toolboxes".  I run the tool by opening the Arc Toolbox button, dropping down to the tool and double-clicking it.

Now I'd like to make a toolbar button from my tool.  I envision running the tool by simply clicking the new button rather than scrolling down through the Arc Toolbox dropdown.  I've read through the book "Getting to Know ArcGIS Modelbuilder", plus a web search, plus the ESRI help documents.  Nothing I've seen helps...

In summary, how do I convert a tool to a toolbar button?


texts: st_intersects of zero area polygon and point returns false

texts:
I have the following three geometries:

PT: Point(0 0)

LINE: LineString(0 0, 0 10)

POLY: Polygon((0 0, 0 0, 0 0, 0 0))

Calling st_intersects on LINE and POLY returns true (as though POLY describes Point(0 0).  But calling st_intersects on PT and POLY returns false.   I could see POLY being considered an erroneous geometry so false might be the right answer but I would expect st_intersects(PT,POLY) == st_intersects(LINE,POLY)


texts: Incomplete Tile Package from ArcGIS Desktop?

texts:
I am running a tile package, but when it finishes the bottom third of the tile is not visible. The only thing that has changed recently is that I have upgraded from version 10.4 to 10.5. A co-worker with version 10.5 runs the same tile on his machine with the same .mdx that I am using and gets a complete tile.

Any ideas on what may be the issue?


texts: Digitized contours on map and interpolate contour line to create surface

texts:
I have digitized the contours on a map and I want to recreate the surface by interpolating (using IDW) using these contour line. When I generate contours (extract) to see how well they match with the input parameters (digitized contours they are quite off.  I set the interpolation for IDW with a P of 20, if I use 2 then bulls eyes form on the raster.
I am not sure if I am using the right approach. should convert the lines to points? Use a different interpolation method etc?


texts: ArcGIS Pro - Basemap in 3D scene is wrong resolution

texts:
I recently started playing around with ArcGIS Pro again and I've noticed that any basemap used with a 3D scene, in local or global view, is the wrong resolution when I run ArcGIS Pro through remote desktop.

It looks as if the zoom level has been shifted a level or two:

compared to ArcGIS Pro running on my local machine:

Zooming and panning makes no difference, and every preset basemap in ArcGIS Pro has the same behavior.

I can't find any differences between my local machine and the machine that I remote into. They are both running ArcGIS Pro 1.2 alongside ArcGIS 10.2.1.

I'm fairly certain this wasn't occurring with ArcGIS Pro 1.1, but I am unable to locate an installer to verify.


texts: ArcPy findAndReplaceWorkspacePath not working with Raster Datasets

texts:
Situation:

I am working on creating a Python script that will replace workspace paths for every layer in a map document in order to fix broken layers after a data migration.  The script works except for when run on MXDs that have a raster dataset in them.  I have isolated the problem to the code below.  Specifically the .findAndReplaceWorkspacePath() method seems to be failing when attempting to perform the operation on the raster layer.

Code:

fullPathToNewWorkspace = r"C:\GIS\New\Data" # for context, not actually used in script
fullPathToOldWorkspace = r"C:\GIS\Old\Data" # for context, not actually used in script
mxd = arcpy.mapping.MapDocument(r"Path\To\MXD")

for df in arcpy.mapping.ListDataFrames(mxd):
    for lyr in arcpy.mapping.ListLayers(mxd, '', df):
        if lyr.supports("DATASOURCE"):
            print lyr.dataSource # added to help troubleshoot
            lyr.findAndReplaceWorkspacePath(r"Public", r"G-DRIVE")
            print lyr.dataSource # added to help troubleshoot

My Troubleshooting:

The troubleshooting print calls are as such:

> C:\GIS\Old\Data\VectorData1.shp
> C:\GIS\New\Data\VectorData1.shp
> C:\GIS\Old\Data\VectorData2.shp
> C:\GIS\New\Data\VectorData2.shp
> C:\GIS\Old\Data\RasterData # ArcGRID format
> # No print output for C:\GIS\New\Data\RasterData

So the script is getting caught up on trying to .findAndReplaceWorkspacePath for the raster.  Furthermore, when I remove the raster dataset from the mxd and rerun the same script it works fine.

Error Message:

This is the traceback I get from running it with the raster dataset in the MXD:

"C:\Users\Chaz\Desktop\TestDataMig\Scripts\DataMigration_v1.2_RasterTesting.py",     line 1, in <module>
    import arcpy
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\utils.py", line 181, in fn_
    return fn(*args, **kw)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\_mapping.py", line 695, in findAndReplaceWorkspacePath
    return convertArcObjectToPythonObject(self._arc_object.findAndReplaceWorkspacePath(*gp_fixargs((find_workspace_path, replace_workspace_path, validate), True)))
ValueError: Layer: Unexpected error

Question:

I was under the impression (ArcPy documentation) that raster datasets worked fine with the .findAndReplaceWorksapcePath method, yet I am still having these issues.  Am I missing something?  Do I need to phrase my arguments differently, or do I need to use a different method to replace a raster path?

Environment:

Windows 7
ArcGIS/ArcPy 10.2
Python 2.7.3


texts: When is a projected length calculated via a pythagorean distance more accurate than an ellipsoidal length?

texts:
I'm using spatialite 3.0.1 with the following test data (in EPSG: 4326):

Name        astext(geometry)
----------  --------------------
Belitsa     POINT(23.563 41.947)
Malko Tarn  POINT(27.533 41.983)
Kotel       POINT(26.4752 42.938

My client demands I provide the projected st_length distance instead of geodesiclength.

(compare the following queries)

select a.name, b.name, st_length(transform(makeline(a.geometry, b.geometry), 32635)) from ProjCist a, ProjCist b;

Belitsa     Malko Tarn  329056.08890108
...

versus

select a.name, b.name, geodesiclength(makeline(a.geometry, b.geometry)) from ProjCist a, ProjCist b;

Belitsa     Malko Tarn  329092.535472995
...

The 50 meter difference in length, for our purposes, is hugely significant. Which is more correct? Why?


texts: Converting case of all values in all fields of table using ArcPy & Python?

texts:
Does anyone know how to use the .upper() in a stand alone python script so it can be used to change all characters in a table into UPPER case?  I don't want to use the field calculator since I have multiple tables, etc...

I currently am using a cursor to correct mistakes in my tables, but can't figure out the correct syntax to use it to convert all the strings to UPPER case!

The above question would apply equally well when the change of case desired is to lower or title.


texts:
How to use SQL API output to create a layer in a map?

I thought it could be achieved with createLayer passing the data directly, i.e. not with URL as second argument, but the data I am getting from SQL API has not the structure of a layer.

Any advice?


texts: converting shp to fileGDB with ogr2ogr changes geometry type?

texts:
I'm converting a shapefile to file geodatabase with GDAL 1.11.3 and the ESRI FileGDB driver.

ogrinfo on the .shp gives:

$ ogrinfo NED1_Index.shp
INFO: Open of `NED1_Index.shp'
  using driver `ESRI Shapefile' successful.
1: NED1_Index (Polygon)

I run this command:

$ ogr2ogr -f "FileGDB" nedtest.gdb NED1_Index.shp

or this one specifying the geometry type:

$ ogr2ogr -f "FileGDB" -nlt POLYGON nedtest.gdb NED1_Index.shp

and then ogrinfo on the output fileGDB gives:

$ ogrinfo nedtest.gdb/
INFO: Open of `nedtest.gdb/'
  using driver `FileGDB' successful.
1: NED1_Index (Multi Polygon)

Why is the geometry changing to MultiPolygon?


texts:
I was checking the best way to cut linestrings by points.

The scenario is: lots of streets, need the segments cut by intersect points, kind of this:

I got

linestrings (full uncuted by points) table
st_intersection points table

I need to get the independent linestring segments cut by the intersection points table.

I'm using PostGIS functions, and found several approaches, but every one of them has giving me some kind of issue.

Here's what I already tested:

1

Line table: 1 row , st_memunion of 1200 lines
Point table: 1700 rows (points)

What's bad: really takes really lot of time and memory flush. Can't create multiple tables at same time cause memory just can't handle it.
And the result is dirty and messy. Instead of giving me the correct row number, and I need to clean it up later (well explained here Split lines at intersection points)

CREATE TABLE lines_with_messy_result AS (
SELECT
((ST_DUMP(ST_SPLIT(a.geom,b.ix))).geom) as geom
FROM st_union_lines a
INNER JOIN lots_of_points b
ON ST_INTERSECTS(a.geom, b.ix)
);

--then need to clean this up
create table lines_segments_cleaned as (
SELECT DISTINCT ON (ST_AsBinary(geom))
geom
FROM
lines_with_messy_result
);

source of this way/approach:  https://stackoverflow.com/questions/25753348/how-do-i-divide-city-streets-by-intersection-using-postgis

2

Same lines/points table.
Still messy results and need to clean this up.
Still lots of time to finish query.

--messy table
Create table messy_lines_segments as
Select
row_number() over() as id,
(st_dump(st_split(input.geom, blade.ix))).geom as geom
from st_union_lines input
inner join lots_of_points blade on st_intersects(input.geom, blade.ix);

--then cleaning the messy table
delete from messy_lines_segments a
where exists
(
select 1 from messy_lines_segments b where a.id != b.id and st_coveredby(b.geom,a.geom)
);

source of this way / approach : Split lines at intersection points

3

I also found this function:
https://github.com/Remi-C/PPPP_utilities/blob/master/postgis/rc_Split_Line_By_Points.sql

which has the good thing that it does not leave a messy_result that then I need to clean it up. But you do need st_memunion from both sides (lines table and points table)

It's kind of:

create table osm.calles_cortadas_03segmentos_sanluis as (
SELECT result.geom
FROM
osm.calles_cortadas_01uniones_sanluis AS line,
osm.calles_cortadas_00intersecciones_sanluis AS point,
rc_split_line_by_points(
input_line:=line.geom
,input_points:=point.ix
,tolerance:=4
) AS result
);

But its also super-long hours for getting the results. And I also tried with longer tables (10k lines, 14k points) and I just get memory out issues.

I also tried Esri's ArcGIS with also bad results...

So, what's the best way of getting this done with PostGIS geom functions?

I mean, without stepping into topology.

Or what's your best recommendation?


texts: TIN to Raster Tool, output elevations are 3.40282e+038 and -3.40282e+038, tried calculate statistics didn't work

texts:
I have an input TIN file that I created using Create TIN. The elevations are all correct. When I convert the TIN to Raster using the tool, my output elevations are 3.40282e+038 and -3.40282e+038

I've checked projections (and probably do it again) and tried calculating statistics, but so far nothing has changed. Someone else used the same file I was using, output it to the same folder I have been, and they got it to work.

EDIT: I edited out all but the end to file names. My S drive is a server drive.

arcpy.CreateTin_3d("C:/...Default.gdb/water_level","PROJCS['NAD_1983_HARN_Adj_MN_Washington_Feet',GEOGCS['GCS_North_American_1983_HARN',DATUM['D_North_American_1983_HARN',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500019.9183995],PARAMETER['False_Northing',100003.9836799],PARAMETER['Central_Meridian',-92.83333333333333],PARAMETER['Scale_Factor',1.000039836799],PARAMETER['Latitude_Of_Origin',44.74583333333334],UNIT['Foot_US',0.3048006096012192]]","B_XS100 Max_Water Hard_Line <None>;Boundary <None> Hard_Clip <None>","DELAUNAY")

# Replace a layer/table view name with a path to a dataset (which can be a layer file) or create the layer/table view within the script
# The following inputs are layers or table views: "water_level"
arcpy.TinRaster_3d("water_level","S:/.../water_level_3.tif","FLOAT","LINEAR","CELLSIZE 100.485944","1")


texts:
Creating a dynamic buffer with PostGIS the following view is working:

CREATE VIEW buffer500m AS
SELECT gid, st_buffer(geom,500)::geometry(polygon,31492) AS geom
FROM point;

Unfortunately this query gives an error using a SpatiaLite database (unrecognized token: ":").

So I've tried

   CREATE VIEW buffer500m AS
   SELECT gid, st_buffer(geom,500) AS geom
   FROM point;

and

   CREATE VIEW buffer500m AS
   SELECT gid, buffer(geom,500) AS geom
   FROM point;

but with both queries the geometry and the CRS are missing in the view. Can anyone help me?


texts:
I have a point dataset wih 911 calls for one year, including a column with the date (day) when the call was made.

I would like to create a graphic that shows in which month how many calls were made, to analyze the temporal variation. I add a picture to illustrate the way the result should look like.

Can this be done in ArcGIS and if yes how? How do the data have to be prepared, is it necessary to add a column indicating the month and can this be done automated (as there are a lot of features present)?


texts: Changing Colour of Leaflet Layer Control Box

texts: Is it possible to reference the box in the CSS to change the colour of the layer selection box? Would be great to change it away from the standard white.


texts: ReactJS TypeError "not a constructor" on running a component with Carto

texts:
I am trying to integrate a CARTO map according to this example:

Importing Carto:

// Import react-leaflet for the map / basemap components
import { Map, TileLayer as Basemap } from 'react-leaflet';

// Import CARTO.js v4 <3
import carto from 'carto.js';

// Import our custom Layer component (it uses carto.js methods internally)
import Layer from './Layer';

// Voyager basemap <3
const CARTO_BASEMAP = 'https://{s}.basemaps.cartocdn.com/rastertiles/voyager_nolabels/{z}/{x}/{y}.png';

This is the constructor of my React Component:

constructor(props) {
    super(props);

    this.state = {
      center: [40.758313915, -3.67774875],
      zoom: 11,
    },

    // Setup the client in the contructor with our user and apiKey
    this.cartoClient = new carto.Client({ apiKey: 'wadus', username: 'rochoa' });

package.json:

{
  "name": "client",
  "version": "0.1.0",
  "private": true,
  "dependencies": {
    "axios": "^0.18.0",
    "babel-polyfill": "^6.26.0",
    "carto.js": "^4.0.1",
    "classnames": "^2.2.6",
    "jsonwebtoken": "^8.3.0",
    "jwt-decode": "^2.2.0",
    "leaflet": "^1.3.3",
    "react": "^16.4.1",
    "react-dom": "^16.4.1",
    "react-leaflet": "^2.0.0",
    "react-redux": "^5.0.7",
    "react-router-dom": "^4.3.1",
    "react-scripts": "1.1.4",
    "react-simple-storage": "^1.2.1",
    "redux": "^4.0.0",
    "redux-connect": "^7.0.0",
    "redux-thunk": "^2.3.0"
  },
  "scripts": {
    "start": "react-scripts start",
    "build": "react-scripts build",
    "test": "react-scripts test --env=jsdom",
    "eject": "react-scripts eject"
  },
  "proxy": "http://localhost:5000"
}

When running, receiving:

TypeError: __WEBPACK_IMPORTED_MODULE_2_carto_js___default.a.Client is not a constructor

Having a hard time to understand what is wrong here...


texts: Export ReduceRegion fixedHistogram toDrive

texts:
I'm trying to export a fixed histogram from reduceRegion as a .csv file. I'm sure im missing the obvious here, but can't find a reasonable example that does this straightforward. I can do it as a FeatureCollection, but the format is difficult. I was thinking 2 columns for values and counts, like when you export your histograms as from the console.

https://code.earthengine.google.com/296846e47fe937897dda456d4552a39b

var landsat = ee.Image(ee.ImageCollection("LANDSAT/LT05/C01/T1_TOA").filterBounds(geometry).first())

var index = 'B1'

function country (image) {var c = image.reduceRegion({
  geometry: image.geometry(),
  reducer: ee.Reducer.fixedHistogram(0, 1, 1000),
  scale: 30,
  maxPixels: 1e12,
}); return ee.Feature(null, c.select(['B1']))}

var hist = landsat.select(index).reduceRegion({
  reducer: ee.Reducer.fixedHistogram(0, 1, 1000),
  geometry: landsat.geometry(),
  scale:100, //.add(upScale),
  bestEffort: true,
  maxPixels: 1e13
});
print(hist)

var featureTest = ee.Feature(null, hist);
print(featureTest)

// Gives a difficult format
Export.table.toDrive({
  collection: ee.FeatureCollection([featureTest]),
  description: 'exportTableExample',
  fileFormat: 'CSV'
});

var list = ee.Array(hist.values().get(0))
print(list, 'list')

var means = list.slice(1, 0, 1).multiply(1000).round().divide(1000).project([0]).toList()
var counts = list.slice(1, 1, 2).project([0]).toList();
print(means, counts)

var fromList = ee.Dictionary.fromLists(means, counts)
print(fromList)

Thanks!


texts: Calculate length of selected features in PyQGIS

texts:
I create a selection with getFeatures(QgsFeatureRequest()) and I want to calculate the length of the result which is a QgsFeatureIterator :

select = QgsExpression(" myexpression ")
selectedfeats = layer.getFeatures(QgsFeatureRequest(select ))

total = 0
for feature in selectedfeats :
    geom = feature.geometry()
    total += geom.length()

But it doesn't work it returns 0 whereas the selectedfeats is not empty.


texts: Getting Raster exported from ArcMap into Google Earth (.kmz) with transparency preserved?

texts:
RasterLayerExportClass will only create a .png with a single value transparency, however there is a bug in Google earth which prevents this from working.

The other export options don't support transparency at all? (jpg, tiff)

I do notice that the 3D analyst toolbox performs this correctly with the "Layer to KML" function, this I have no source for.


texts: How to extract a buffer layer from a vast amount of individual buffers around polygons using ArcMap?

texts:
I have a vector layer containing nearly 8000 singlepart polygons. It looks like this:

I want statistics on the area that is further than 2,5m away from any border. I thought I could do this by creating a buffer of 2,5m around all polygons, then extracting the buffer area from the original dataset. I manage to create the buffer, but it is not a separate layer, and the polygons that are created are all multi-part with double lines.

When I select a polygon, I see the lines are overlapping, and I need the buffer lines to be separate polygons. I have tried using the multipart to singlepart tool, but it does not solve the problem. Neither does creating a union between the original layer and the buffer layer.

This is one example - here I have selected one polygon in the buffer layer(the original layer behind with grey lines, the buffer layer with see through background ad red lines, selected polygon marked with blue lines)::


texts: netCDF import problem in QGIS (XY issue)

texts:
I working with a netCDF file related to precipitation over time in QGIS. When I import the file, QGIS reads the data indicating y-coordinate, instead of timestep, as the dynamic dimension to associate raster values with.

When I use netCDF browser to select slices of the data, what I am getting is bands plotting x-coordinate and time, separated for each value of the y-coordinate.

What I would like to get is bands plotting precipitation for each x and y coordinate separated for each value of time.

Is there a way to change which dimensions QGIS uses for spatial reference for a netCDF file, either during importation or afterwards?


texts: Sensor Orientation (XYZF) to Lat/Long? triangulation required?

texts:
I have IAGA2002 data regarding Magnetic North, as read from Sensor Orientation (XYZF), documented here: https://www.ngdc.noaa.gov/IAGA/vdat/IAGA2002/iaga2002format.html

Example data is:

 Format                 IAGA-2002                                    |
 Source of Data         IPGP                                         |
 Station Name           ADDIS ABABA                                  |
 IAGA CODE              AAE                                          |
 Geodetic Latitude      09.033                                       |
 Geodetic Longitude     038.767                                      |
 Elevation              2442                                         |
 Reported               XYZF                                         |
 Sensor Orientation     HDZF                                         |
 Digital Sampling       1 hour                                       |
 Data Interval Type     1-hour (00:30 - 01:29)                       |
 Data Type              definitive                                   |
DATE       TIME         DOY     AAEX      AAEY      AAEZ      AAEF   |
2004-01-01 00:00:00.000 001     36034.50    859.60   1318.10  36068.40
2004-01-01 01:00:00.000 001     36037.10    861.60   1318.80  36071.00
2004-01-01 02:00:00.000 001     36042.20    861.70   1319.10  36076.10
2004-01-01 03:00:00.000 001     36035.20    864.40   1321.10  36069.20
2004-01-01 04:00:00.000 001     36018.20    859.30   1322.60  36052.20
2004-01-01 05:00:00.000 001     36017.70    867.10   1330.10  36052.10
2004-01-01 06:00:00.000 001     36041.60    880.30   1332.20  36076.40

Given the (station-name)X, Y, Z, F parameters, and perhaps the source location (Geodetic Lat/Long listed in the headers), how might I derive a specific Lat/Long.

And, will I need to triangulate with multiple stations (which I have many of)?


texts: Transfer data with 'geometry' to Elasticsearch

texts:
I would like to transfer my data which contains 'geometry' to Elasticsearch but I don't know how to change my 'body'. 'geometry' is the last column. What should I change?
Code:
es = Elasticsearch()
bulk = Bulk(es, index_name ='gis-test' , max_size = 5)

budynki = geopandas.read_file(r'/home/gis/budynkipolska/bdot10k_bubd_a.gpkg', rows=slice(0, 50),  encoding = "cp1250")
adresy = geopandas.read_file(r'/home/gis/adresypolska/PRG_PunktyAdresowe_POLSKA.shp', rows=slice(0, 50),  encoding = "cp1250")
budynki = pd.DataFrame(budynki)

def export_elastic(powiat,lokalnyid,wersjaid,status_bdot,nazwa,kategoria_bdot,funkcja_ogolna_budynku,funkcja_szczegolowa_budynku,liczba_kondygnacji,zabytek,skrot_karto,aktualnosc_geometrii,aktualnosc_atrybutow,kod_kst,geometry):
    body2 = {
        'classification': {
                'powiat': powiat,
                'lokalnyid': lokalnyid,
                'wersjaid': wersjaid,
                'status_bdot': status_bdot,
                'nazwa': nazwa,
                'kategoria_bdot': kategoria_bdot,
                'funkcja_ogolna_budynku': funkcja_ogolna_budynku,
                'funkcja_szczegolowa_budynku': funkcja_szczegolowa_budynku,
                'liczba_kondygnacji': liczba_kondygnacji,
                'zabytek': zabytek,
                'skrot_karto': skrot_karto,
                'aktualnosc_geometrii': aktualnosc_geometrii,
                'aktualnosc_atrybutow': aktualnosc_atrybutow,
                'kod_kst': kod_kst,
                'geometry': geometry, #<----- HERE
            }
        }
    return(body2)

for index, row in budynki.iterrows():
    print(row['geometry'])
    powiat = row['powiat']
    lokalnyid = row['lokalnyid']
    wersjaid = row['wersjaid']
    status_bdot = row['status_bdot']
    nazwa = row['nazwa']
    kategoria_bdot = row['kategoria_bdot']
    funkcja_ogolna_budynku = row['funkcja_ogolna_budynku']
    funkcja_szczegolowa_budynku = row['funkcja_szczegolowa_budynku']
    liczba_kondygnacji = row['liczba_kondygnacji']
    zabytek = row['zabytek']
    skrot_karto = row['skrot_karto']
    aktualnosc_geometrii = row['aktualnosc_geometrii']
    aktualnosc_atrybutow = row['aktualnosc_atrybutow']
    kod_kst = row['kod_kst']
    geometry = row['geometry']
    body2 = export_elastic(powiat,lokalnyid,wersjaid,status_bdot,nazwa,kategoria_bdot,funkcja_ogolna_budynku,funkcja_szczegolowa_budynku,liczba_kondygnacji,zabytek,skrot_karto,aktualnosc_geometrii,aktualnosc_atrybutow,kod_kst,geometry)
    res = es.index(index = 'gis-test', body = body2, doc_type = '_doc')
    bulk.append_body(body = body2, id_doc = index, method='update')
bulk.commit()

Screen:


texts: How to color GeoTIFF if its values are from -1 to 1?

texts:
There is a GeoTIFF, 1 channel. How do I color it? Namely, values that were from -1 to 0.02 color in red, from 0.02 to 0.5 color in yellow, and values that are greater than 0.5 color in green.

GeoTIFF shows the calculated NDVI by fields.
I read about GDAL COLOR TABLE, but you can only put integer values there

How can I then apply a vector to this colored GEOTIFF?

Here's how it should end up
I need to do this in Python.
UPD
I uploaded a pseudo color:
INTERPOLATION:DISCRETE
0.02,252,0,4,255,<=0.02
0.5,200,180,0,255,0.02 - 0.05
1,2,177,26,255,0.5 - 1

How to add to GeoTIFF from Python?


texts:
Any advice on which map projection to use for a map of Russia?

I am trying to project my shapefile (WGS84) but not sure which one to select. I would like to maintain equal areas.


texts: How to import country administrative boundaries from OSM Planet to PostGIS Polygons

texts:
I need to import all world countries as polygons or multipolygons into PostGIS. Specifically, I need the administrative boundaries that are like this orange bubble around the actual land borders.

I've already gone through all the motions that I thought would get me the result I wanted: a single polygon or multipolygon for each country. But once I get the data into PostGIS, the polygons just aren't there in the right places, and there are many more rows than I expected.

Here's how I download the planet-latest.pbf file, convert to o5m format, run osmfilter to get data with type=boundary and boundary=administrative and admin_level=2.

$ wget http://planet.openstreetmap.org/pbf/planet-latest.osm.pbf
$ osmconvert planet-latest.osm.pbf -o=planet-latest.o5m
$ time ./osmfilter planet-latest.o5m --keep="type=boundary and boundary=administrative and admin_level=2" -v -o=country_boundaries_3.o5m
    (ran 22m 14.162s)
$ time osmconvert country_boundaries_3.o5m --out-pbf >country_boundaries_3.pbf
    (ran 0m 26.276s)

Then I import into PostGIS:

# create database geocell_countries_3 with owner=geocell tablespace=postgres_huge;
# \c geocell_countries_3
# create extension postgis;
# create extension hstore;
# \q
$ time osm2pgsql --create --verbose --proj 4326 --hstore-all --database geocell_countries_3 /mnt/planet-latest/country_boundaries_3.pbf
    (ran: 2m 27.034s)

I copy the data into a new table:

$ psql geocell_countries_3
# create table country_boundaries (like planet_osm_polygon);
# insert into country_boundaries select * from planet_osm_polygon where admin_level = '2';
# alter table country_boundaries set schema g;
# select AddGeometryColumn ('g', 'country_boundaries', 'multipoly', 0, 'MULTIPOLYGON', 2);
# update g.country_boundaries set multipoly = ST_Multi(way);
# \q

Now there are 48 records just for USA but I just want one.

# select count(*) from g.country_boundaries where (tags -> 'ISO3166-1:alpha3') = 'USA';
 count
-------
    48
(1 row)


texts: QGIS Atlas I want a second map as Overview but not have it on every other page when i print it to pdf

texts:
I'm working with QGIS Atlas and I'm showing in my "map1" a polygon layer that jumps through 24 different polygons.
Now I want to add a second or third map (maybe more depending on the polygon layer) that shows an overview of the polygon layer.
My problem is that when I print or export to PDF it prints the second map on every other page.
I was wondering if there was another option than doing a second atlas and combine the two PDFs.


texts:
Can anyone advise how do I do a hyperlink that directly link to a bookmarked area of a pmf or mxd? Something like what we do to hyperlink to a specific page in pdf.

For example, I have a mxd call BD_WORLD with Regional_View as one of the bookmark.

SO something like this will be my hyperlink:
\\serverABC\GIS\WORLD_MXD\BD_WORLD.mxd#Regional _View

Preferably, I would like to do it in attribute table, if possible.

PS: using 9.3.1, with only Arcview, Spatial Analyst and Publisher.


texts: Switch QGIS proxy settings programatically

texts:
I am using QGIS at two offices with different proxy-servers.

At the moment I have to change the proxy settings manually every time i switch the office.

I was wondering if its possible to change the proxy-settings programatically with PyQGIS? then I could write a plugin to switch between the proxy-settings.

EDIT1:

In the meantime i found a way to change the proxy settings of QGIS but still it's not working.

With this code I can change the settings:

from PyQt4.QtCore import QUrl, QSettings
from PyQt4.QtNetwork import QNetworkRequest, QNetworkProxy
from qgis.core import QgsNetworkAccessManager

my_settings={"Proxy enabled": u'proxy/proxyEnabled', "Proxy Host ": u'proxy/proxyHost', "Proxy Port": u'proxy/proxyPort'}
fiddler={"Proxy enabled": True, "Proxy Host ": "localhost", "Proxy Port": 8888}
freiburg={"Proxy enabled": True, "Proxy Host ": "192.168.95.165", "Proxy Port": 8080}
aus={"Proxy enabled": False, "Proxy Host ": "192.168.95.165", "Proxy Port": 8080}

current_choice=aus

s = QSettings() #getting proxy from qgis options settings

for key, val in my_settings.iteritems():
    #print str(key)+":"+str(val)
    settings_key=key
    #print str(settings_key)
    # Get user defined current setting
    for key2, val2 in current_choice.iteritems():
        if key2==settings_key:
            #print key
            #print val
            settings_val=val2
    current_setting = s.value(str(val).decode('unicode-escape'))
    #print str(val).decode('unicode-escape')
    #print str(key)+": "+str(current_setting)
    s.setValue(unicode(str(val)), settings_val)
s.sync()

# procedure to set proxy if needed

proxyEnabled = s.value("proxy/proxyEnabled", "")
proxyType = s.value("proxy/proxyType", "" )
proxyHost = s.value("proxy/proxyHost", "" )
proxyPort = s.value("proxy/proxyPort", "" )
proxyUser = s.value("proxy/proxyUser", "" )
proxyPassword = s.value("proxy/proxyPassword", "" )
if proxyEnabled == "true": # test if there are proxy settings
   proxy = QNetworkProxy()
   if proxyType == "DefaultProxy":
       proxy.setType(QNetworkProxy.DefaultProxy)
   elif proxyType == "Socks5Proxy":
       proxy.setType(QNetworkProxy.Socks5Proxy)
   elif proxyType == "HttpProxy":
       proxy.setType(QNetworkProxy.HttpProxy)
   elif proxyType == "HttpCachingProxy":
       proxy.setType(QNetworkProxy.HttpCachingProxy)
   elif proxyType == "FtpCachingProxy":
       proxy.setType(QNetworkProxy.FtpCachingProxy)
   proxy.setHostName(proxyHost)
   proxy.setPort(int(proxyPort))
   proxy.setUser(proxyUser)
   proxy.setPassword(proxyPassword)
   QNetworkProxy.setApplicationProxy(proxy)

This works so far as I can see the changed settings in the QGIS UI (settings->options).

The settings are also written to the windows registry but the changes won't have any effect until i click the OK button in the QGIS settings dialog.

You can test this by setting the proxy programmatically to some proxy-settings that should prevent QGIS from accessing the internet (e.g. localhost:98765) and try to load and pan through a wms-layer.
Any idea what's missing?

Edit2:
I just piped the output from qgis to a file and had a look at what is going on when I change the proxy-settings using the GUI:

src/core/qgsnetworkaccessmanager.cpp: 364: (setupDefaultProxyAndCache) [9134ms] setting proxy 3 192.168.95.165:8080 /
src/core/qgsnetworkaccessmanager.cpp: 167: (setFallbackProxyAndExcludes) [0ms] proxy settings: (type:HttpProxy host: 192.168.95.165:8080, user:, password:not set

So there are two functions called (setupDefaultProxyAndCache and setFallbackProxyAndExcludes). Perhaps something like that has to be done when using pyQGIS to change the settings?


texts: leaflet - Drawing line following the road


texts:
I just now downloaded a plugin named "XY tools"from fetch plugins option . I am a user of Qgis v 1.9.0.It shows me the following error:
An error has occured while executing Python code:

Traceback (most recent call last):
  File "C:\Documents and Settings\15491/.qgis2/python/plugins\xytools\xytools.py", line 123, in xyToolClick
    if self.layerInfo.has_key(self.layer) or self.getXyColumns(self.layer):
  File "C:\Documents and Settings\15491/.qgis2/python/plugins\xytools\xytools.py", line 296, in getXyColumns
    dlg = XyToolsDialog(self.XY_COMBO_MSG, layer.dataProvider().fields(), layer.name())
  File "C:\Documents and Settings\15491/.qgis2/python/plugins\xytools\xytoolsdialog.py", line 35, in __init__
    for field in fieldMap.values():
AttributeError: 'QgsFields' object has no attribute 'values'

Python version:
2.7.2 (default, Jun 12 2011, 15:08:59) [MSC v.1500 32 bit (Intel)]

QGIS version:
1.9.0-Master Master, 1a723b4

Python path: ['C:\\Documents and Settings\\15491\\.qgis2\\python\\plugins\\GeoCoding', 'C:\\OSGeo4W\\apps\\qgis-dev\\python\\plugins\\sextante', 'C:\\Documents and Settings\\15491\\.qgis2\\python\\plugins\\GeoCoding', 'C:\\Documents and Settings\\15491/.qgis2/python/plugins\\sextante_taudem', 'C:\\Documents and Settings\\15491/.qgis2/python/plugins\\sextante_animove', 'C:\\Documents and Settings\\15491/.qgis2/python/plugins\\sextantelwgeomprovider', 'C:\\Documents and Settings\\15491/.qgis2/python/plugins\\sextanteexampleprovider', 'C:\\Documents and Settings\\15491/.qgis2/python/plugins\\permaclim', 'C:\\Documents and Settings\\15491\\.qgis2\\python\\plugins\\GeoCoding', 'C:/OSGeo4W/apps/qgis-dev/./python', 'C:\\Documents and Settings\\15491/.qgis2/python', 'C:\\Documents and Settings\\15491/.qgis2/python/plugins', 'C:/OSGeo4W/apps/qgis-dev/./python/plugins', 'C:\\OSGeo4W\\bin\\python27.zip', 'C:\\OSGeo4W\\apps\\Python27\\DLLs', 'C:\\OSGeo4W\\apps\\Python27\\lib', 'C:\\OSGeo4W\\apps\\Python27\\lib\\plat-win', 'C:\\OSGeo4W\\apps\\Python27\\lib\\lib-tk', 'C:\\OSGeo4W\\bin', 'C:\\OSGeo4W\\apps\\Python27', 'C:\\OSGeo4W\\apps\\Python27\\lib\\site-packages', 'C:\\OSGeo4W\\apps\\Python27\\lib\\site-packages\\PIL', 'C:\\OSGeo4W\\apps\\Python27\\lib\\site-packages\\win32', 'C:\\OSGeo4W\\apps\\Python27\\lib\\site-packages\\win32\\lib', 'C:\\OSGeo4W\\apps\\Python27\\lib\\site-packages\\Pythonwin', 'C:\\OSGeo4W\\apps\\Python27\\lib\\site-packages\\wx-2.8-msw-unicode', 'C:\\Documents and Settings\\15491\\.qgis2\\python\\plugins\\mmqgis/forms', 'C:\\Documents and Settings\\15491\\.qgis2\\python\\plugins\\pdokbaggeocoder/forms', 'C:\\OSGeo4W\\apps\\qgis-dev\\python\\plugins\\fTools\\tools']

what is the problem?


texts:
is there a way in Arcpy (arcgis 10.2.1) to buffer lines by selecting thier objectid?

arcpy.Buffer_analysis("MY_JUNK","d:/data/ArcGIS/Default.gdb/MY_JUNK_Buffer1","10 Meters","FULL","ROUND","NONE","#")

Currently the copout is to manually select the objects but they are changing over time and it is easier to select by objectid if possible.


texts: SAGA Slope Stability algorithms missing in QGIS 2.18.11

texts:
I have install QGIS 2.18.11 via OSGEO4W 64 bit (that programs and versions install with QGIS SAGA version 2.3.2 GRASS version 7.2.1).
I can use SAGA algorithms without problem in QGIS using processing toolbar but not all algorithms from this SAGA version.
If I use SAGA stand-alone, I can use the terrain analysis Slope Stability toolbar. That toolbar has six SAGA algorithms.

Why are algorithms from this toolbar missing from the Processing toolbox?

Is it possible to add these algorithms manually?


texts: QGIS 3 shows and saves raw default-values from spatialite-tables

texts:
QGIS shows default values from sqlite tables as string in the edit-form when i try to create a new row. After save, the raw default values are stored in the table. How can i edit in QGIS LTR 3.4.7 a sqlite-database without messing up the default values from the tables?

QGIS-form with new feature and 'wrong' default-values. Feature from sql-insert with id 1 in background.

QGIS-table after save with pointless default-values.

tried "solutions"

with "sql-insert" everything works as exepected
same behaviour on different workstations
same behaviour for sqlite- and spatialite-databases
same behaviour with or whitout geometry
same behaviour in QGIS 3.4.7 (LTR) and 3.10
same behaviour in Windows and Linux
same behaviour with default-values in a geopackage
in spatialite InvalidateLayerStatistics(); and UpdateLayerStatistics(); does not help.
Project Properties -> Data Source -> Check "Evaluate default values
on provider side".
(maybe) relevant google-search results :

https://qgis.org/en/site/forusers/visualchangelog30/index.html#feature-detect-literal-default-values-for-spatialite-provider
https://github.com/qgis/QGIS/commit/0ae610c5e3dd364385c2fd91f901e706a04101f1

Create statement for the table in the sqlite-database

CREATE TABLE my_table (
    id integer primary key autoincrement,
    comment text,
    created_at_01 text DEFAULT (datetime('now','localtime')),
    created_at_02 text DEFAULT CURRENT_TIMESTAMP
);

INSERT INTO my_table (comment) VALUES ('hello world');

SELECT * FROM my_table;

id          comment      created_at_01        created_at_02
----------  -----------  -------------------  -------------------
1           hello world  2019-12-09 18:00:12  2019-12-09 17:00:12


texts: Using Python to extract flood hazard info from National Flood Hazard Layer (NFHL) using any given lat/long point?

texts:
Can Python be used to extract flood hazard classifications from National Flood Hazard Layer (NFHL) site using any given lat/long point?

Specifically the following seven classifications:

1% Annual Chance Flood Hazard
Regulatory Floodway
Special Floodway
Area of Undetermined Flood Hazard
0.2% Annual Chance Flood Hazard
Future Conditions 1% Annual Chance Flood Hazard
Area with Reduced Risk Due to Levee

Flood Hazard Zones (28) appears to contain these seven classifications.

Additional information:

NFHL GIS Services user guide

Definitions of FEMA Flood Zone Designations

some notes:

nfhl_layers = fiona.open(county_nfhl_file_path) nfhl_layer_1 = shape(nfhl_layers['geometry'] lat = 47.661694 long= -122.192574

-- Given the lat/long, calculate the seven flood hazards for that point

This function returns true if that point is in the layer and false if it is not.

def is_lat_long_in_nfhl_layer(lat, long, nfhl_layer): pass print is_lat_long_in_nfhl_layer(lat,long, nfhl_layer_1) def get_nfhl_predicates(lat,long): ...... return {"1%_annual_chance" : percent_chance_1, "0.2%_annual_chance": percent_chance_0_point_2,}


texts: How to display a "Save As" dialog for raster files?

texts:
I'm writing a QGIS plugin in python. I need to show a "Save as" dialog window in order to allow the user to select the path of the output raster file. Do I need to create a QFileDialog and manually add all the raster formats supported by QGIS or does it exist a method to load these formats from somewhere? Or, better, is there a method to show a "Save As" dialog window ready to manage raster files?

Thanks!


texts: How to transpose reclassified field values from one polygon shapefile to identical polygons of another shapefile?

texts:
I used vector output from modeling software as clip features for ecoregion input.  I reclassified and dissolved the clipped layers into 3 subunits, and now would like to transpose those simplified values to the source ecoregion layer in areas of intersection, such that the same 3 field values from the clipped layers are mapped onto the source ecoregion layer.

I need to do this for about 200 pairs of shapefiles, and it seems like there must be some at least partially automated way to approach this.

In the image below, I've symbolized the 3 field values semi-transparently over the source ecoregion layer.  You can see that a single value polygon in the top layer overlaps multiple polygons in the source layer, but the outer boundaries of the individual polygons of the source layer exactly match the simplified polygons of the top layer.

In this image, I've reclassified the source ecoregion layer to match the field values of the top layer, such that the field values from any polygons from the top layer that intersect with polygons in the source layer are copied to the attribute table of the source layer.

I've been doing this via manual selection and the attribute table field calculator, but this will take forever for 200 sets of vector data.

Anyone know of either a canned tool in arcGIS that could help with this, or python coding approach that could help speed this up?

Thanks in advance for your help!

EDIT in hopes of clarifying:
For polygons in the reclassified Layer A and Layer B (ecoregions) that intersect, I'd like to populate the subunit field in Layer B with the same value as the subunit field in layer A.

This is the attribute table of the reclassified Layer A:

This is the attribute table of the source ecoregion Layer B after I've manually populated the subunit fields based on intersection with Layer A:


texts: How to create points at shortest euclidean distance to other points?

texts:
So, I have a lot of points, and I need to site three new points so they will be at the shortest euclidean distance between all of the other points.

In other words, these three new points that I need to site will minimize the sum distance that the large group of points will be from any new point.  I need to find the most efficient positioning the new points so the group of points will have the smallest total distance from a new point.

I'd like to be able to do this in Model Builder, or at least develop some kind of process that I can translate into Model Builder and possibly a python script later on (though I'm really bad at python).  I assume this will involve some kind of iterative procedure, I just don't really know where to start.

Anyone have any awesome ideas?

P.S. This is in ArcGIS desktop 10.0

Edit:  The idea is to site three new points that would be distributed to minimize total sum difference, thus they would, I guess, be taking a subset of the total point set to find the median center of them. It would look something like this (where the orange points would be something like the theoretical new locations):

Edit 2: Okay, this has been confusing, but I was going on the description I was given for this particular problem.  Yes, the minimized sum distance is the centroid.  Mathematically, this is the best solution I've come up with so far, where the goal is the minimize Dist, but it just gets the total centroid:

I've now talked to a spatial statistics instructor at an Ivy league university and two Ph.D.'s in math and we still can't figure out an appropriate answer that satisfies what I've been asked.

The only other thing that I can think to do is to somehow evenly divide the original points into three spatially equal parts and find their centroids, but I don't know how to do that.


texts: Why does QGIS provide [N / A] next to the categories?

texts:
What are the reasons why when I ask to show the total number of elements with a categorized stylization it shows me [N / A] next to each category?

The style is well constructed and the total number of elements for the layer works. I insert an image


texts: Creating polygon with attributes using ArcGIS Desktop?

texts:
I'm new of the use of arcGIS.
I have xy data that descibe a polygon. But when I use the tool "point to line" I lost all attributes. There is a way to create a polygon with attributes. I try to explain with an example

Example
x (25 125 55 36),
y (25 336 12 12),
turtle: 3

x and y are the cordinates that describe a sea portion sea. Inside this part there were three sightings of turtles. How can I create a polygon that with the number of sightings in other to create a map of density?

I'm using ArcGIS 10.


texts: How to be notified of a change that is made on an alternative in another session in GE Smallworld?

texts:
I have two Smallworld sessions (Core 4.x) both are on same alternative, the first one has made the alternative writable and the other one is obviously in readonly mode.

I want my readonly session to be notified and refresh the map each time the writable session commits a change.

If I am using a short transaction, I need to start and end a dummy short transaction on readonly session, and if there I am in a long transaction mode I need to switch(:write)or rollforward() so it updates the database and map respectively on my readonly session. But in both cases the readonly session is not being notified of changes and it has to check for changes regularly through a light weight background thread. I have already tried observer and added one to gis_ds_view, yet it does not work when changes that are made in another session.

Is there any way that I could be aware of changes in another session?


texts: How do I get that polygon from Javascript ArcGIS API?

texts: I would like to have a user draw a Polygon (or Line) onto a web map with the Javascript ArcGIS API drawing tools. For the life of me I can't figure out how to pull that polygon from the client and get it to my server. It baffles me. Help?


texts: CREATE view in Geoserver based on current timestamp

texts:
I am trying to use the "Create NEW SQL view" function in Geoserver to create a new view, in which I wish to show a subset of records from a table in MSSQL based on a value in another table.
I have two layers in the same MSSQL DATABASE one with some geometries ("LER_REQUEST_MAP") and another with some needed information ("LER_REQUEST_V2"). One of the informations in the later is an ESRI based timestamp field called "DIG_PERIOD_FROM".

I wish to use my view to:

a) join these two tables (INNER JOIN?), so that I get the information on DIG_Period_FROM to my geometry layer (LER_REQUEST_MAP) without compromising the original table.

b) SELECT * from "the new joined table" where the field DIG_PERIOD_FROM > CURRENT_DATE

Is this possible?


texts:
I use Quantum GIS 2.2.0 Valmiera

I have a polyline layer. I can drag&drop polyline nodes by using "Node Tool". Is it possible to drag&drop line precisely 6 meters to right or left?


texts: 3D variography: rotation matrices for multiple angles support

texts:
I'm currently working on a 3D variogram algorithm based on geostatsmodels (which only handles 2D datasets). I'm having a hard time understanding how to have rotation along the z-axis for azimuth direction handling to work.

From the sources available online that implement 2D directional variograms (scikit-gstat, pyGeoStatistics), I understand that I need to apply rotation matrices to the found paired vectors.

Currently, I'm computing the initial vector (difference of two points in space) then applying the 3D rotation matrices (Rotation Matrix on Wikipedia), then computing the new angles from the rotated coordinates.

But I'm still facing differences in my variogram compared to what I compute in R, when changing the azimuth angle along the z-axis (for exemple east direction, ie 90 degrees). Probably because the new angles are incorrect.

For the moment I only use the "yaw" matrix for horizontal rotation:

I wanted to know if this approach seems correct, and if yes, if I was rightly considering the good matrix ?


texts: How to calculate specific rows in layer field- ArcGIS 10.3

texts:
I try to calculate specific rows (OBJECTID 23-29) in "field_name" field with field calculator- i try to fill in the value "rt8". This is the attribute table:

but i can't find any option to calculate specific rows in it- i can calculate for the whole filed-and that isn't what i'm seeking. I must mention that i want to calculate without selecting features.