texts: Storing two kinds of point geometry in shapefile or creating multipoint shapefiles?

texts:
Following How many interior and exterior rings can a polygon have in a standard ArcGIS shapefile?, can I have a mixture of two kinds of geometry in a shapefile of type point. Just like I have about shapefiles of type polygon or I should always create a shapefile of type multipoint in order to support multipoints in my layer?

My problem is as you see at How can I use ogrinfo to reach information about a .shp?, I can't use ogrinfo in order to understand that the geometry type of a multipoint shapefile will be known as wkbMultipoint or wkbPointwhen I use this line of code:

OGRwkbGeometryType GeometryType = poLayer ->GetGeomType();

for a layer of type multipoint before programming.

Will OGR know a layer of type multipoint as wkbpoint or
wkbmultipoint?
When we create a shapefile in ArcCatalog, we have point and
multipoint, but not multipolyline or multipolygon for polyline and
polygon.Why?


texts:
I have read that Geoserver fully supports WFS-Transactions. But for example here somebody also mentioned that shapefile-datastores would not be the best choice for using Transactions in Geoserver, and that it would be better to do the transactions in a database-store, such as PostGIS.
Things may have changed over the years, so my question is:

Does Geoserver has any limitations in using WFS-T with Shapefiles?


texts:
I'm new to R and to lidR.
I have LAS files that are already ground classified and I'm attempting to use the method below to classify the vegetation and then write a new LAS file retaining the ground and vegetation classifications. The resulting LAS files end up with every point classified as vegetation.
What am I doing wrong?
Here's my script:
library(lidR)
las <- readLAS("D:/tree-segmentation/test.las", select = "xyzc", filter = "-drop_z_below 0")
las <- segment_trees(las, li2012(R = 3, speed_up = 5))

plot(las, color = "treeID")

vegetation_class = function(Z)
{
  mz = max(Z)
  if (mz < 10)
    return(LASLOWVEGETATION)
  else if (mz >= 10 & mz < 22)
    return(LASMEDIUMVEGETATION)
  else
    return(LASHIGHVEGETATION)
}

las@data[!is.na(treeID), Classification := vegetation_class(Z), by = treeID]
plot(las, color = "Classification")
writeLAS(las, "D:/tree-segmentation/test-veg-classified.las", index = FALSE)


texts: Semi-circle representation in legend - ArcGIS!

texts:
I have a little problem that needs solving.
How to represent the legend in ArcMap as in the picture above?

Or as semi-circle instead of complete circles?


texts: Extract By Mask gives different result (and WARNING 000622) in ArcPy than from tool dialog?

texts:
I'm using ArcMap 10.3.1

When I use the gp tool to extract by mask, I get the correct result with any area outside the mask set as NoData.

When I do the same thing in a script, I end up with cells with other values outside the mask.

Note the values within the mask are correct, the stretching is different because the black cells outside the mask have a value of -32768 (the minimum value in a 16 bit signed raster).

It looks like something to do with tiling?

Note when using the gp tool, I get the message "WARNING 000632: Datum conflict between map and output." I don't know if this is relevant?

Below is my script.

Any ideas?

import arcpy
from arcpy import sa as sa
ws = 'c:/temp/test/'
arcpy.env.workspace = ws
outMask = sa.ExtractByMask('band2.tif', 'Clip109071.shp')
outMask.save('band2clipped.tif')


texts: Intersecting a raster with a polygon using PostGIS - artefact error

texts:
I'm using PostGIS2.0 to do some raster/polygon intersections. I'm having difficulty understanding which operation I should use, and what the quickest way of performing this is. My problem is as follows:

I have a polygon and a raster
I want to find all the pixels that fall within the polygon, and get the sum of the pixel value
And (updated problem): I'm getting massive values for some pixels that do not exist in the original raster when I perform the query

I'm having difficulty understanding whether I should use ST_Intersects() or ST_Intersection(). I also don't know what the best approach for summing my pixels is. Here is the first approach that I have tried (#1):

SELECT
    r.rast
FROM
    raster as r,
    polygon as p
WHERE
    ST_Intersects(r.rast, p.geom)

This returns a list of rast values, which I'm not sure what to do with. I tried calculating the summary statistics using ST_SummaryStats() but I'm not certain if this is the weighted-sum of all the pixels that lie within the polygon.

SELECT
        (result).count,
        (result).sum
FROM (
         SELECT
            ST_SummaryStats(r.rast) As result
         FROM
            raster As r,
            polygon As p
         WHERE
            ST_Intersects(r.rast, p.geom)
    ) As tmp

The other approach which I have tried (#2) uses ST_Intersection():

 SELECT
        (gv).geom,
        (gv).val
 FROM
 (
    SELECT
        ST_Intersection(r.rast, p.geom) AS gv
    FROM
        raster as r,
        polygon as p
    WHERE
        ST_Intersects(r.rast, p.geom)

      ) as foo;

This returns a list of geometries that I analyze further, but I assume this is less efficient.

I'm unclear on which is the quickest order of operation also. Should I always choose raster, polygon or polygon, raster, or convert the polygon into a raster so that it is raster, raster ?

EDIT: I updated approach #2 with some details from R.K.'s link.

Using approach #2, I have noticed the following error in the results which is part of the reason why I did not understand the output. Here is the image of my original raster, and an outline of the polygon that is being used to intersect it, overlaid on top:

And here is the result of the intersection using PostGIS:

The problem with the result is that there are values of 21474836 being returned, which are not in the original raster. I've no idea why this is occurring. I suspect it is related to small numbers somewhere (dividing by almost 0), but it returns the wrong result.


texts: Convertion of a .dem file to coordinates and elevation

texts:
I have a .dem-file with elevation data of an area and I am looking to programmatically convert this data to coordinates and elevation in C++. After a lot of searching, it looks like GDALTranslate might do the trick.

Is .xyz a good format when coordinates (northing, easting) and elevation is what I'm looking for?

Also where can I find C++ example for that scenario ?


texts:
I have a managed raster catalog in a file geodatabase that is comprised of 12 images.  The screenshot shows the raster catalog footprint with the names of the tiff images that make up the raster catalog.  How can I programmatically list the names of the items in a raster catalog?  For example

["tiles_1000_overlap_1_1.tif", "tiles_1000_overlap_1_2.tif"...]

The following approach seems to access the items, although does not list the names.  Similarly, using the arcpy.ListRasters() method returns an empty list.

import arcpy

arcpy.env.workspace = r'C:\temp\myFGDB.gdb\rascat'

rasters = arcpy.ListDatasets()

>>> rasters
[u'rascat\\Raster.OBJECTID = 1',
 u'rascat\\Raster.OBJECTID = 2',
 u'rascat\\Raster.OBJECTID = 3',
 u'rascat\\Raster.OBJECTID = 4',
 ....]


texts: In [Actions] How do I refer to any file with the .qgs ending within the county folder? (rather than need to specify a specific file .qgs?)


texts:
I administer a PostGIS server that was originally intended for internal use. Then, along the way, we open up access to external institutions.
Internally, I therefore created views according to the needs and users. Then for internal webmapping, mapserver allowed me to produce WFS and WMS flows.
Now that we're opening up to the outside world, without much thought, we've started to share mostly WFS streams. Data is visualized using, mostly, GIS clients like QGis and Arcgis.
Before continuing this deployment, I'll try to give you an update on the subject:

Apart from the technological question, what is better to distribute: views or WFS streams?
Are some better than others?
In which case?
In short, I am interested in feedback or reflections more accomplished than mine.


texts: Importing multiple GPX files without selecting waypoints/tracks/... for each file separately using QGIS?

texts:
I've got several GPX files.

I'd like to import them into a QGIS project.

If I'm going with >add Vectorlayer QGIS asks one for each file, which specific parts (waypoints/tracks/...) to load.

I'm looking for a better option.


texts: ST_Intersects not working with ST_Expand from a point

texts:
I would like to do a real simple thing:
I have a bounding box B1 (ST_MakeEnvelope) and I want to see if a square of 100 meters around a lat/long position intersects the other bounding box (B1).

I have tried to do this :

SELECT ST_Intersects(
        ST_MakeEnvelope(2.3400235176086426, 48.87195488537595, 2.3521900177001953, 48.86635145037327, 4326)::geometry,
            ST_Expand(ST_SetSRID(ST_MakePoint(1.888275146484375, 47.89788468044446), 4326)::geometry, 100));

This should return false, because the point position if more than 100 km from my B1 bounding box ...

But I always get :

 st_intersects
---------------
 t

Where is my mistake?

After @zimmi's answer I did :

SELECT ST_Intersects(
            ST_Transform(ST_MakeEnvelope(2.3400235176086426, 48.87195488537595, 2.3521900177001953, 48.86635145037327, 4326)::geometry, 3857),
                ST_Expand(ST_Transform(ST_SetSRID(ST_MakePoint(1.888275146484375, 47.89788468044446), 4326)::geometry, 3857), 100));

This works perfectly.


texts: Re-projection from BNG to WSG84 consistently fails using geopandas' .to_crs

texts:
I am trying to use a shapefile of Great Britain to display some heatmap results on top of it. My data is in EPSG 4326 (longitude, latitude), but the shapefile is in EPSG 27700. I tried re-projecting it, but somehow this goes wrong for England. I have only just started working with geopandas (and spatial data), so maybe this is an obvious mistake.
Shapefile from https://data.gov.uk/dataset/3fd8d2d2-b591-42ff-b333-c53a6a513e96/countries-december-2017-full-clipped-boundaries-in-great-britain
Code:
import pandas as pd
import numpy as np
import geopandas as gpd

UK_map = gpd.read_file(r'C:\Users\Downloads\Countries__December_2017__Boundaries-shp\Countries__December_2017__Boundaries.shp')

UK_map = UK_map.to_crs(epsg=4326)

When I check UK_map with UK_map.crs and UK_map.plot() before transformation, everything seems fine.
But after transformation, England has lots of "inf" in geometry:

Can anyone help?


texts:
I simply want to label a column in the attributes table as: "line 1", "line 2", "line 3" etc. I am sure there's an easy fix but couldn't find the exact solution from other posts.

(using QGIS 2.18.7)

Each row represents a different line feature (GPS track)


texts: OpenLayers 5 forEachFeatureIntersectingExtent Method Not Correctly Selecting Features in Circle or Polygon?

texts:
Using OpenLayers 5.3.3

Draw a circle
Get the extent of the circle in the drawend event:
  const extent = event.feature.getGeometry().getExtent();
The try to determine the features in the extent using forEachFeatureInExtent
The collection of features returned include features that are outside the boundaries of the circle.

The same behavior occurs for an irregularly shaped polygon. However, for a drawn box it appears to work correctly.

Why is the forEachFeatureInExtent not correctly returning the features in the specified extent? What has to be done differently for a circle or irregularly shaped polygon?


texts: Use sf (R) to return the overlapping area between two polygons by rows?

texts:
I have two sfc data frames with the same number of rows (individual IDs), but each data frame contains a polygon assigned to the individual. For example, dataframe1 may have a circle assigned neared an individuals' home coordinate, while dataframe2 has a square assigned near their home coordinate.

In some cases, the square and circle overlap for an individual. In some cases, they do not. I want to know, when they do overlap, what is the area of that overlapping space.

At present I've been using st_intersection, but this, of course, returns the area overlapped by all circles (y)  for an individual's square (x). I've considered a loop as well, but it seems that st_intersection won't work in cases where x and y do not overlap.

Perhaps a loop with a skip-error function built in? Or a more simpler way...?


texts: How to make highlighted polygons when click on Leaflet map from qgis2web plugin?

texts:
I have generated the Leaflet web map using qgis2web plugin in qgis software. The map contain polygons. Now I want the polygons to be highlighted when I click on each and to stay highlighted until I click on another one. I want to know which part of code I should change. I have tried to turn on option "Highlight on hover" when I run qgis2web plugin but that polygon is highlighted just on hover not on click.

Note that I am using qgis 3.0 and the last version of qgis2web


texts: Distorted exported georeferenced tif from far north

texts:
I'm trying to export a tif file for R, to plot my data on top of.

I have geo-referenced this tiff file: "background_changed" from "Background". But when i export it, the file is very compressed in height (because it's from far north) and I do not know how to get around this?


texts: Create a Dynamic Table in your Map Layout

texts:
I am trying to create a graphic table in the layout below are ESRI instruction but I don't have an option to insert the Graphic Table Element (see attached photo

You can only add this table in layout view.

Start ArcMap.
On the main menu, click View > Layout View.
On the main menu, click Insert > Graphic Table Element.

The New Table dialog box appears.
Select Interactive format


texts:
I am working on polygon shapefile and I want to display all the labels in the map. But, I am unable to see the labels of few polygons. I checked my geometry, I almost tried every single permutation and combination mentioned in Layer->Labeling.

Please tell me how to resolve this problem?


texts: Seeing maps created by others in organization in ArcGIS Collector?

texts:
A follow-up question for Editing of layers on an ArcGIS server in ArcGIS Collector

I know it's in ArcGIS Collector's FAQ, but we have met all the requirements specified there and still we can't see others maps within our organization. My co-worker's map is shared with everyone and I can edit it in ArcGIS Online even outside our organization's admin account. However, I can't see it among "all maps" in ArcGIS Collector. The editable layer on his map is shared from a feature service on our organization's server; I have the same layer on one of my maps as well. Saving the layer as a hosted one on that invisible map doesn't make it appear.

So, what may we be doing wrong and how to make all our organization's maps visible and editable by everyone at least in our organization?


texts: GeoDjango Order by distance for Model that is linked to User who has a location on their user profile

texts:
I need to order my query model by distance, however the point field in question is not on the model class but on the user profile and the i am getting an error when attempting the following.

g = geocoders.Google(domain='maps.google.co.uk')
place, (lat, lng) = g.geocode('SW1')
MyModel.objects.distance(Point(lat, lng), field_name='user__userprofile__location').order_by('distance')

The error i get is "ST_Distance output only available on GeometryFields"

Both MyModel and UserProfile have a GeoManager but obviously User does not.

Filtering works fine using

MyModel.objects.filter(user__userprofile__location__distance_lte=(Point(lat, lng), D(mi=distance)))

But I can't understand why sorting is not working.


texts: Geodatabase is registered with the server; new feature class is not

texts: I have an SDE geodatabase with several feature datasets. The geodatabase is registered with the server, and I can publish services without copying the data to the server. I added a new feature dataset, and it is not registered with the server. An ESRI support guy at the UC last year told me that when you add a new feature dataset, it is automatically registered. I have tried adding them as from a database server connection, as a DBO connection. I have tried re-registering the geodatabase, to no avail. Is there any way to register a new feature dataset?


texts:
I'm working through python examples in a geospatial book using jupyter notebook and have run into trouble trying to import gdal_array from my gdal installation.  I've seen several old questions and answers that supposedly fix this with a series of conda uninstalls and re-installs but following the answers gets me nowhere.  I'm running on windows 10 and 64 bit if that makes a difference.  I'm also in Python 3.7 and gdal version 2.3.3.  I've already uninstalled gdal and it's dependencies (numpy etc.) and reinstalled but to no avail.

All the examples use gdal_array.


texts:
I have a button with class search_ and when clicked i want to change the map view to a specific extent. I am using ZoomToExtend API.

Here is what i try, without anything working:

   var roadlayer = new ol.layer.Tile({
    source: new ol.source.OSM()
});
var map = new ol.Map({
    target: 'map',
    layers: [roadlayer],
    controls: ol.control.defaults().extend([new ol.control.ScaleLine()]),
    view: new ol.View({
        center: ol.proj.fromLonLat([33.33386, 35.14710]),
        zoom: 14,
        maxZoom: 19
    })
});
var layerExtent = roadlayer.getExtent();

$('.search_').on('click', function() {
    var ZoomToExtend = function() {
        map.getView().fit(layerExtent, map.getSize());
    }

    try {
        ZoomToExtend();
    } catch (err) {
        alert(err);
    }

});

Errors:

Error 1: Uncaught TypeError: Cannot read property 'getSimplifiedGeometry' of undefined

roadlayer.getExtent(); is undefined. Note that it's the base layer.

Update:
Mike's comment helped but getExtent() only works for the baselayer. What should i use for the dynamically added layer below?

var wmsLayer_Parcel = new ol.layer.Image({
  source: new ol.source.ImageWMS({
    url: NicURL,
    params: {
      'LAYERS': 'parcel_layer'
    },
    serverType: 'geoserver',
   // crossOrigin: 'anonymous',
  })
});

This won't work:

var layerExtent = wmsLayer_Parcel.getExtent() || wmsLayer_Parcel.getSource().getTileGrid().getExtent();

var ZoomToExtent = function(){
map.getView().fit(layerExtent,{
  duration: 1600,
  nearest: true
});


texts: How to create polygons layer from datasheet in QGIS

texts:
I would like to create a polygon shp file from a large dataset. The header looks like this: "Name   MINx    MAXx    MINy    MAXy"
I expect the layer to look like a grid.


texts: Using IF operator with Python parser of ArcGIS Field Calculator?

texts:
I want to assign numbers to a certain types of building and I need to use field calculator. I tried many different versions of this and I still can't find what's wrong.  Could you tell me where do I make a mistake?

]1


texts: Creating symbology on multiple attributes using ArcMap

texts:
I am using ArcMap and want to style datasets based on a number of attributes fields to communicate visually several aspects of the feature.  An example is below, where up to 4 attributes are used to define colour, linetype and symbol.  I do not have the option to split the layer tables into smaller parts due to the use of a defined schema.
Is this possible in ArcMap?


texts: ArcGIS - Python - Extend XY Line to Z Value

texts:
Long time reader, first time poster. So here's my goal:

Working in ArcGIS, trying to work with Model Builder and Python Scripts.

I have a series of lines and a single raster DEM.

I'd like to extend my line segments until their Z-component has changed by at least 20 metres. My initial thoughts on this were:

Write Z_Min and Z_Max for the line
Calculate Z_Diff = Z_Max - Z_Min, for a check
If not, extend the line 30m (random number), check, if not, extend, check, etc.

I've managed to get some code that will run through the process but only stops once ALL features are satisfied (not each individually). I'd like to iterate through each feature and do the check (or check all, extend those that don't satisfy, etc.)

My attempts at querying the datasets via using a Where statement within UpdateCursor function statement were unsucessful - I'm fairly certain I was messing with the created tuples of coordinates (the tuple of OID@ and SHAPE@XY) but I'm not certain.

Current Code - loops through just once and despite my best efforts - more or less resembles some sample script already posted. If anybody has suggestions/samples, that would VERY much be appreciated.

Is this a simple "While" Loop?

#Code from <http://gis.stackexchange.com/questions/71645/a-tool-or-way-to-extend-line-by-specified-distance> by Paul.

from math import hypot
import collections
from operator import add
import arcpy

# Check out any necessary licenses
arcpy.CheckOutExtension("3D")

#Local variables:
DEM = "DEM"

#Need to set these as imported Script for ArcGIS
layer = arcpy.GetParameterAsText(0)
distance = float(arcpy.GetParameterAsText(1))

# Process: Add Surface Information
arcpy.AddSurfaceInformation_3d(layer, EnfDEM_Bird, "Z_MIN;Z_MAX;SURFACE_LENGTH", "BILINEAR", "", "1", "0", "NO_FILTER")
# Process: Calculate Field
arcpy.CalculateField_management(layer, "Z_Diff", "[Z_Max] - [Z_Min]", "VB", "")

#Computes new coordinates x3,y3 at a specified distance
#along the prolongation of the line from x1,y1 to x2,y2
def newcoord(coords, dist):
    (x1,y1),(x2,y2) = coords
    dx = x2 - x1
    dy = y2 - y1
    linelen = hypot(dx, dy)

    x3 = x2 + dx/linelen * dist
    y3 = y2 + dy/linelen * dist
    return x3, y3

#accumulate([1,2,3,4,5]) --> 1 3 6 10 15
#Equivalent to itertools.accumulate() which isn't present in Python 2.7
def accumulate(iterable):
    it = iter(iterable)
    total = next(it)
    yield total
    for element in it:
        total = add(total, element)
        yield total

#OID is needed to determine how to break up flat list of data by feature.
coordinates = [[row[0], row[1]] for row in
           arcpy.da.SearchCursor(layer, ["OID@", "SHAPE@XY"], explode_to_points=True)]

oid,vert = zip(*coordinates)

#Construct list of numbers that mark the start of a new feature class.
#This is created by counting OIDS and then accumulating the values.
vertcounts = list(accumulate(collections.Counter(oid).values()))

#Grab the last two vertices of each feature
lastpoint = [point for x,point in enumerate(vert) if x+1 in vertcounts or x+2 in vertcounts]

#Convert flat list of tuples to list of lists of tuples.
#Obtain list of tuples of new end coordinates.
newvert = [newcoord(y, distance) for y in zip(*[iter(lastpoint)]*2)]

j = 0

with arcpy.da.UpdateCursor(layer, ["SHAPE@XY"], explode_to_points=True) as rows:
    for i,row in enumerate(rows):
        if i+1 in vertcounts:
            row[0] = newvert[j]
            j+=1
            rows.updateRow(row)
            # Process: Add Surface Information
            arcpy.AddSurfaceInformation_3d(layer, EnfDEM_Bird, "Z_MIN;Z_MAX;SURFACE_LENGTH", "BILINEAR", "", "1", "0", "NO_FILTER")
            # Process: Calculate Field
            arcpy.CalculateField_management(layer, "Z_Diff", "[Z_Max] - [Z_Min]", "VB", "")


texts: Loading tiles of EXTENT in OpenLayers locally (localStorage)?

texts:
I have created a simple openlayers 3 map HERE , the map initialization code looks like so:

var styles = [
    'Road',
    'Aerial',
    'AerialWithLabels',
    'collinsBart',
    'ordnanceSurvey'];
var layers = [];
var i, ii;

var coordinate1, coordinate2;
coordinate1 = ol.proj.fromLonLat([77.6013, 12.9044]);
coordinate2 = ol.proj.fromLonLat([77.6885, 12.9536]);

for (i = 0, ii = styles.length; i < ii; ++i) {
    layers.push(new ol.layer.Tile({
        visible: false,
        preload: Infinity,
        extent: [coordinate1[0], coordinate1[1], coordinate2[0], coordinate2[1]],
        myattribute: styles[i],
        source: new ol.source.BingMaps({
            key: 'Ak-dzM4wZjSqTlzveKz5u0d4IQ4bRzVI309GxmkgSVr1ewS6iPSrOvOKhA-CJlm3',
            imagerySet: styles[i]
            // use maxZoom 19 to see stretched tiles instead of the BingMaps
            // "no photos at this zoom level" tiles
            // maxZoom: 19
        })
    }));
}
var map = new ol.Map({
    layers: layers,
    // Improve user experience by loading tiles while dragging/zooming. Will make
    // zooming choppy on mobile or slow devices.
    loadTilesWhileInteracting: true,
    target: 'map',
    view: new ol.View({
        center: ol.proj.fromLonLat([77.6244159, 12.9370699]),
        maxZoom: 20,
        minZoom: 13,
        zoom: 13
    })
});

var select = document.getElementById('layer-select');

function onChange() {
    var style = select.value;
    for (var i = 0, ii = layers.length; i < ii; ++i) {
        layers[i].setVisible(styles[i] === style);
        console.log(layers[i].get('extent'));
    }
}
select.addEventListener('change', onChange);
onChange();

Now my map has a special feature if you noticed , I.E. it has an extent set , I.E. tiles beyond a certain extent won't load, now what i really want to acomplish is, i want the tiles for this extent to be loaded from locally I.E. when the user visits my website the second time around , the tiles sound load form his localstorage , How can this be accomplished ? or can this be accomplished at all , this is a requirement for my project and i have been stuck for some days now on this , can anybody shed some light on this as to how is this possible?


texts:
I have Apache 2.2.22 running on Ubuntu Server 12.04.5 LTS. Following the pycsw install instructions I issued these commands:

sudo add-apt-repository ppa:pycsw/stable
sudo apt-get update
sudo apt-get install python-pycsw pycsw-cgi

After that, if I try to access the address http://localhost/pycsw I get a 403 error in return.

The install instructions seem to be incomplete, what else is needed to get pycsw running on Ubuntu?

Edit: This was just a matter of pointing to the wrong addresses. The main page seems to be in the tests folder:

http://localhost/tests

And the service itself is executed directly through the csw.py file:

http://localhost/csw.py


texts: How to get the data in dbf file while creating shapefile using ArcGIS

texts:
I am using ArcObject DLL to create shapefile. I am successfully able to create the shapefile but I am having some issues with populating the shapefile data into .dbf file. At present the .dbf file just contains the column names but the dataset is empty.

I adding the data in db file as follows

 public IFields CreateFieldsCollectionForFeatureClass(ISpatialReference spatialReference, string ViewType, DataSet ds)
    {
        //IFeatureLayer featureLayer = (IFeatureLayer)annotationLayer;
        IFeatureClassDescription fcDesc = new FeatureClassDescriptionClass();
        IObjectClassDescription ocDesc = (IObjectClassDescription)fcDesc;

        // Create the fields using the required fields method.
        IFields fields = ocDesc.RequiredFields;

        // Locate the shape field with the name from the feature class description.
        int shapeFieldIndex = fields.FindField(fcDesc.ShapeFieldName);
        IField shapeField = fields.get_Field(shapeFieldIndex);

        // Modify the GeometryDef object before using the fields collection to create a new feature class.
        IGeometryDef geometryDef = shapeField.GeometryDef;
        IGeometryDefEdit geometryDefEdit = (IGeometryDefEdit)geometryDef;

        // Alter the feature class geometry type to lines (default is polygons).
        if (ViewType == ViewsType.Powerplant.ToString())
        {
            geometryDefEdit.GeometryType_2 = esriGeometryType.esriGeometryMultipoint;
        }
        else if (ViewType == ViewsType.Pipeline.ToString())
        {
            geometryDefEdit.GeometryType_2 = esriGeometryType.esriGeometryPolyline;
        }
        else if (ViewType == ViewsType.Facility.ToString())
        {
            geometryDefEdit.GeometryType_2 = esriGeometryType.esriGeometryMultipoint;
        }
        geometryDefEdit.HasM_2 = true;
        geometryDefEdit.GridCount_2 = 1;

        // Set the first grid size to zero and allow ArcGIS to determine a valid grid size.
        geometryDefEdit.set_GridSize(0, 0);
        geometryDefEdit.SpatialReference_2 = spatialReference;

        // Because the fields collection already exists, the AddField method on the IFieldsEdit interface
        // will be used to add a field that is not required to the fields collection.
        IFieldsEdit fieldsEdit = (IFieldsEdit)fields;
        foreach (DataColumn col in ds.Tables[0].Columns)
        {
            IField filed = new FieldClass();
            IFieldEdit editFiled = (IFieldEdit)filed;
            editFiled.Type_2 = esriFieldType.esriFieldTypeString;
            editFiled.Name_2 = col.ColumnName;
            editFiled.DefaultValue_2 = "test";//ds.Tables[0].Rows;
            //foreach (DataRow d in ds.Tables[0].Rows)
            //{
            //    editFiled.DefaultValue_2 = d["ID"];
            //}

            //editFiled.DefaultValue_2 = ds.Tables[0].Rows[0].ItemArray[0].ToString();
            editFiled.Length_2 = 50;
            fieldsEdit.AddField(editFiled);
        }
        //fields = (ESRI.ArcGIS.Geodatabase.IFields)fieldsEdit;

        return fields;
    }

How can I get the data of that field in .dbf file.


texts: Can Project Home variable be used to save default outputs of QGIS graphical modeler

texts:
Is it possible to use the @project_home variable in QGIS to set a default file path location for an output in the graphical modeler?

Ideally I would like the default to be the project home folder which contains a sub folder called Shapefiles and then the file name to be the output parameter name from the model.

  @project_home/Shapefiles/[Model output parameter name]


texts: Making legends in Carto Builder reflect custom class breaks for graduated symbols?

texts:
I am trying to build a map using Carto Builder that compares the average price of single family housing and multi-family housing by neighborhood. The data are offset centroids for each neighborhood, sized by the average cost of each type of housing. It is important that each data set use the same class breaks so that a visual comparison is possible. The data sets have different min/max values, so I manually set the class breaks for each using CartoCSS. However, the legend does not reflect the class breaks I set, and returns 'No data available' in the legend box. How do I fix this?

Here is the CartoCSS I am using:

#layer {
  marker-width: ramp([avg_sf], (6,12,18,24,30),
(300000,400000,600000,800000,1000000));
  marker-fill: #e19c0c;
  marker-fill-opacity: 0.79;
  marker-allow-overlap: true;
  marker-line-width: 1;
  marker-line-color: #FFF;
  marker-line-opacity: 1;
}

And a screenshot showing the empty legend.


texts: Joining nodes to network lines ( and not other points) using GRASS

texts:
I'm trying to do a shortest path analysis on my network. In order to prepare the network and attach polygon centroids to the roads, I used GRASS' V.NET operation=connect. the problem with this approach is that instead of connecting points to nearest line, it connect nodes to nearest line OR node which in turn results in unrealistic paths. For example, as you can see in the screen shot, nodes in the central region instead of being connected to the nearest road, are connected to the line that is three time further away.
Is there a way to fix this?
I'm using Grass 6.4 on Mac.


texts: How to convert result of dissolved buffer into individual polygons?

texts: I am wondering if there is a way to easily split results of a disolved buffer into individual polygons? E.g in the screenshot below all of these buffers are one item, can I split them into individual items? I have Arc, MapInfo and QGIS.


texts:
I'm trying to make a style based on rules, and I want to use a boolean field I have created for that.

But I do the following:

"field" like 'TRUE'
"field" like TRUE
"field" like 1
"field" like '1'

and nothing is selected.

How can I select true values in a boolean field then? I'm using QGIS reading from PostGIS database.


texts: Do Different ArcGIS for Desktop License Levels Process Faster?

texts:
At my job, I have access to all three levels of ArcGIS licensing on my computer.

Does the license you are using in ArcMap make a difference on the speed of ArcMap?

I'm pretty sure that the drawing of the map isn't affected, based on tests I've done myself. But, it really does seem like some of the ArcToolbox tools do work faster in ArcInfo. I'm just trying to either validate my suspicious or prove I am going a bit crazy.


texts: External GPS feed for HTML5 geolocation API in django-leaflet?

texts:
My goal is the following:

I'm developing a Django project. The goal is to locate places with a form. I'm using a CreateView generic view (maybe I should create it manually for more customization.

3 optionally ways to locate :

by hand
HTML 5 geolocation API
External GPS (that gives the position to HTML 5 geolocation API or directly in Django)

The first is working well.

For the last point. Would I like to know the best way and HOW?

I read location from serial with pyserial and pynmea
I send to a Linux or cross-platform bridge to give it HTML 5 Geolocation API (external to Django project)

OR

An activation button that switches from HTML 5 Geolocation API to external GPS. Latlong feeds then the leaflet-Django instance.

I know this relevant post for windows but I need make it working on Linux.

use the web geolocation from HTML 5 with the plugin leaflet.location.

Django 2.1 - GeoDJango - Django-leaflet
GPS: BT-747 (serial GPS)


texts: Is there a south-up setting for a projection in QGIS?


texts: Cannot edit ArcGIS data in an SQL Server enterprise geodatabase

texts:
I have created a geodatabase in SQL Server 14. I import some layers to the geodatabase and everything goes ok. I use ArcGIS 10.6.1 Advanced.

The problem is when I try to EDIT these layers. It doesn't let me and the following error appears:

  Not editable layers
  The workspace containing this data cannot be edited.

What is the problem?


texts: Save As Copy from Data Driven Page in ArcGIS Desktop 10.5 to Separate MXD per page in ArcGIS Desktop 10.0?

texts:
I am a beginner in ArcGIS and Python.

I have MXD document created in ArcGIS 10.5. I had created data driven page (DDP) using index shapefile consisted of 35 polygons. I had set name field in DDP from that shapefile with attribute information such as percil-1, percil-2,...percil-35.

I want to save new MXD with the name using the field name attribute but in previous version arcgis 10.0 or 10.1 not in 10.5 MXD version, and I want the final result something like this: percil-1.mxd, percil-2.mxd, percil-3.mxd....percil-35.mxd.

Can I do that in Python?


texts: Backup and restore GWC data with minimal downtime

texts: I have a geowebcache of for an entire country for 20 zoom level have more than 85 Laksh tiles which is saved in C drive ,I need to move the same file into another system with minimum downtime , traditional copy and paste method take more than 2 days and Re caching is also not possible , is there any way to do that using any specific tool


texts: Generating rectangle centered by point in QGIS

texts:
I have a layer which contains points. The layer has height and width attributes.
I need rectangles around those points, but a point should be in the center of the related rectangle whose height and width value comes from the point attributes.
Is there any tool?

Example data for three points:
height   width
20       12
6        10
4        14


texts: Where are raster pixel values stored in PostGIS?

texts:
I have several raster images of temperature, they are stored as geotiffs and I import them to a PostgreSQL database using raster2pgsql.
raster2pgsql -s $SRID -t auto -I $file $SCHEMA.$(basename $file .tiff) | psql -U $USERNAME -d $DATABASE -h $HOST -p $PORT

But when I look at the data using pgadmin I see only an ID and rast column and no temperature values.
Using QGIS I can display original geotiff without any problems.
Any idea what is going wrong?


texts: Calculate distance between several points?

texts:
I'm working on a project that involves calculating distances between apartments and public transport, which I have as two different point datasets.

For this I have used Point Distance tool in ArcMap. It works if I have only one public transport point, i.e. apartment point data as input and only one public transport point as near feature.

However, since there are so many public transport points I want to be able to use the full point dataset in one go.

I'm looking for individual distances from apartments to every unique public transport point, not an aggregated number. I.e. from apartment A there is X meters to point A and Y meters to point B etc.

There must be some way to perform this at the same time, so I won't have to repeat the same thing over and over.

I also tried the Near tools but I don't think they are suited for the task either.


texts: Getting coordinates of user clicked point in current MXD with ArcPy

texts:
I would like to know if it is possible to get the coordinates from a user clicked point in the current map document - arcpy.mapping.MapDocument("CURRENT")?

I am using ArcGIS Desktop 10.0.


texts: How to generate a buffer outside polygons having a given pixel value?

texts:
For example, I would like to build a 20-m buffer (represented in yellow in the image below) outside polygons having a value of 2 (represented in green) and to assign the value of 3 to buffer cells.

I tested the function buffer (package raster) but all buffer cells (including those of original polygons) have the same value.

Here is a reproducible example to build the raster with cell polygons.

pol1 <- rbind(c(-159,-19), c(-158,3), c(-81, 3), c(-148,-75), c(-180,-32))
pol2 <- rbind(c(40,0), c(121,52), c(151,0), c(121,-58), c(81,0))
polys <- SpatialPolygons(list(Polygons(list(Polygon(pol1)), 1),
                              Polygons(list(Polygon(pol2)), 2)))
proj4string(polys) = CRS("+init=epsg:32188") ## NAD83 / MTM zone 8
plot(polys)
r <- raster(ncol=180, nrow=180)
extent(r) <- extent(polys) + 100
projection(r) <- CRS("+init=epsg:32188")
res(r) <- 10
rp <- rasterize(polys, r)
rp[is.na(rp)] <- 0
plot(rp)


texts: Convert x, y, z local points to WGS84 having only 3 coordinates known

texts:
I have a set of (x, y, z) points on a local coordinates system cartesian plan.
For only 3 of these N points, I have the corresponding WGS84 coordinates (latitude, longitude and altitude).
How can I convert each of the other points in lat/lng coordinates?
EDIT:
Here an example of dataset:

Point #
x, y, z
lat, lng, h

P1
0, 0, 0
43.88072176381581, 11.072238294111774, 53

P2
0, 12.24, 0
43.88099406334439, 11.072485923398222, 53

P3
32.42, 6.79, 0
43.88080644977896, 11.072808964924704, 53

P4
-12.4, 1, 0.5
?

P5
55, 8.02, 0.2
?

P6
0.4, 3, 0
?

P7
1.03, -187.4, 5
?


texts:
I wrote a simple script bind to a button

import arcpy
import pythonaddins

class ButtonClass3(object):
    """Implementation for python_add_in_addin.button (Button)"""
    def __init__(self):
        self.enabled = True
        self.checked = False
    def onClick(self):
        import arcpy
        mxd = arcpy.mapping.MapDocument("CURRENT")
        for lyr in arcpy.mapping.ListLayers(mxd):
            if lyr.supports("DATASOURCE"):
            arcpy.AddMessage("Layer: {0} Source: {1} ".format(lyr.name, lyr.dataSource))

But when I click the button, nothing happens. I wonder where could I see the message. I am using ArcMap 10.2.2.


texts:
This may have been asked before, but
A) what software/libraries etc do I need to
B) load multiple rasters into PostGIS?

I have Python 2.7.5, QGIS 1.9 (latest Weekly Build), PostGIS (whatever the latest version is) and I've downloaded the latest batch of GDAL codecs)

Any help would be greatly appreciated

Thanks

Steve


texts: QGIS: Combine proportionte Symbology with graduate fill based on two different attributes

texts:
I would like to style a layer where the fill of the symbols should be graduated e.g. based on an attribute showing the "proportion of the female population" and the size of the symbols based on another attribute e.g. "total population size".

In -> Layer Properties -> Symbology -> Graduated I can choose as a "method" either "size" or "color" (QGIS 3.2.2). However, I would need both.


texts: GeoTag and Import Photo Plugin ERROR what Caused?

texts:
An error occurred while executing the following Python code:

KeyError: 'USERPROFILE'
Traceback (most recent call last):
  File "/home/hermanguse/.qgis2/python/plugins/ImportPhotos/ImportPhotos.py", line 234, in toolButtonImport
    os.path.join(os.path.join(os.environ['USERPROFILE']),
  File "/usr/lib/python2.7/UserDict.py", line 40, in getitem
    raise KeyError(key)
KeyError: 'USERPROFILE'

Python version: 2.7.14 (default, Sep 23 2017, 22:06:14) [GCC 7.2.0]
QGIS version: 2.14.18-Essen Essen, exported

Python Path:
/usr/share/qgis/python/plugins/processing
/usr/share/qgis/python
/home/hermanguse/.qgis2/python
/home/hermanguse/.qgis2/python/plugins
/usr/share/qgis/python/plugins
/usr/lib/python2.7
/usr/lib/python2.7/plat-x86_64-linux-gnu
/usr/lib/python2.7/lib-tk
/usr/lib/python2.7/lib-old
/usr/lib/python2.7/lib-dynload
/usr/local/lib/python2.7/dist-packages
/usr/lib/python2.7/dist-packages
/usr/lib/python2.7/dist-packages/PILcompat
/usr/lib/python2.7/dist-packages/gtk-2.0
/usr/lib/python2.7/dist-packages/wx-3.0-gtk2
/home/hermanguse/.qgis2//python
.
/home/hermanguse/.qgis2/python/plugins/gearthview/ext-libs
/usr/share/qgis/python/plugins/fTools/tools


texts: What is this system of coordinates and how to convert to lat lon with R?

texts:
I am not familiar with geographic coordinates systems.

Initial data

I found these data, and I know that the position is in Paris, but I don't know how to "read" them.

603536, 2429320

Systems of coordinates

With https://epsg.io, I found the approximate location of this point with two very common systems:

EPSG:32633 WGS 84 / UTM zone 33N X: -424216.403492 Y: 5489351.63476
EPSG:4326 WGS 84 : longitude: 2.3845325°, latitude: 48.8622053°

Questions

My questions are :

How to read this: 603536, 2429320 ?
How to convert this to lat long ?


texts: How to divide a feature class into ranges?

texts:
I have a shapefile consisting of around 4000 polygons. The polygons are created from buffering around a similar number of islands somewhere in the Norwegian archipelago. Due to an error in the dissolve functionality using arcpy dissolve fails. I am assuming this is because arcpy.Dissolve_management can not handle too many overlapping features.

To try to work around this I plan on splitting my initial file into ten separate files consisting of around 400 features which are then dissolved separately. These ten files are then merged together and a final dissolve is performed.

How can I do such a selection? The arcpy.select_analysis allows for a where clause. How can I use this to select the n first features, then the n+1 to n*2 features, then the n*2+1 to n*3 features etc?

My final goal is to make a small arcpy function which could be a preliminary fix for us mere mortals relying on ESRIs bug fixing.


texts: Exporting Processing Modeler as Python script in QGIS 3?

texts:
How do I export my model from Processing Modeler as .py so that I can edit its code?

This does not seem as intuitive as in previous QGIS versions.


texts: Point displacement, rule based styling and symbol levels

texts:
So basically I have a point layer with double symbology that looks like this:
type = 1 (big dot red)
 ----status = 1 (small red)
 ----status = 2 (small blue)
 ----status = 3 (small green)
type = 2 (big dot blue)
 ----status = 1 (small red)
 ----status = 2 (small blue)
 ----status = 3 (small green)
type = 3 (big dot green)
 ----status = 1 (small red)
 ----status = 2 (small blue)
 ----status = 3 (small green)

Using symbol levels, it works fine showing the small dot in front of the big dot. The problem is that sometimes people can place multiple type dots overlapping on the same location (or same type). In other words in a single location there can be 3 type3 dots and 2 type1 dots
That's where point displacement comes in. Technically it should work perfectly if I just duplicate the layer using the same symbology. But it does not, even though the symbols levels are already set up properly, the point displacement seems to ignore this and do the following ignoring the inner small dot

The symbol levels are as follows

How to make my double symbology work with point displacement?


texts:
I am trying to calculate an average width of a large set of polygons.

I found the script below researching the topic.
The function is:

def add_average_width(featureClass, averageWidthField='Width'):
    '''
    (str, [str]) -> str

    Calculate the average width of each feature in the feature class. The width
        is reported in units of the feature class' projected coordinate systems'
        linear unit.

    Returns the name of the field that is populated with the feature widths.
    '''
    import arcpy
    from math import pi

    # Add the width field, if necessary
    fns = [i.name.lower() for i in arcpy.ListFields(featureClass)]
    if averageWidthField.lower() not in fns:
        arcpy.AddField_management(featureClass, averageWidthField, 'DOUBLE')

    fnsCur = ['SHAPE@LENGTH', 'SHAPE@AREA', averageWidthField]
    with arcpy.da.UpdateCursor(featureClass, fnsCur) as cur:
        for row in cur:
            perim, area, width = row
            row[-1] = ((perim/pi) * area) / (perim**2 / (4 * pi))
            cur.updateRow(row)

    return averageWidthField

I replace 'featureclass' with the feature I am trying to calculate widths on, and am lost on how to actually execute the script.

I added the script to a toolbox, but am unsure if I am setting the parameters up correctly.

Forgive my Python ignorance - new to programming.


texts: How to draw an exact polygon with exact messurements

texts:
I want to draw a polygons on a Google satellite map. I made a Shapefile to work in, but I can't find the right tool to make a polygon that is exact 20 x 20 or 40 x 10 or some other measure. As I start drawing now its just a messy polygon...

Can somebody give me some advice?


texts:
I have this model:
The tool dissolve has this precondition of the Calculate Tool in order not to run with empty layers:
However, when it runs with a empty layer it returns the following error:

  ERROR 000953: Extent of the overlay operation will be empty. Failed to
  execute (Dissolve).

I put the extent: union of imputs in every tool and in the model.

How can I run the model in this serie of layers (Layers_A),some of them with no features?


texts:
I am just about done with my shapefile in QGIS, but I needed help saving my shapefile as a .prj file. (If that is possible). If not, how do I obtain the .prj file for my shapefile?

Please advise.


texts: How to use Leaflet with Proj4.js when GeoServer is serving tiles in EPSG:404000?

texts:
I have a GeoTIFF served via a GeoServer WMS service and referenced in the EPSG:404000 CRS (Wildcard 2D cartesian plane in metric unit).

I would like to display the tiles in my Leaflet client and keep its original Cartesian coordinates.

I know that I need to use the Proj4Leaflet plugin to use system that are not natively supported by Leaflet.

But the problem is that EPSG:404000 seems to have no Proj4 definition.

How can I specify Proj4Leaflet to use this spatial reference?


texts:
I try to load a WFS layer but

  Analysis of DescribeFeatureType response failed for url  restrictToRequestBBOX='1' srsname='EPSG:3877' typename='akaava:Kaava'
  url='https://opaskartta.turku.fi/TeklaOGCWeb/WFS.ashx' version='auto'
  table="" sql=: it is probably a schema for Complex Features

The URL contains open source information of a city I'm making a map of. The WMS services work just fine but the WFS doesn't work for some reason.


texts:
I am going through File size limitation for shapefile ESRI support article (this article applicable for ArcGIS 9.X version not for 10.X ) I know may be same article will applicable for 10.x version

So any update for 10.1 version??

Purpose: Above question asked from client, side and I explained and given support article just wanted to know updated information for Windows OS (win 7/8) & Linux.Looking in white paper also (10.1)


texts: Clip raster with shapefile won't work without reprojection (Linux)

texts:
So as the title suggests, I am trying to clip an image (UTM) using the OSM land polygon shapefile (WGS 84).

I use gdal_rasterize to do that and when I do it on a windows machine it works flawlessly regardless of the fact that they don't share the same SRS. When I try to do that on Linux I get an error saying that on-the-fly reprojection is not supported and the image is not clipped. This can be solved by reprojecting either of the files, but I can't reproject the image because it will change its dimensions (and I don't want that) and I also can't reproject the shapefile because there are a lot of errors which lead to an empty file.

I really can't find anything on the internet mentioning this issue.

PS. I use CentOS and Windows 10. GDAL in both machines is the newest version.
EDIT: I used the following command to clip the raster with the shapefile:

gdal_rasterize -burn 0 shapefile.shp raster.tif


texts: How to add 2 layers in qgis and display them at the same time

texts:
I added 2 layers in qgis using python but only one layer is shown on canvas. After clicking (anywhere on the canvas) both layers are displayed.

layer =  QgsVectorLayer("Point?field=gid:integer","items","memory")

maplayer=QgsVectorLayer("IND_adm1.shp","IND_adm1.shp","ogr")

QgsMapLayerRegistry.instance().addMapLayers([maplayer,layer])
feature = QgsFeature()
feature.setGeometry( QgsGeometry.fromPoint(QgsPoint(X,Y)) )
feature.setAttributes(values)
layer.startEditing()
layer.addFeature(feature, True)
layer.commitChanges()
layers=[]
layers.append(maplayer)
layers.append(layer)
canvas = QgsMapCanvas()
canvas.setLayerSet([QgsMapCanvasLayer(layers[0])])
canvas.setLayerSet([QgsMapCanvasLayer(layers[1])])
canvas.refresh()
canvas.show()

when it run on qgis it will shows only second layer and when i click on canvas it shows both layer
how could I avoided this one click and both layer should be displayed directly using python.


texts: How to change the color of the layer after searching

texts:
My task is searching a single house from a layer and then zoom to the center. I did this part but i can't clearly identify this house. I am try to do this by adding a fillcolor to the searched area. but i can't.

var point1 = new OpenLayers.Geometry.Point(mins.x, mins.y);
            point1.transform(new OpenLayers.Projection("EPSG:4326"), new OpenLayers.Projection("EPSG:900913"));
            point1.style={fillColor: '#FF0000'};// coloring the point
var point2 = new OpenLayers.Geometry.Point(maxs.x, maxs.y);
            point2.transform(new OpenLayers.Projection("EPSG:4326"), new OpenLayers.Projection("EPSG:900913"));
            point2.style={fillColor: '#FF0000'};// coloring the point
var bounds = new OpenLayers.Bounds();
            bounds.extend(point1);
            bounds.extend(point2);
            bounds.toBBOX();

            map.zoomToExtent(bounds);

but it can't change the previous color
eg: see this page https://gis.fatih.bel.tr/webgis/default.aspx


texts: Seeking USA electric transmission powerlines as shapefiles?

texts:
Does anyone have a good source for USA electric transmission powerlines as shapefiles?

I have some old (1993) data from the US National Renewable Energy Laboratory (NREL) and Cubewerx' WMS (Cubewerx are a spatial data provider); neither of which are very accurate.

Pipeline shapefiles other than the 10 year old NPMS data would also be helpful.


texts:
I am trying to get distance between two polygons and somehow I get different distance, depending on how I call distanceTo

a.distanceTo(b) != b.distanceTo(a) should it be like that?

layers[0].features[0].geometry.distanceTo(
    layers[0].features[1].geometry,
    {details:true,edge:true}
);
Object {distance: 10.48962511638198, x0: 2684839.287184146, y0: 7750111.194373621, x1: 2684849.7767009, y1: 7750111.146694}

and

layers[0].features[1].geometry.distanceTo(
    layers[0].features[0].geometry,
    {details:true, edge:true}
);
Object {distance: 19.05365206557588, x0: 2684857.7419808176, y0: 7749989.676176046, x1: 2684838.7291617, y1: 7749988.4294338}

Some additional information...

Main map object:

var map = new OpenLayers.Map('map', {
    units: 'm',
        projection: 'EPSG:3857',
        layers: [
            new OpenLayers.Layer.Google(
                "Google Roadmap",
                {type: google.maps.MapTypeId.ROADMAP, numZoomLevels: 22}
            ),
            conflictLayer,markerLayer,
        ],
        center: new OpenLayers.LonLat(1787495.895083,  5755860.4398742),
        zoom: 18
    });

later I just add more layers and keep these layers in array, so 'layers[0]' is one of them.
code to generate new layer:

  var newLayer = new OpenLayers.Layer.Vector(o.name, {
      styleMap: styleMap,
      renderers: renderer,
      zone_id:o.id
    });
    map.addLayer(newLayer);
    layers.push(newLayer);


texts: combine two OCG filters in Openlayers : ERROR: invalid input syntax for integer

texts:
i need to combine two OCG filters. each one works fine separably but if i combine them using Filter.Logical.OR, the first one works only.

var c_filter = new OpenLayers.Filter.Logical({
type: OpenLayers.Filter.Logical.OR,
filters: [
new OpenLayers.Filter.Comparison({
    type: OpenLayers.Filter.Comparison.EQUAL_TO,
    matchCase:false,
    property: 'id',
    value:(document.getElementById("Text")).value
    }),
new OpenLayers.Filter.Comparison({
    type: OpenLayers.Filter.Comparison.LIKE,
    matchCase:false,
    property: 'name',
    value: "*" +(document.getElementById("Text")).value + "*"
    })
    ]
});

as you can see, the first one expect integer and the other one a string.
also i tried:

var id_filter = new OpenLayers.Filter.Comparison({
    type: OpenLayers.Filter.Comparison.EQUAL_TO,
    matchCase:false,
    property: 'id',
    value:(document.getElementById("searchText")).value
    });
var name_filter = new OpenLayers.Filter.Comparison({
    type: OpenLayers.Filter.Comparison.LIKE,
    matchCase:false,
    property: 'name',
    value: "*" +(document.getElementById("searchText")).value + "*"
    });

var orFilter = new Array();
orFilter.push(id_filter);
orFilter.push(name_filter);

wfsProtocol.read({
    filter: new OpenLayers.Filter.Logical({
                 type: OpenLayers.Filter.Logical.OR,
                 filters: orFilter
             }),
    callback: processTheQuery,
    scope: strategy
         })

and still, the id_filter works fine. but the name_filter i got the error:
ERROR: invalid input syntax for integer: "test name"
i don't understand why it should be integer !!
and if i try name_filter alone, it works fine !


texts: Drawing contour lines with exact distance?

texts: I am trying to draw contours with exact distances around my samples (points).


texts: Concatenating fields while also replacing values within certain stings

texts:
I am trying to use python to concatenate multiple fields in a polygon layer. Specifically, I am attempting to combine the values of four fields into one:

EXPCODE =

!MATRIX! & !PLANT! & !PCT! & !SITE!

This function seems straightforward to me, however, I would like to replace characters withing the MATRIX field. A few examples of a MATRIX value would be:

H3C
H4A
CS4B
HS4D

After concatenating the EXPCODE field, I need to then write a replace function where all MATRIX "A" and "B" values are combined (H4B becomes H4AB).

Can anyone help me write this as one fluent function? I know I posted this as a python question, but if you are more comfortable with helping me using VBScript I am fine with that.


texts:
I try to use MapServer 7.4 with MapCache 1.6.1 on a Suse SLES
15-system.
MapServer is delivering a WMS very well.
I installed MapCache successfully (I think).
In my mapcache.xml I have set:

<source name="kanaldaten" type="wms">
<getmap>
<params>
<FORMAT>image/jpeg</FORMAT>
<LAYERS>schaechte</LAYERS>
</params>
</getmap>
<http>
<url>http://blablabla/cgi-bin/mapserv?map=/home/blablabla/wms.map</url>
</http>
</source>

   <tileset name="kachelnkanal">
      <source>kanaldaten</source>
      <cache>disk</cache>
      <grid>WGS84</grid>
      <format>PNG</format>
      <metatile>5 5</metatile>
      <metabuffer>10</metabuffer>
      <expires>3600</expires>
   </tileset>

If I use the WMS-demo application:

http://blablabla/mapcache/demo/wms everything works fine.
A couple of tiles are produced on the server.

But: If I use my WMS e.g. in QuantumGIS or at a simple
...GetMap-Request in an URL
no tiles at all are produced on the server??

Am I right that the tile-producing only works, if the mapclient has
specific entries?


texts:
I want to symbolise based on the values, so for example I have many value fields (see pic)
The map is drawn so that all these values have the same symbology and are not distinct. I tried to add all values but this only adds all the attributes of a single field. I want to add all the value fields from NVC_TYPE1 up to RSPB....


texts: Clipping raster is taking too long in 3.16.0-Hannover

texts:
Okay so I wanted to clip [this][1] raster with this [JPN_adm0.shp][2] file. It's been more than 6 hours and it's still on "Processing algorithm part."
I have a decent PC that has 12GB of RAM 1TB NVME SSD and 6GB of DDR6 GPU. I'm pretty sure that I'm doing something wrong on my part.
Here's the log file:
QGIS version: 3.16.0-Hannover
QGIS code revision: 43b64b13f3
Qt version: 5.9.5
GDAL version: 2.2.3
GEOS version: 3.6.2-CAPI-1.10.2 4d2925d6
PROJ version: 493
Processing algorithm…
Algorithm 'Clip raster by mask layer' starting…
Input parameters:
{ 'ALPHA_BAND' : False, 'CROP_TO_CUTLINE' : True, 'DATA_TYPE' : 0, 'EXTRA' : '', 'INPUT' : '/home/rez/Downloads/Japan Merged.tif', 'KEEP_RESOLUTION' : True, 'MASK' : '/home/rez/Downloads/JPN_adm/JPN_adm0.shp', 'MULTITHREADING' : False, 'NODATA' : 0, 'OPTIONS' : '', 'OUTPUT' : '/home/rez/Desktop/Japan Clipped.tif', 'SET_RESOLUTION' : False, 'SOURCE_CRS' : QgsCoordinateReferenceSystem('EPSG:4326'), 'TARGET_CRS' : QgsCoordinateReferenceSystem('EPSG:4326'), 'X_RESOLUTION' : None, 'Y_RESOLUTION' : None }

GDAL command:
gdalwarp -s_srs EPSG:4326 -t_srs EPSG:4326 -of GTiff -tr 0.00027777777777777783 -0.0002777777777777778 -tap -cutline /home/rez/Downloads/JPN_adm/JPN_adm0.shp -cl JPN_adm0 -crop_to_cutline -dstnodata 0.0 "/home/rez/Downloads/Japan Merged.tif" "/home/rez/Desktop/Japan Clipped.tif"
GDAL command output:
Creating output file that is 111794P x 77320L.

Processing input file /home/rez/Downloads/Japan Merged.tif.

The input file is 1.3 GB but the output file is currently 4.3 GB and still processing. What should I do to improve the clip processing time? Should I try using GDAL on cmd?
[1]: https://mega.nz/file/bygC2ZrY#f9qVhTLILM6B_0QxIOSV7z3v-_gPhC09pz1dH18iLfQ
[2]: https://mega.nz/folder/n6gklJDQ#tAlU-G6a5ak6J4ZaVlyYhA


texts: Calculate distance between point and linestring and how to represent that on QGIS?

texts:
I have loaded my data with osm2pgrouting and I have my network with ways (Linestring) and nodes (Point). The problem is : I have a database of doctors addresses (point) and I should calculate the distance between these points and my ways (Linestring).
 : 1
I am using the following query :
SELECT ST_Distance (
ST_GeomFromText('POINT(55.4504545699355 -20.8752577759738)',4326),
ST_GeomFromText('LINESTRING(55.4748349 -27.1342723, 55.4747091 -21.1338699)',4326) )
As ST_DISTANCE;

But I don't think that works. I need to calculate the distance from all the doctors points and represent that in QGIS as a linestring.
Any help or advice?


texts: Installing CartoDB on Raspberry Pi Model B?

texts:
Has anyone gotten CartoDB to work on a Raspberry Pi?

I was able to install everything (minus the PPAs and the two node apps) on my Pi that are in the README. I'm wondering if anyone else has done this or if they can help me with the remaining things. I doubt I'm missing anything with the PPAs, but the two node apps (which basically said things were not okay, couldn't find another app or something, don't recall) worry me.

As I noted in the comments, I'm installing from Github and using the README from there. I'm running on a Model B (512mb RAM). I want to have a development environment to play around with and figured my Pi would work alright (sure it might be a little slow, but that's okay). I'm also running raspbian built mid December 2013.


texts:
Is it possible to add multiple/all layers from a WMS service with PyQGIS?

It is possible to add a WMS layer with PyQGIS if layer name to be added is exactly known. (which is 'MA-ALUS' in the code below)

urlWithParams = 'url=http://kaart.maaamet.ee/wms/alus&format=image/png&layers=MA-ALUS&styles=&crs=EPSG:3301'
rlayer = QgsRasterLayer(urlWithParams, 'MA-ALUS', 'wms')
QgsMapLayerRegistry.instance().addMapLayer(rlayer)

Assuming that if want to add all of the layers in the WMS service, what should be the proper PyQGIS code/implementation?

e.g. code below gives a list of all layers in the WMS if it is executed under QGIS Python Console (so that one can choose and add multiple layers).

qgis.utils.iface.addRasterLayer("https://tucbs-public-api.csb.gov.tr/csb_cdp_abi_wms?&service=wms")

My question is; is it possible to implement similar operation with PyQGIS?


texts: Feature to Raster Conversion: Missing Pixels Near Border

texts:
I am trying to convert a feature to a raster in ArcGIS 10.5. When I do that, it only keeps the pixels where the centroid falls within the boundary of the featureclass. However, I would like to keep any pixel that touches the boundary. I have tried changing environment settings but nothing seems to help.

In the following image the issue can be visualized:

Edit: Here are the settings I am using in the Polygon to Raster tool:


texts: How to set data-driven symbol levels in QGIS?

texts:
Is it possible to drive symbol levels based on an attribute or data element such as area in a QGIS layer?

We have a set of boundaries that include very large and very small sizes.
Manually setting the symbol level to the known order is ok for a single map, but when dozens of maps are to be made (atlas), driving the symbol levels based on an attribute or data driven expression would be valuable.


texts:
I am exporting a raster file (tiff) in ArcMap as a png/pdf. But no matter how much I increase the dpi, the exported map looks smoothed version of the original images in my map.

Any help?


texts: Apply volumes to points on a raster DEM in QGIS

texts:
I am trying to create a flood inundation map in QGIS from known flood volumes at specific points across a raster DEM. I could potentially do this manually for each individual point using contours surrounding the point and estimating the volume of the contours until the contour volume equals the flood volume. However, these would be very tedious, so I'm looking for a tool that would do this automatically.

The attached picture describes what I'm trying to do, except that I am assuming that the water surface elevation is constant for simplicity. I'd apply a volume at a point on a DEM of the ground surface, and would like to extract the flood extent / depths based on that volume.


texts: Python script to check node / edge coincidence

texts:
I have written a script with arcpy to test node and edge coincidence. Each edge (line) should have a node (point) at either end. If this is true, the script does not write anything. If there is node at one end it is marked as closed; if there is not a node at the other end it is marked as open. The script works, but I'd like to optimize it. For 100 edges it is fine, but my dataset has 89,000 edges and it takes over an hour to run:

import arcpy
import os

node_featureclass = arcpy.GetParameterAsText(0)
edge_featureclass = arcpy.GetParameterAsText(1)
folder = arcpy.GetParameterAsText(2)
filename = arcpy.GetParameterAsText(3)

openfile = open(folder + (os.sep) + filename + '.csv', 'a')
openfile.write('PIPEobjectID, ERRORtype1, ERRORtype2' + '\n')

nodeDict = {}
for row in arcpy.da.SearchCursor(node_featureclass, ["SHAPE@XY"]):
    nodeXY = row[1]
    nodeDict[nodeXY] = row[0]

for row in arcpy.da.SearchCursor(edge_featureclass,["OID@", "SHAPE@"]):
    edgeOID = row[0]
    firstPt = row[1].firstPoint
    lastPt = row[1].lastPoint
    firstPointXY = firstPt.X, firstPt.Y
    lastPointXY = lastPt.X, lastPt.Y

for key in nodeDict:
    if nodeDict.has_key(firstPointXY):
        upStreamMH = "CLOSED"
    else:
        upStreamMH = 'OPEN'

    if nodeDict.has_key(lastPointXY):
        downStreamMH = "CLOSED"
    else:
        downStreamMH = 'OPEN'

if upStreamMH == "CLOSED" and downStreamMH == "CLOSED":
   pass
else:
    openfile.write('%i'% edgeOID + ',' + upStreamMH + ',' + downStreamMH + '\n')

openfile.close()

Any suggestions on how it might be re-written to make it run faster?


texts:
As a first step of a pyqgis plugin I'm trying to import several shape files at once. To get them aligned, using the GUI, I'd enable on the fly reprojections and set the project crs to WGS84. I'd then set the layer CRS to WGS84 for each layer resulting in all layers being in the same layer extent.

The following code tries to replicate the GUI steps. However, the end result is that there is still no overlay between the layers, unlike when I do the steps via GUI. What am I missing out?

qgis.utils.iface.mapCanvas().mapRenderer().setProjectionsEnabled(True) # Enable on the fly reprojections
qgis.utils.iface.mapCanvas().mapRenderer().setDestinationCrs(QgsCoordinateReferenceSystem(4326, QgsCoordinateReferenceSystem.PostgisCrsId)) # Set project coordinate reference system to WGS84

for path, name in ShapeFileList:
    vector = QgsVectorLayer(path, name, "ogr")
    if not vector.isValid():
        print "Layer failed to load!"
    vector.setCrs(QgsCoordinateReferenceSystem(4326, QgsCoordinateReferenceSystem.PostgisCrsId)) # Set layer coordinate reference system to WGS84
    QgsMapLayerRegistry.instance().addMapLayer(vector)


texts: What Algorithm is used by ArcGIS Watershed tool?

texts:
Does anyone know what type of algorithm is used in the ArcGIS Watershed tool (in the Spatial Analyst package)?

Very little information given on Esri's website ... but I suspect it may be some kind of depth/breadth search.

I've looked at these ArcGIS Online Help pages:

Watershed (Spatial Analyst)
How Watershed works

So yes, it uses the flow direction raster, but what algorithm is it using to traverse the raster?

Please note, I'm not looking for answers along the lines of 'it uses D8..'...D8 is not really an algorithm, but a model to help define the algorithm you would use. I.E you could implement the D8 scheme within a depth-first search algorithm and/or a breadth-first search algorithm


texts: ESPG:27700 is displaying incorrect coordiantes

texts: I am new to GIS, so forgive me. I have set up a project using Open Street Map as my background, I have set the CRS to ESPG:27700 British National Grid, however the status bar is not displaying the coordinates correctly, and still says 000,000 is located on the intersection of Greenwich meridian and the equator...any thought how to resolve this?


texts: Separating list with coordinates in two different coordinate systems using ArcMap?

texts:
I have an excel list with approximately 10.000 poles which I want to export to a Point layer to view in ArcMap. However, the coordinates are a mix between SWEREF99 and RT90 and I am not sure how to know which of the poles are in which coordinate system.

Is there anyway to know which coordinate system is which? Or is it possible to create the Point layer and correct it in ArcMap?


texts: GeoNode Layer overview Legend not displaying

texts:
The Geonode documentation indicates that once a layer has been uploaded a Legend can be seen, however I am unable to see the legend in the layer overview page nor can I find a reason/error for it not to display.

Is there a setting that needs to be configured during setup in order to see the Legend or can someone maybe point me in the direction of the server log(s) where I might be able to see why this isn't displaying as expected.

Screenshot of Geonode documentation showing the legend

Screenshort of my test server not showing the legend


texts: When is an address in a new subdivision available for mail delivery?

texts:
I have customers that are complaining that some addresses validate on the USPS validation site (https://tools.usps.com/go/ZipLookupAction!input.action) but my geocoding vendor can't geocode the address.   Google Maps shows recent construction in the area, so I'm guessing that the underlying data for the geocoding vendor is out of date.

However, I'm wondering why the postal address is valid.  What is the process by which the USPS validates addresses in new subdivisions?


texts: Spatial shifting of Point feature using Vector Bender-Reg

texts:
I am the new to this QGIS software. I wanted to spatially shift the vector files(Point, Line and Polygon). In that connection, I found this Vector bender plugin and I do. But when I approach Polygon and line features its works perfectly. At the same time, Point features are not working and displaying Python error as "Attribute Error: type object 'Qgis Geometry' has no attribute 'from Point".

How do I do this?


texts: Merging shapefiles to one, rearrange columns?

texts:
I am trying to merge two shapefiles which both contain polygons representing properties in QGIS 2.2.

One problem is that the names of some of the columns in the attribute table differ between the two files. They contain the same information, but two of the column names are different.

I tried making two new columns in one of the files, and copied the information from the ones with the wrong names. Then I deleted the column with the wrong names.

The resault was two attribute tables with the same column names but in different order.
When I then tried merging them, I lost information in the attribute table.

What can I do differently to prevent the loss of data when merging shapefiles?


texts: What would be 'normal' processing time for a large postgis query

texts:
I have two datasets which I want to spatially join (with an intersect).
The first dataset consists of 9 million addresses. The second dataset consists of 12.000 neighbourhoods. I want the add the neighbourhood to the address.

Some information about the system:
Microsoft Windows Server 2003 (SP-2)
Inter(R) Xeon(R) CPU
X5650 @ 2.67 GHz
2.66 GHz, 4,00 GB of RAM

I log in via Remote Desktop

The processing time for this query is about 50 minutes. Quite long I would say.
Can someone tell what would be the general processing time?

Both geometry columns have an index. First I add a column

ALTER TABLE temp.test ADD COLUMN buurtnr INTEGER;

Then I want to fill the column by doing a spatial join.

UPDATE temp.adresses a SET buurtnr = n.buurtnr FROM temp.neighb n WHERE st_intersects (a.geometrie, n.the_geom)

I have to run the query monthly, together with similar query. The whole process now takes about a day. I want to make it faster.


texts: Labeling points numerically in ArcGIS Desktop?

texts:
I just created a points Shapefile, with 300 points.

Using ArcGIS Desktop 9.3 I now need to label all of them numerically.

Is there a way of labelling all of them, instead of me going to the attribute table and label them manually?


texts: Searching geodetic databases for historical landmarks in Canada by Station names (or numbers)

texts:
I am interested in historical landmarks within Canada.  I have a Station name ... Station 08U1209.  I would like to read what information is available about it from our Canadian Geodatic database, but I have tried for hours now to access this database.  I cannot seem to find it.

If it exists, then how do I access it?


texts: Extracting coordinates of placemark from KML using FastKML?

texts:
I would like to parse a KML document in order to put the information (names of placemarks, coordinates) into a python array for further processing. This is what I have so far:

from fastkml import kml

with open('evalPoints.kml') as myfile:
    doc=myfile.read()
k = kml.KML()
k.from_string(doc)

outerFeature = list(k.features())
innerFeature = list(outerFeature[0].features())

placemarks = list(innerFeature[0].features())

for p in placemarks:
    coords = p.coordinates #this does not work

Unfortunately, parsing the coordinates does not work, while parsing the name of the placemarks works like a charm with p.name. What am I doing wrong?


texts: QGIS how can I create a legend with only one set of categories (in common) and two different set of symbols?

texts:
In QGIS I would like to create a legend with two sets of symbols (for the two overlapping layers which are shown in the map) sharing the same categories. Something like that:

Is there a way to do it in QGIS in an organized way without copy and paste?


texts: How to extract long and lat from geometry column in GeoPandas?

texts:
I read a shapefile with GeoPandas and I get a nice table with a column "geometry". This column contains shapely.geometry.polygon.Polygon or shapely.geometry.multipolygon.MultiPolygon object.

import geopandas as gpd
gdf = gpd.GeoDataFrame.from_file("ashapefile.shp")

Now I would like to use this geoDataFrame to make such as a chloropeth plot or something else. Thus I would like to get the (x, y) coordinates, (latitude, longitude) of these polygons.

For example if I do :

for boundary in gdf.boundary:
    print(boundary.xy)

Thus boundary are shapely.geometry.linestring.LineString or MultiLineString. boundary.xy are what I am looking for. But for MultilineString I have to do a loop on all the line.

Is there the way to do that ? Or is there a more simple or straightforward way ?

EDIT :
Issue with unary_union:


texts: I am looking for a Python method to Split the list of attributes in a layer

texts:
I have a list of attributes in QGIS:

[137.0, u'0101S', u'MUNICIPIO', u'28004', 19849.4753178, 22000935.5836, u'\Alamo, El']

And I am trying to obtain each element in the list, using pyqgis.

I can't use a simple loop , like

for feature in features:
    print feature

because it returns the object id.

And if I use this:

for feature in features:
print '\n'.join(feature)

It returns the following error:

TypeError: sequence item 0: expected string, float found

*For example, my aim is:

"value1" = 137.0

"value2" = 0101S

and so on.


texts: QGIS "Controlled by atlas/Fixed scale" makes random scale for features

texts:
I have a point layer as control point in an atlas. I set map to be controlled by atlas, and there's only one option for point layer: fixed scale.

I render it, and every feature gets it's own scale, some of which are completely useless.

Is there a way to fix this thing?

QGIS version 2.18.19

Feature 1:

Feature 2: far away, barely usable

Feature 3: USELESS!!!

Every time same features get the same scale. These three are in the same city, impossible that WGS/Google projection has different scale there within 30 km.

What is going on?!

addition: Just to try, I added a field "sc" and pasted the correct scale there, and set scale in atlas to be used from this field. Did not help at all.

addition 2: just tried moving the map a little bit in composer, and scale changes by several times!


texts: Overpass API returns points instead of polygons

texts:
I'm able to retrieve geojson using the following call to the Overpass API for boroughs in London:

import overpass
api = overpass.API()
response = api.get('relation(1811626);>;out;')

However the result for each geo is individual coordinates, not a polygon:

{"features": [{"geometry": {"coordinates": [25.292698, 56.7083821],"type": "Point"}, "id": 1347523073, "properties": {}, "type": "Feature"}, {"geometry": {"coordinates": [25.292627, 56.7080704], "type": "Point"} ... }}

Is it possible to have Overpass return polygons instead of coordinates when querying relations? Or will we need to convert the coordinates to polygons.


texts:
I am new to ArcGIS, and currently developing a small app that uses ArcGIS JavaScript API.

I have encountered a couple of problems while working on this.

I have about 5000 points that I am returning from a database, parsing them, and displaying on my map (a state in the USA). When my app finishes rendering, the entire state is full of points and you can't make anything out of it. Now, so my understanding, simplemakersymbols, when declared have pixel size, so you can't increase the size dynamically when zooming in.

Is there any alternative to this?

I tried first declaring a smaller size, and then once zoomed in to a certain scale, I redeclare all the points with a larger pixel size value. Hence all this rendering and re-rendering takes a while before all the points are displayed again.

So is there a better way to accomplish this task?

I am using ArcGIS 10.1 and create points based on data I receive from the database on the fly.


texts: Import point cloud (.ply) into PostgreSQL ( PostGIS )

texts:
I am trying to import some point cloud ( .ply file format ) into PostgreSQL Database with PostGIS extension.
How can I do that? Any Idea? Any Importer?


texts: Displaying features by multi-attributes using proportional symbols

texts:
I want to display my features by using multi-attributes. I dont want to use catogories for the symbol size, I would prefer proportional symbols. I know it works pretty well for classifying the symbol size, but I want my second value to be displayed by proportional symbols. Is this even possible or do I just have to be content with categories?
I am using ArcGIS Desktop 10.5.1.
Thanks for the help!


texts: GDAL in C++: Wrong Raster Values After Cropping geoTiff Dataset

texts:
I'm trying to transfer a subset (box (0,0)-(300,300)) of a panchromatic geotiff to a new file.
Although the georeference and the set are transferring correctly, the raster values seem off.

Loading the raster files in qGIS give me at the original .tif values 7614-15959 and at the resulting .tif 0-14530.

Just to clarify: I have no other problems (with projection, georeference or the crop result's borders).

(datasets are GDALDataset and bands are GDALRasterBand)

/*GET PANCHROMATIC BAND*/
    band_pan = dataset_pan->GetRasterBand(1);
/*GET IMAGE/BAND DATATYPE*/
    bufDataType=band_pan->GetRasterDataType();
/*CREATE NEW IMAGE WITH SAME DATATYPE*/
    char **cropOptions = NULL;
    dataset_crop = poDriver->Create(fname_crop,512,512,1,bufDataType,cropOptions);
/*TRANSFER GEOREFERENCE INFO*/
    double adfGeoTransform[6];
    dataset_pan->GetGeoTransform( adfGeoTransform )
    dataset_crop->SetGeoTransform(adfGeoTransform);
    dataset_crop->SetProjection( dataset_pan->GetProjectionRef() )
/*READ IMAGE SUBSET*/
    band_pan->RasterIO(GF_Read,0,0,300,300,buffer,bufsizeX,bufsizeY,bufDataType,0,0);
/*GET NEW IMAGE'S BAND*/
    band_out = dataset_crop->GetRasterBand(1);
/*WRITE SUBSET ON NEW IMAGE*/
    band_out->RasterIO(GF_Write,0,0,300,300,buffer,bufsizeX,bufsizeY,bufDataType,0,0);

Am I missing something?


texts:
Possible Duplicate:
  How can I measure area from geographic coordinates?

I created my spatial tables in WGS-84 with SRID:4326, I want to transform them to UTM Projection to measure length and areas.Plz Explain..
Is it possible to measure area lengths of spatial features in WGS-1984?
another my doubt is in postgis  can we  store coordinates only in decimal degreed system ? is there way to store in utm projection?


texts:
How to write subscript or superscript labels in QGIS 3.4?

Is there any function to do it without any external involvement like Character Map or etc?

I'm EXPLAINING that this is not a duplicate question and I want to use those scripts not only in MAP LAYOUT but everywhere!

For example like this below ▼


texts:
I'm trying to convert a shape geometry to WKB within a arcpy.da.cursor for import it in a Postgis database. When I analyze the result I saw that arcpy give me a 2d geoemtry when I convert to WKB.
Is there anyone who knows this problem??

Am I doing something wrong??


texts: I have Sentinel-2 data (13 bands, spatial resolution 10m/20m/60m) and want to add 30 m SRTM DEM to increase accuracy while classification. The composite prepared should be kept at what spatial resolution?


texts: Forcing label placement to be always inside polygon using QGIS?

texts:
How do I force a label to be ALWAYS INSIDE a polygon with placement method "Using perimeter (curved)"?

I experimented with "Left of line", "Right of line", but some labels are still placed outside the polygon.


texts:
I have a water-polygon grid shape file as shown (Each grid is 1latx1long):

I have created this in QGIS and files (coastal lines and land polygons) selected from http://openstreetmapdata.com/data
Now if I have given a source grid and a destination grid, I have to come up with different routing ways (shortest path, dijkstra or at-least one of them). pgrouting is the way to go, if I am interested in all the different algorithms. These are my questions:
 1. Is there any other simpler way than pgrouting
 2. How can I create the road topology for this grid structure.


texts:
Since Google pulled the plug on Google Maps Classic, I'm reluctantly moving to its Google Maps New application.

However, I didn't find how to export a route to a GPX or KML file so it can be copied onto my smartphone.

Can it do this? If not, is there a third-party solution?


texts:
I have a SQL database that has our local address database on it - it has no geometry within it, just a list of addresses and X and Y co-ordinates. I want to use this data to enable people to search for an address within QGIS and find it on the map?

Any help would be greatly appreciated... otherwise it's a crash-course in python for me!


texts:
I'm a bit lost about NoSQL and applications using GIS.

My aim, to sum it up quickly, is to store geospatial object into a MongoDB database and display them into a map via OpenLayers.
Is it possible to retrieve geospatial data directly from MongoDB and display them into OpenLayers after someone requested it into the application, or have I to use inevitably a server than can handle geospatial data like GeoServer ?

It seems that MongoDB and GeoServer are quite hard to plug, and even if I found a driver it needs further development and I have no time.

If you have some examples, I'll be glad to take a look.
Or If you have some documentation for the driver mentioned above, I'll be glad too :-).


texts: How to get spatial reference in ArcGIS Server identify response

texts:
I'm using the ArcGIS API for JavaScript to query features in a web map. After sending an identify request to my ArcGIS Server 10.2.2 Map Service (added as ArcGISDynamicMapServiceLayer, the response is missing a spatialReference object, which in my application causes the parsed feature to default to WGS 1984 (EPSG: 4326).

This is my request:

http://my.server.ip:6080/arcgis/rest/services/MyMapService/MapServer/identify?f=json&geometry=%7B%22x%22:1452552.552270479,%22y%22:7513080.65512973,%22spatialReference%22:%7B%22wkid%22:102100%7D%7D&tolerance=5&returnGeometry=true&mapExtent=%7B%22xmin%22:1400269.6249234888,%22ymin%22:7501309.352773828,%22xmax%22:1547028.7192308295,%22ymax%22:7531884.164087858,%22spatialReference%22:%7B%22wkid%22:102100%7D%7D&imageDisplay=400,400,96&geometryType=esriGeometryPoint&sr=102100&layers=all:1

Notably, is specifyes using the Spherical Mercator projection (Google/OSM etc) in these three parameters:

sr: 102100,
geometry:{"x":1452552.552270479,"y":7513080.65512973,"spatialReference":{"wkid":102100}}
mapExtent:{"xmin":1400269.6249234888,"ymin":7501309.352773828,"xmax":1547028.7192308295,"ymax":7531884.164087858,"spatialReference":{"wkid":102100}}

The response comes back like this (indented for legibility):

{
  "results": [{
    "layerId":1,
    "layerName":"MyLayer",
    "displayFieldName":"objectTitle",
    "value":"TestTitle",
    "attributes": {
      "geom":"Polygon",
      "objectTitle":"TestTitle",
      "ESRI_OID":"5"
    },
    "geometryType":"esriGeometryPolygon",
    "geometry": {
      "rings": [
        [
          [1454668.8173357714,7517201.4526984785],
          [1455184.6049274686,7515021.5271282298],
          [1454315.3831670529,7513419.0967461951],
          [1453618.9498499723,7512120.8890253156],
          [1452935.2001523944,7510514.1489123777],
          [1452177.4970599019,7510191.5673684804],
          [1450715.0293152332,7515680.0804733364],
          [1453744.8112331014,7517255.1012338968],
          [1454668.8173357714,7517201.4526984785]
        ]
      ]
    }
  }]
}

The same request sent to my ArcGIS Server 10.2.1 server returns the same response but with a spatialReference at the end of the feature:

{
  "results": [{
    "layerId":1,
    "layerName":"MyLayer",
    "displayFieldName":"objectTitle",
    "value":"TestTitle",
    "attributes": {
      "geom":"Polygon",
      "objectTitle":"TestTitle",
      "ESRI_OID":"5"
    },
    "geometryType":"esriGeometryPolygon",
    "geometry": {
      "rings": [
        [
          [1454668.8173357714,7517201.4526984785],
          [1455184.6049274686,7515021.5271282298],
          [1454315.3831670529,7513419.0967461951],
          [1453618.9498499723,7512120.8890253156],
          [1452935.2001523944,7510514.1489123777],
          [1452177.4970599019,7510191.5673684804],
          [1450715.0293152332,7515680.0804733364],
          [1453744.8112331014,7517255.1012338968],
          [1454668.8173357714,7517201.4526984785]
        ]
      ],
      "spatialReference":{"wkid":102100,"latestWkid":3857}
    }
  }]
}

Could there be something wrong with my request or is there a server/service option I've missed?


texts:
I have a growing database of routes. Each route is represented by a polyline with say 200 points on average.

What I'm trying to achieve is automatic tagging of these routes with names of places and areas which each route passes through.

I already have a local installation of stripped-down GeoNames.org database which works well for tagging with administrative areas.

What I am looking for now is the ability to tag them also with geographic features they pass through - like valleys, hills, mountains, forests, natural parks etc.

The problem here is that such features can have quite big surface, so I can't really search for them using their proximity to each route point - I need to know their bounding boxes, which GeoNames does not provide.

My question then is - does anyone know of an API or, even better, downloadable database which:

allows to reverse geocode lat/lon to geographic features,
covers whole world (or at least Europe),
is free to use for non-commercial purposes,
isn't very restrictive in terms of quota limits?


texts: convert dbf file to txt file using arcpy

texts:
I have a series of dBase files and I'd like to convert them to text files  Anyone can help me? I've tried table to table, doesn't seem to work.

The idea is to get ASCII files so I can do further analysis with R script. My current way is to run exportXYv_stats tool, I need a way to avoid manually adding fields as parameter.


texts: Adding new field and populating with shapefile name

texts:
I have a set of shapefiles that need a new field adding to them called "Policy_Ref". I then need to update this field with the shapefile name. However when I do, it uses the full title, including the .shp. However I do not want the .shp in the fields.
This is my code so far:

import arcpy
arcpy.workspace.env = r"teamGIS/shp"
fcs = arcpy.ListFeatureClasses()
for fc in fcs:
    arcpy.AddField_management(fc, "Policy_Ref", "TEXT")
    arcpy.CalculateField_management(fc, "Policy_Ref", "'" + fc + "'")
    arcpy.Merge_management(fcs, "new_shp")

How do I remove the .shp either before or after the code is executed?


texts: Convert discrete data (e.g. geological units) to percentage within pixel?

texts:
I am looking for a way to convert discrete data in a way that they can be used in ecological niche models.  I try to avoid factor data so I thought about transforming the data to %.

I have a shapefile with geological units. I'd like to convert specific units (eg. limestone) to rasters representing the percentage of the respective unit within each pixel. My intention is to use this as an environmental variable in species distribution models.
Here is a similar question I found. However, this results in sharp transition zones. I am wondering if there might be other ways of doing things like that?

I also tried a mooving window approach with focal() in R to smooth the transition a little bit. Additionally, this would give bigger units more "weight" compared to very small isolated ones. This, however,  sometimes results in strange artefacts caused by the window.

Another approach could be to convert the resulting raster from the "aggregate" solution (see link) to a data.frame and use it as input for an interpolation.

What do you think, how could a general workflow look like?

I am working with R and QGIS.

EDIT:

Here are some maps for visualisation.

original shapefile (blue is limestone)

shp rasterized to 10 m and % within 50 m spatial resolution. The values are in %. The transition zone is sharp.

shp rasterized to 50 m spatial resolution and smoothed with mooving window (89*99 pixel). The values are in %, the transition zone is smoothened. However the shape of the mooving window ist visible.

Nearest neighbour interpolation (10 neighbour pixel)

Inverse distance interpolation (IDW) Power:2


texts: Is there any way I can get a list of all the graphics currently on the map? I was thinking that there would be something along the same lines as adding a graphic to a map such as map.graphics.getGraphics(). I am using the ArcGIS API for javascript.


texts: arcpy FromWKT to functioning in_memory layer to be used in geoprocessing tools

texts:
As part of a larger script I am extracting postgis geometries from WKT(got that part down) then I am trying to create an in_memory polygon geometry to be used in arcpy tools.

so far:

    postgis_layer="select ST_Astext(st_force2d(shape)) from {}".format(fc2)
    print postgis_layer
    cur.execute(postgis_layer)
    #good until here!
    in_mry=arcpy.Polygon(arcpy.FromWKT((x for x in cur)))
    print in_mry
    #arcpy.Union_analysis(fc1,in_mry,output)

how do I go from arcpy.FromWKT to the arcpy.Polygon? then can I use that in the union tool? I do not want to save the geometries from postgis as feature classes

    postgis_layer="select ST_Astext(st_force2d(shape)) from {}".format(fc2)
    print postgis_layer
    cur.execute(postgis_layer)
    geoms=arcpy.Array()
    for geom in cur:
        print geom[0]
        geoms.add(arcpy.FromWKT(geom[0]))
    in_mry=arcpy.Polygon(geoms)
    print in_mry
    #arcpy.Union_analysis(fc1,in_mry,output)

select ST_Astext(st_force2d(shape)) from rmp_indicators.total_forest
POLYGON((535836.248384263 740131.124948196,535934.437164262 740343.250836283,535909.998236761 740114.438629866,535836.248384263 740131.124948196))
Traceback (most recent call last):
  File "R:\Data\LUCZ_2017\Python_Scripts\Unclassified_Areas.py", line 53, in <module>
    p.union_fcs('blah','rmp_indicators.total_forest','b')
  File "R:\Data\LUCZ_2017\Python_Scripts\Unclassified_Areas.py", line 43, in union_fcs
    geoms.add(arcpy.FromWKT(geom[0]))
  File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\arcobjects\arcobjects.py", line 48, in add
    return convertArcObjectToPythonObject(self._arc_object.Add(*gp_fixargs(args)))
ValueError: Array: Add input not point nor array object

so I got until here:

    polys=[arcpy.FromWKT(geom[0]) for geom in cur]
    for x in polys:
        print x
    poly=arcpy.Polygon(polys)

polys works, but I do not know how to pass polys into the arcpy.Polygon

x prints out

<geoprocessing describe geometry object object at 0x24FAC2E0>
<geoprocessing describe geometry object object at 0x24FAC320>
<geoprocessing describe geometry object object at 0x24FAC360>
<geoprocessing describe geometry object object at 0x24FAC3A0>
<geoprocessing describe geometry object object at 0x24FAC3E0>
<geoprocessing describe geometry object object at 0x24FAC420>
<geoprocessing describe geometry object object at 0x24FAC460>
<geoprocessing describe geometry object object at 0x24FAC4A0>
<geoprocessing describe geometry object object at 0x24FAC4E0>
<geoprocessing describe geometry object object at 0x24FAC520>

so clearly the geometries are there, I just do not know how to pass them into the Polygon geometry. I am going to just try and use polys as the input to the union tool, maybe the arcpy.FromWKT automatically turns it into a functioning Polygon because it is one


texts: Locating and copying millions of images in QGIS?

texts:
I have image data collected across the US and Canada from 2016 and 2017.  I need to find all images that are within 5 meters of each other in order to compare changes in landscape over the year.  Eg, an image on I-80 East taken in 2016 and a similar image within 5 meters taken in 2017.  Once the images are found I need to copy them into separate folders for further analysis.

All the images are geolocated and there are approximately 20 million to go through.

I am considering using photo2shape to get the location then using a 5 meter buffer to find the images however processing time would be huge considering it takes about 10 minutes to do 100,000 images.  There is also the problem of copying the selected images to another folder which may be the hardest part.  I don't think the output file from photo2shape gives you the ability to copy the images with out opening them in qgis.

Has anybody done anything like this in QGIS?


texts: How can I transfer waypoints and tracks from Base Camp into QGIS2.6?

texts:
I'm working on transferring all my field data from Garmin Base Camp onto QGIS. This is an attempt to overcome the learning curve from basic trail GPS mapping programs to the real thing.

I'd like to transfer large groups of waypoints and tracks into their own vector layers. What's the simplest, most intuitive procedure to do this?


texts: Animating a heatmap over time using Google Maps API

texts:
I'm currently building a web application using flask and google maps api.

I was wondering if it would be possible to create a page where user can select the time using a slider bar and the google maps api will display a animated heatmap of that duration to see how the heatmap clusters move over that period.


texts: OpenLayers 2 WMS CQL_Filter with diacritic

texts:
I´m using GeoServer where I can query WMS services to get image from a CQL Filter query. The URL for the service looks like:

http://myserver/ows/LIMITS_ADM/wms?LAYERS=LIMADM_MUNICIPI&TRANSPARENT=true&FORMAT=image%2Fpng&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&STYLES=&SRS=EPSG%3A25831&CQL_FILTER=(NOM_MUNI%20ILIKE%20%27%25%C3%80ger%25%27%20OR%20NOM_MUNI%20IS%20NULL%20)&BBOX=153185.99759615,4418250,640814.00240385,4818750&WIDTH=1519&HEIGHT=1248

If I run debug on Chrome I get the following:

1. Headers Headersview source:

1.1.    Request URL:
http://myserver/ows/LIMITS_ADM/wms?LAYERS=LIMADM_MUNICIPI&TRANSPARENT=true&FORMAT=image%2Fpng&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&STYLES=&SRS=EPSG%3A25831&CQL_FILTER=(NOM_MUNI%20ILIKE%20%27%25%C3%80ger%25%27%20OR%20NOM_MUNI%20IS%20NULL%20)&BBOX=153185.99759615,4418250,640814.00240385,4818750&WIDTH=1519&HEIGHT=1248

1.2.    Request Method: GET

1.3.    Status Code: 200 OK

1.4.    Remote Address: myserver:80

2.  Response Headersview source

2.1.    Connection: close

2.2.    Content-Disposition: inline; filename=LIMITS_ADM-LIMADM_MUNICIPI.png

2.3.    Content-Type: image/png

2.4.    Date: Wed, 22 Mar 2017 11:52:07 GMT

2.5.    HTTP_GICAR: CODIINTERN=NRDRJN0001;NIF=46335340T;EMAIL=nuria.majo@gencat.cat;UNITAT_MAJOR=CTTI;UNITAT_MENOR=CTTI Qualitat

2.6.    Server: Apache-Coyote/1.1

2.7.    Transfer-Encoding: chunked

2.8.    Via: 1.1 giscorp.sitep.com

3.  Request Headersview source

3.1.    Accept: text/html,application/xhtml+xml,application/xml;q=0.9,image/webp,*/*;q=0.8

3.2.    Accept-Encoding: gzip, deflate, sdch

3.3.    Accept-Language: es-ES,es;q=0.8

3.4.    Authorization: Digest username="XXXXXXX", realm="XXXXXXXX", nonce="MgeAXFBLBQA=1e5891f56cbbe98d13a8861bfca9c2389a1aaaaa", uri="/ows/LIMITS_ADM/wms?LAYERS=LIMADM_MUNICIPI&TRANSPARENT=true&FORMAT=image%2Fpng&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&STYLES=&SRS=EPSG%3A25831&CQL_FILTER=(NOM_MUNI%20ILIKE%20%27%25%C3%80ger%25%27%20OR%20NOM_MUNI%20IS%20NULL%20)&BBOX=153185.99759615,4418250,640814.00240385,4818750&WIDTH=1519&HEIGHT=1248", algorithm=MD5, response="4ca4bf17a155b05c43d2427dd315fd35", qop=auth, nc=00000002, cnonce="72bfb7cba5bac55b"

3.5.    Cache-Control: max-age=0

3.6.    Connection: keep-alive

3.7.    Cookie: XXXXXX; XXXXXX=ca-ES

3.8.    Host: myservice

3.9.    Upgrade-Insecure-Requests:1

3.10.   User-Agent: Mozilla/5.0 (Windows NT 6.1; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/56.0.2924.87 Safari/537.36

4.  Query String Parametersview sourceview URL encoded

4.1.    LAYERS: LIMADM_MUNICIPI

4.2.    TRANSPARENT: true

4.3.    FORMAT: image/png

4.4.    SERVICE: WMS

4.5.    VERSION: 1.1.1

4.6.    REQUEST: GetMap

4.7.    STYLES:

4.8.    SRS: EPSG:25831

4.9.    CQL_FILTER: (NOM_MUNI ILIKE '%Àger%' OR NOM_MUNI IS NULL )

4.10.   BBOX: 153185.99759615,4418250,640814.00240385,4818750

4.11.   WIDTH: 1519

4.12.   HEIGHT: 1248

The response returns an empty image. If I change the cql filter with no diacritic character, it does return an image.

My question is:

Is there a way that I can use diacritic on cql filter?
Is there a way that I can include a header with content type = application/x-www-form-urlencoded;charset=UTF-8?