texts: Disabling custom form fields with python in QGIS

texts:
I have a custom input form and a script to disable the fields that are currently not needed.

from PyQt4.QtCore import *
from PyQt4.QtGui import *

def formOpen(dialog,layerid,featureid):
    e_kon_mang = dialog.findChild(QComboBox,"e_kon_mang")
    e_kon_mang.setEnabled(False)

When I start the editor and then call the input form it works as desired. But if I call the form first and then start the editor, the disabled fields are active again.

As a workaround I even tried to deactivate the fields in the QT Designer, so I have to activate them by script first. But this setting is ignored by QGIS.

Does anyone know what I missed?


texts: Python in field calculator parsing error

texts:
I'm trying to perform a selection in an attribute table based on two columns.
The intital step I'm intending to perform is to select fields between values, e.g. (a > 0 and a < 10) and from this selection then in the second step find fields in another column, e.g. (if b == 1, return "OK" elif b == 2, return "NOT OK"). I think I'm safe with the logic behind it, but the Python scripting seems to be my drawback. What Am I doing wrong since I get parsing error everytime I try to execute the code?

def Check ( a, b ):
 for (a > 0) and (a < 10):
  if b == 2:
   return "NO"
  elif b == 1:
   return "YES"
  elif b == 0:
   return "MAYBE"
 for (a => 10) and (a < 15):
  if b == 2:
   return "NO"
  elif b == 1:
   return "YES"
  elif b == 0:
   return "MAYBE"
 for (a > 15):
  if b == 2:
   return "YES"
  elif b == 1:
   return "NO"
  elif b == 0:
   return "MAYBE"

Result (Check)


texts: Interaction of leaflet.browser.print and esri.geosearch

texts:
I am using leaflet.browser.print:
 L.control.browserPrint({
     documentTitle:'Verify Me!',
     closePopupsOnPrint:false
 }).addTo(map);

with an esri geocode search:
var searchControl = L.esri.Geocoding.geosearch({
providers: [
  L.esri.Geocoding.featureLayerProvider({
    url: 'https://maps.spartanburgcounty.org/server/rest/services/address_points/MapServer/0',
    zoomToResult: false,
    useMapBounds: false,
    searchFields: ['FullName'],
    label: 'Search a Spaco Address',
    formatSuggestion: function (feature) {
      return feature.properties.FullName + ' <br> ' + feature.properties.MSAGComm;
    }
  })

Then I call the result of the search and successfully generate a popup:
 searchControl.on('results', function(data) {
    var popup = L.popup({closeOnClick:false}).setLatLng(data.results[0].latlng).setContent(data.results[0].text).openOn(map);

But even with closePopupsOnPrint set to false, the popup closes when I go to print. I need to it stay visible during the print. I figure there must be some issue between the two leaflet plugins, but I cannot as of yet figure out why. Hoping someone else does.
Here is my pen:
https://codepen.io/pen/?template=jOWmaQb


texts: Hydrological Enforcement of culverts in QGIS

texts:
I work on a farm and we are in the process of developing a surface drainage program for individual fields. I am using QGIS 3.16 and the publicly available 2020 QL2 LIDAR flyover of Indiana, and my goal is to create bluespot maps to help visualize the problem areas, with as much free and open source as possible. I've been able to teach myself quite a lot and have produced maps that are quite usable.
However, these maps are highly affected by culverts, which are plentiful in the flat lands around where the drainage is most important. I have been able to find references to hydrological enforcement and burning streams, but I have not had any luck finding a clear set of instructions on how to burn short culverts into a DEM using QGIS.
For instance, here is an image of a field with a known culvert at the southeast corner running under the road, which I have marked with a green vector (from a previous attempt at burning streams). As you can imagine, the large wet spot in the center of the image doesn't actually exist. I would like to create a culvert bottom in my DEM so that when I calculate the bluespots the algorithm knows to not treat this as a sink.

I have marked quite a few culverts in my dataset that directly affect the drainage we're looking at, so if there's a way to do this as a batch with a shapefile with multiple vectors that would be ideal.
EDIT:
I am currently only working with a DEM mosaic but I should be able to get access to the raw LIDAR download if that is what is necessary.
"Bluespot" is the term that I found to describe a map of basins that do not have outlets, and thus will become wet spots in a field or wetlands/lakes/etc in a catchment area. I am using the process from this video to create these maps, which creates a map of the difference between a filled DEM and the original DEM. The color gradient indicates the depth below the sink border, or in practical terms the severity of the hole. I apologize for the unclear terminology; I am just learning GIS so I have a very weak handle on the jargon.


texts: Activating Data Interoperability extension in ArcGIS for Desktop?

texts: How to solve this problem. I have installed arcgis 10.3 and installed the extensions as well but the extension does not get activated and shows following error.


texts: Geoserver How to call existing functions inside new WPS process

texts:
I am pretty new to Geoserver and programming processes for it. I having trouble to find the information I need.
Is it possible to call a existing function inside my own WPS class? I want to clip some polygons but don't want to import the whole clipping sourcecode into my project.
Do I need some special imports or is it possible to access other WPS processes directly? Would be some sort of process chaining I guess. I can only find examples for this using the WPS builder and XML files but I don't want to go this way, insted do this directly in Java.


texts: How to improve PostGIS ST_Union performance?

texts:
I need to dissolve some large datasets using PostGIS but it takes a lot of time.

What I have already tried:

Lower the precision of the dataset to 6 digits.
Simplify the features
Cluster the table based on spatial index
Adjusted PostgreSQL based on this article: https://postgis.net/docs/performance_tips.html

My datasets are some RF coverage with 5M of features and 2 ~ 5GB of data.

My postgis version is

POSTGIS="2.4.4 r16526" PGSQL="100" GEOS="3.6.2-CAPI-1.10.2 4d2925d6" PROJ="Rel. 4.9.3, 15 August 2016" GDAL="GDAL 1.11.4, released 2016/01/25" LIBXML="2.9.1" LIBJSON="0.11" TOPOLOGY RASTER

And my PostgreSQL version is

PostgreSQL 10.3 on x86_64-pc-linux-gnu, compiled by gcc (GCC) 4.8.5 20150623 (Red Hat 4.8.5-16), 64-bit

Is there a way to improve st_union performance?


texts: Calculating difference of 2 layers elevation

texts: I want to calculate an aquifer thickness. I have bed rock elevation and also I have elevation of surface of earth. How can I calculate difference of elevation of 2 layers=aquifer thickness by ArcMap?


texts:
I want to remove the "Leaflet.com" that appears on the right bottom of the screen
here. How can I do that?


texts: Time estimates on building overviews for mosaic datasets

texts:
I'm trying to create some time estimates for imagery caching in ArcGIS 10.  Does anyone know of any forumlas or ratios to calculate the time it takes to generate overviews for a mosaic dataset?

Does it have anything to do with the area, resolution, file size etc of the actual imagery, or is there even a way to calculate a time estimate for overview generation?

Thanks

UPDATE:
Does anyone have experience in building overviews on massive mosaic dataset (upwards of 9,000 images) in arcGIS 10?

I'm wondering if there is a way to speed up the process ... right now it looks like arcGIS is only using 1 core to build the overviews.  I have a 12 core machine.  Obviously things would go much faster if we can utilize more cores.

My original estimates are blown out of the water ... I'm looking at over 3 weeks to complete this process.

Any ideas?


texts: Difference in the maximum values of a Raster band when using numpy array and when using gdal function directly

texts:
import numpy as np
from gdalconst import*
band = 'band1.tif'
band1_12 = dataset1_12.GetRasterBand(1)
max = band1_12.GetMaximum()
print(max) #this gives 2425.0

data1 = np.array(band1_12.ReadAsArray()).astype(np.float)
data1 = data1.flatten()
max = np.amax(data1)
print(max) #this gives 1977.0

What is the possible cause of error?
Also if I use gdalinfo band1.tif , I see the max value as 2425

I am using gdal version- '1100100' and python 3 on Jupyter Notebook.

Link to file band1.tif -here


texts: Calculating projected area of polygons with field calculator in ArcMap?

texts:
I've got a polygon featurelayer in ArcMap 10.

The coordinate system of the featurelayer's featureclass is web mercator.

The dataframe of the map uses an equal-area projected coordinate system.

I've got a field on the featureclass called Acreage.

I'd like to calculate Acreage based on the polygon's re-projected area (into that of the dataframe's projection).  I don't want to create a new featureclass to do this.

Is there an easy way to do this in ArcMap 10?


texts: How to store route names of a stop belonging to multiple routes

texts:
I've been reading the ArcGIS help files on route analysis to gather more information on what attributes are stored with Stops. Maybe I've stared at the screen too long, but I cannot seem to understand the concept below.

When I Load Locations into the Stops layer, the default value in RouteName is null, as expected. I do not understand what is meant by the second sentence:

  By using this property, stops within one route analysis layer can be
  assigned to multiple routes.

How is this done, if RouteName can only contain one route's name?

I tried to understand the theory by using an example. Two buses BusA and BusB travel in different sections of town, on completely different routes. The only bus stop they both stop at is Stop4, so if someone was on BusA and needed to change routes, they could only get onto BusB at Stop4. When the network dataset was set up, the two bus lines were assigned separate connectivity groups, while the bus stops was assigned to both.

When using the Route solver for this scenario, how is Stop4 assigned a route name, since it is common to both bus lines? Will it only be assigned the route name of one of the buses, but not both?


texts: How to generate a vector map of change from subtracting two rasters

texts:
I am trying to create an erosion map by subtracting two lidar based DEMs from one another and then turning that output into a vector that I can query elevation change from.

One of my first problems is I seem unable to clip the larger raster down to size of the smaller on. I have tried clip by masking where I created a vector from the smaller raster to use as a mask but for some reason the output looked identical to the original large raster.

I had previously subtracted the two rasters and was going to go ahead and vectorize the output, but the processing tool seems to just freeze on processing algorith with the progress bar solid blue before even reaching a % done estimate. I can't tell what is causing this as I have made a few vectors from rasters already.

Assuming I even can get a vector from the subtraction output, I'm not sure how I can get the vector to contain the difference in height as an attribute.

Sorry this is a rambling mess, I am rather new to using qgis and am at a complete loss as to how to solve these issues.


texts: How to convert shapefile to mapinfo file preserving special characters?

texts:
I have to convert a bunch of shapefiles to mapinfo files on a machine without a MapInfo license, so I am using ogr2ogr to perform this task.

The problem is that when opening the mapinfo files in MapInfo Professional, the table's Portuguese characters like Ã are not there any longer. E.g.:

SÃO PAOLO = SÃfO PAULO
BRASÍLIA = BRASÃLIA
GOIÂNIA = GOIÃ,NIA

When opening the same input shapefile through MapInfo Professional (File|Open|Choosing Windows US & Western Europe (ANSI)), I do get the correct result.

I tried lots of different combinations in ogr2ogr using:

SHAPE_ENCODING
-lco ENCODING
intermediate formats (conversion in 2 go's)
etc.

I cannot seem to get it right.
Does anyone have THE solution? :-)


texts: Get Feature Service syncCapabilities with ArcPy

texts:
I'm working on an ArcPy script to test layers. I have either an arcpy mapping module layer object, or a feature service URL as a text string. I wish to get the feature service's capabilities, in particular the syncCapabilities to ascertain if it's sync enabled.

I'm not sure what API to use. I assume you can import JSON in Python, and use it to get the capabilities? The REST API is not very day 1 noob friendly. I could not find a real working example in my context.

  "capabilities": "Create,Delete,Query,Sync,Update,Uploads,Editing"

My hack attempts to integrate it into Python results in failure

ArcGIS REST API: Sync-enabled feature services

import arcpy
import json

URL = "https://services.myserver.com/ERmEceOGq5cHrItq/ArcGIS/rest/services/USA/FeatureServer/jobs/jf90d6386f7494cd59c024f749773fe7b"
print URL
{
    "syncCapabilities" : {
        "supportsSyncEnabled"
}

How do I get the true or false this back from this?


texts:
I have installed QGIS and I would like to upload a .xls-file (or Access database) with 2 GPS points (starting point and ending point), a road name and some distances from the starting point, for instance, every 500m from the starting point until the ending one.

I would like to have a route in QGIS and to give a colour to every single segment (in the example each 500m of the route coloured in a different way). Is it possible? How the .xls-file should be to allow the upload?


texts: What can cause Geoserver to deliver a WMS tile so slowly?

texts:
I have a PostGIS database with a table containing about 300 000 simple polygons (table has spatial index).
I have a Geoserver installed on the same server delivering this table with WMS.

I have a response time of about 6 sec to receive a single WMS tile from that dataset.
(I'm using firebug to get the getMap query url and parameters including the BBOX for a single tile)

When querying that postgis database directly for extracting the same data using the same BBOX with the following query I receive the result in less than 1sec...

SELECT
    a.*,
    ST_Transform(a.geom, 3857)
FROM
    role.my_table a
WHERE
    a.geom && ST_Transform(ST_MakeEnvelope(-8079088.141629988,5721158.693088872,-8077865.1491774265,5722381.685541436, 3857), 4326);

So my question is what can be causing that slow response from Geoserver while PostGIS can extract and transform the data very fastly? Geoserver need 5 sec just to render it? It may worth mentionning that the same Geoserver deliver tiles very fastly from others tables with fewer rows...

Thanks a lot


texts: Modeler - output of first algorithm not going to second one and so on

texts:
I've a model with subsequent steps of LAStools algorithms. Individually they work well, but in the modeler I get only the first output, and no more. It seems like the second algorithm is not taking the output of the previous one.
The log doesn't show any errors. I've only seen that in input (of second step) it shows " " something empty (maybe is here the problem...).
The source is a .laz file, and also the subsequent.
I disabled step 3 and beyond to simplify.


texts: Translate PostgreSQL column names in QGIS

texts:
I have connected to PostgreSQL and imported the layer to QGIS.

The column names are written using English letters.

Can the names of columns be automatically translated to another language?


texts: How to transform coordinates to location in x-y grid in R?

texts:
I have a polygon with a 40m diameter and several points that lay inside it. Now I want to reformat the point locations (projected in DHDN/GK-3, but could vary) to x-y coordinates ranging from zero to the maximal diameter of the polygon (in this case 40m).

I am working in R and therefore I need a scripted solution (this is just a testing subset and many more follow, therefore "click"-solutions in Qgis/arcgis won't suffice.)

Pictures say more than words:
The triangle is the polygons centroid. The green dots the point locations.

I appreciate any coding assistance :-)


texts: Getting the maximum value of a vector field per each unique value of another vector field using PyQGIS?

texts:
I have a vector layer with two columns: Code (string) and Area (float/double). The column Code includes some repeated values. I'd like to evaluate for each unique value of Code the maximum value of the Area values which are associated. Similar questions but with different approaches, not suitable with my need, are at Selecting Maximum value for each unique ID ; How can I create a numeric field (not ObjectID) with a unique number for each record?.
I made some trials but they don't work at all...any suggestions?

1ST STEP
First of all, based on other posts in the forum (Maximum value of vector layer field in PyQGIS), I tried  to implement a script for just calculating the maximum value of the field "Area". It returns 3071.96126727

values = []
    for feat in layer.getFeatures():
       attrs = feat.attributes()
       values.append(attrs[area_index])
    print max(values)

2ND STEP (CODE TRIAL)
Then, I tried to insert the previous script into another one for identifying and iterating on unique values. The problem is that it returns a list of the unique values associated only to the maximum value of the whole field Area i.e. [u'aa1014551', u'ab1014552', etc...] 3071.96126727

uniquevalues=layer.uniqueValues(code_index,limit=10000)
for uv in uniquevalues:
    values = []
    for feat in layer.getFeatures():
        attrs = feat.attributes()
        values.append(attrs[area_index])
    maximum=max(values) #it also has such a temporaly stuck here
    print uniquevalues, maximum


texts: Is it possible to edit or make layers selectable using data from ArcGIS online?

texts:
I am new to using data from ArcGIS online. I am trying to run some analyses using data I obtained from ArcGIS online, specifically the EPA - Ecoregions file of North America listed below. (Note: I have ArcGIS 10.1):

EPA Ecoregions

I cannot edit, select or run tools on the data. The data is listed as "Not Selectable" in the Table of Contents tab.

When using the Set Selectable Layers command, it doesn't appear and, thus, can't change it.

Is it possible to change this setting so I can run analyses using this layer from ArcGIS online? Or is it possible to save this data as a shapefile? Data -> Export Data is not available when right-clicking.

(Note: Access and Use of Contents is listed as NONE on the website, so I believe it's usable)


texts: I got a shapefile in which a column in the attribute table has a specific font (they send me this font file). How can I add this font to ArcGIS to display this column correctly?


texts: Distribution of Human Population with Altitude

texts:
I recently found Joel E. Cohen and Christopher Small's work on hypsographic demography, i.e. the distribution of human population with altitude.

Does anyone know of any other research in this area? I am looking for ideas on how to visualize the relationship between elevation and population in a map and would like to see what has been done already. I could not find any maps beyond those created by Cohen and Small.


texts: How to create images from symbol in ArcObjects?

texts:
I'm trying to create images from symbols so that I can create a custom symbol dialog. Found this sample code on the ESRI forums and it workes great on almost all symbols. I'm having problems with a IMultiLayerFillSymbol that has a IPictureFillSymbol on one of it's layers. It's throwing an exception on symbol.Draw(geometry);

Does anyone know if there is any additional steps that has to taken to get it to work with a IMultiLayerFillSymbol?

static public Image ImageFromSymbol(ISymbol symbol, Int32 width, Int32 height)
{
    try
    {
        Bitmap bitmap = new Bitmap(width, height);
        Graphics g = Graphics.FromImage(bitmap);
        IntPtr hdc = g.GetHdc();

        IGeometry geometry = null;
        if (symbol is IMarkerSymbol)
        {
            IPoint point = new PointClass();
            point.X = (width / 2);
            point.Y = (height / 2);
            geometry = point;
        }
        else if (symbol is ILineSymbol)
        {
            IPolyline line = new PolylineClass();
            IPoint ptFrom = new PointClass();
            IPoint ptTo = new PointClass();
            ptFrom.X = 3; ptFrom.Y = (height / 2);
            ptTo.X = (width - 3); ptTo.Y = ptFrom.Y;
            line.FromPoint = ptFrom;
            line.ToPoint = ptTo;
            geometry = line;
        }
        else
        {
            IEnvelope bounds = new EnvelopeClass();
            bounds.XMin = 1; bounds.XMax = width - 1;
            bounds.YMin = 1; bounds.YMax = height - 1;
            geometry = bounds;
        }

        symbol.SetupDC(hdc.ToInt32(), null);
        symbol.Draw(geometry);
        symbol.ResetDC();
        g.ReleaseHdc(hdc);
        g.Dispose();

        return bitmap as Image;
    }
    catch (Exception ex)
    {
        Debug.Print(ex.Message);
        return null;
    }
}

Update, Added the following code to make it work:

Updated the  ImageFromSymbol metod with:

ITransformation pTransformation = CreateTransFromDC(
                       hDC, size.Width, size.Height);
pSymbol.SetupDC(hDC.ToInt32(), pTransformation);

Added:

 private const int LOGPIXELSY = 90;

[DllImport("GDI32.dll", EntryPoint = "GetDeviceCaps",
           ExactSpelling = true, SetLastError = true)]
       private static extern int GetDeviceCaps(IntPtr hdc, int nIndex);

private static ITransformation CreateTransFromDC(IntPtr hDC,
           int lWidth, int lHeight) {
           // Calculate the parameters for the new transformation,
           // based on the dimensions passed to this function.
           try {
               IEnvelope pBoundsEnvelope = new EnvelopeClass();
               pBoundsEnvelope.PutCoords(0.0, 0.0, (double)lWidth,
                   (double)lHeight);

               tagRECT deviceRect;
               deviceRect.left = 0;
               deviceRect.top = 0;
               deviceRect.right = lWidth;
               deviceRect.bottom = lHeight;

               int dpi = GetDeviceCaps(hDC, LOGPIXELSY);
               if (dpi == 0) {
                   throw new Exception(
                     "Could not retrieve Resolution from device context.");
               }

               // Create a new display transformation and set its properties
               IDisplayTransformation newTrans =
                   new DisplayTransformationClass();
               newTrans.VisibleBounds = pBoundsEnvelope;
               newTrans.Bounds = pBoundsEnvelope;
               newTrans.set_DeviceFrame(ref deviceRect);
               newTrans.Resolution = dpi;

               return newTrans;
           }
           catch {
               return null;
           }
       }


texts: Using custom mapbox studio styles with mapnik

texts:
I created my own tile server according to this instruction and my own layout/style with Mapbox Studio. Following the mentioned instruction I currently use OSM Bright as layout/style.

Question:
How can I use my style on my own tile server (with mapnik)?


texts: Developing "Action" that shows changed value of form-field using PyQGIS?

texts:
for a vector-layer I would like to produce an action, that shows the
changed value in a formular.

E.G. I am editing a vector-point-layer, click on a point, the fomular
is opening. I change a value (lets say
the formular-field "z" had a value "test" and I change it to "test2").

Then I start my developed action.

I tried for example:

from qgis.PyQt import QtWidgets

QtWidgets.QMessageBox.information(None, "Feature id", "feature id is [
$z ] '%Z%' '%z%' $z [z] $Z [Z] [%z%] [%Z%]")

but I always get the "old" value "test" and not the formular-changed
value "test2".

Has anybody an idea what I can do?


texts: Batch clip feature class to preserve area field

texts:
I have a large land cover map in a geodatabase, that I want to batch clip by defined catchments (about 63 polygons). The output files are named based on a field in the "clipper" (sys.argv[2]) called Watershed (currently a hard coded portion of my script).

The script below works if I run it on a shapefile, and point the outputs to an empty folder.
Ideally, I want to run this on a feature class, and point the outputs to a geodatabase using ArcMap9.3, however, I get the following error code: 000210: Cannot create output .

Does anyone have any suggestions?  I have full permissions on my machine, and path directories are correct, so I'm not sure why I'm getting this particular error.

# Parameter list for ArcToolbox
#                                   Parameter Properties
#           Display Name            Data type           Type        Direction   MultiValue
#   argv[1] Output workspace        Workspace           Required    Input       No
#   argv[2] Layer used to Clip      Feature Class       Required    Input       No
#   argv[3] Layer to be clipped     Feature Class       Required    Input       No
#   argv[4] prefix for file name    String          Required    Input       No
#-----------------------------------------------------------------------------

#Import system modules
import sys
import string
import os
import arcgisscripting

# Create the Geoprocessor object
gp = arcgisscripting.create(9.3)

# Load required toolboxes...
gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Analysis Tools.tbx")

# Designate a workspace for cliped vector data
gp.workspace = sys.argv[1]

# Assign variables to be processed...
clipper = sys.argv[2] # feature used to clip vector data
vector = sys.argv[3] # vector data to be clipped and then rasterized
suffix = sys.argv[4] # file naming convention prefix

# Search cursor sequence...
rows = gp.SearchCursor(clipper)
row = rows.Next()
feat = row.Shape # used as the clip_features in the Clip_analysis (clipper input must have a field called "Index_Key", which is used as the file name of each output)
while row:
    n = str(row.Watershed) # assign a variable for the processing message based on a field
    print "clipping to: " +n # tells you what grid was clipped
    # Clip_analysis <in_features> <clip_features> <out_feature_class> {cluster_tolerance}
    gp.Clip_analysis(vector, feat, str(row.Watershed)+"_"+str(suffix)+".shp", "") #update str if necessary
    row=rows.next()

# reset the array...
del rows
del gp


texts: Is it possible to create a distance matrix in parts

texts: Is it possible to create a distance matrix, using Pgrouting in parts? My computer doesn't complete the total nodes(crashes) , about 10.000 ....  I am looking for a way to use parts of the total number of nodes, to complete it later..


texts: Python osgeo/OGR not returning all features in PBF

texts:
I have a PBF file that is a limited export of the Nepal area. I'm trying to pull out all the adminstrative boundaries.

If I use osmconvert to turn my PBF file into a OSM file, sure enough I see the outer national boundary for Nepal in it:
  https://www.openstreetmap.org/relation/184633#map=7/28.416/84.128

<relation id="184633" version="428" timestamp="2020-01-11T00:26:39Z" changeset="79439437" uid="4892967" user="XD346">
        <member type="way" ref="744544230" role="outer"/>
        <member type="way" ref="35550176" role="outer"/>
        <member type="way" ref="45325045" role="outer"/>
        <member type="way" ref="331271373" role="outer"/>
        <member type="way" ref="201928163" role="outer"/>
        <member type="way" ref="331310291" role="outer"/>
        <member type="way" ref="201928153" role="outer"/>
        <member type="way" ref="393979068" role="outer"/>
        ....

However if I build a simple Python script to echo out all the features from the same source PBF file, that relation is not present:

from osgeo import ogr

datasource = ogr.Open("/Users/nicp/nyaruka/posm/nepal/admin_levels.pbf")

layerCount = datasource.GetLayerCount()

for i in range(layerCount):
    layer = datasource.GetLayerByIndex(i)
    print(f"{i}: {layer.GetName()}")

    feat = layer.GetNextFeature()
    while feat is not None:
        print(feat.GetField("osm_id"))
        feat = layer.GetNextFeature()

What's going on here? Is there something I'm missing?

PBF File: https://nyaruka.com/sitestatic/nepal.pbf

Pip Freeze:

GDAL==2.4.4
numpy==1.18.1
psycopg2==2.8.4
PyYAML==5.3
Rtree==0.9.4
Shapely==1.7.0
six==1.12.0


texts:
i am having some trouble connecting to Oracle12c via ArcObjects and EZCONNECT.

Results so far:

•   Can connect via ArcCatalog 10.2.2 using an *.sde file. (EZCONNECT)

•   Can connect via ArcObjects using an *.sde file. (WorkspaceFactory2.OpenFromFile)

•   Cannot connect via ArcObjects using a connection property set. (WorkspaceFactory2.Open)

•   Cannot connect via ArcObjects using a connection string. (WorkspaceFactory2.OpenFromString)

Can someone offer any pointers?

Supporting info:

•   Oracle FAQ on EZCONNECT

http://www.orafaq.com/wiki/EZCONNECT

•   ESRI Blog on Easy Connect

http://blogs.esri.com/esri/supportcenter/2012/10/10/easy-connect-connecting-to-arcsde-just-got-easier/

EDIT - No errors are thrown, this .NET add-in just shuts down ArcMap with no warning or error. Not having luck debugging in VS2012 either.

UPDATE - Well, I think we got it mostly resolved; a misuse of connection properties.
The code snippets below should shed the light.

So can now programmatically pass all of the db connection properties necessary, including the Easy Connect string.

HOWEVER, I still have an issue where the Database Connection dialog appears and must be manually “OK’d” to initiate connection.

The Database Connection dialog was not appearing previously when using the “traditional” connection properties to connect to Oracle11g.

So if you have any thoughts, I am all ears.

/// CODE SNIPPET BEFORE
/// Using "traditional" connection properties.
/// Create an SDE connection property set and populate it with connection properties.
IPropertySet propertySet = new PropertySetClass();
propertySet.SetProperty("SERVER", "server");
propertySet.SetProperty("INSTANCE", "instance"); //this is where i was trying to pass the Oracle Easy Connect string
propertySet.SetProperty("DATABASE", "");
propertySet.SetProperty("USER", "user");
propertySet.SetProperty("PASSWORD", "pass");
propertySet.SetProperty("VERSION", "sde.DEFAULT");
/// Open the geodatabase using the property set.
IWorkspace workspace = workspaceFactory.Open(propertySet, 0);

/// CODE SNIPPET AFTER
/// Using Oracle “Easy Connect” connection string.
/// Create a DB connection property set and populate it with connection properties.
IPropertySet propertySet = new PropertySetClass();
propertySet.SetProperty("DATABASE_CONNECTION_PROPERTIES", "user/pass@server:port/service"); //this is the Oracle Easy Connect string
//propertySet.SetProperty("DBCLIENT", "oracle");
//propertySet.SetProperty("USER", "user");
propertySet.SetProperty("PASSWORD", "pass");
/// Open the database using the property set.
IWorkspace workspace = workspaceFactory.Open(propertySet, 0);


texts:
The question Are Online Certs worth their weight? made me wonder if there are any open source training materials for GIS.

By this I don't mean training specifically for open source software - of which I know there are some - but a set of general training materials for GIS and geographic data that are developed as an open source project.

i.e. a set of training materials that are community developed, freely available, and with an open license. That anyone could take and use for general GIS training.

If there isn't, does anyone think this would be a worthwhile project?

Please add license info for the material you're suggesting.


texts:
I am trying to calculate a con function for a couple hundred rasters (Meridional and Zonal raster data for wind, each being their own raster file). Is there a way to use model builder with the raster calculator and do a loop to go through all the files? I know I need to start off with the iterate raster input; but am lost after that. Also, the files are different: Zonal_04_01_2011_0000 and Meridional_04_01_2011_0000 (for each hour through 2300). each con function uses the zonal and meridional files. If there is a way to do a batch in model builder, when the function continues onto the next files, will it automatically use the files for the 0100 hour meridional and zonal? I hope all this makes sense. I have also looked up ways to do this and a lot of them use python, however, I do not know any python. So if there is a way to compute this without python, that would be great.

using this con function within rater calculator:

Con(("Zonal_4_01_2011_0000.tif" < 0) & ("Meridional_4_01_2011_0000.tif" >= 0), 180- ((ATan(("Meridional_4_01_2011_0000.tif")/("Zonal_4_01_2011_0000.tif"))*57.296), Con(bla bla bla with the same files))

The last four digits in the files are the hour: 0000 is midnight and the files go to 2300 hours.

I need the model builder to go through each of the hours and output a raster file. So after the function goes through the "Zonal_4_01_2011_0000.tif" and "Meridional_4_01_2011_0000.tif" files, the next ones to be used at the same time are "Zonal_4_01_2011_0100.tif" and "Meridional_4_01_2011_0000.tif".


texts: Using lat/long NAD83 in XYtool in Arcgis

texts:
I have loaded a raster map that uses WGS 84. I want to locate stations/points on this map by using lat/long coordinated from an online map that uses NAD83.

Can I just copy and paste the lat/long coordinates from the online map into the XYtool in arcmap?
Does the XY tool directly convert the lat/long to the correct decimal degree position on the map?

I tried a few samples, the stations look to be in proper positions but I am not sure because the difference in WGS84 and NAD83.


texts: Summing values of smaller polygons within larger ones using ArcMap

texts:
I am trying to find the population for Neighborhoods in Chicago. I have a shapefile of the Neighborhoods and a shapefile of census tracts with the population for each tract.
The neighborhoods are larger than the census tracts, so how can I estimate the population of each neighborhood?
I am using ArcMap 10.8.


texts: Add the attributes of selected features to my new layer in PyQGIS

texts:
I am trying to create a layer from the selected points coordinates. I am drawing a circle to the coordinate of selected points. I also need to add initial attributes of my selected features to my new layer. I tried some examples but every time I failed. Can anyone help me with this. Need to add something before "provider.addFeatures( [feature] )"

def draw_circle_to_selected_coordinates(self):
    """
    Gets the coordinate of selected features and draws the circles for desired radius in qgis
    """

    #set active layer as layer
    layer = utils.iface.activeLayer()
    selected_features = layer.selectedFeatures()

    #creating layer structure using selected layer fields.
    fields = layer.pendingFields()
    fields_in="Polygon?"
    field_counter=1

    for field in fields:
        fields_in=fields_in+'field='+str(field.name())+ ':'+str(field.typeName())+'&'
        field_counter+=field_counter

    fields_in= fields_in+('field=Band:Double')

    #creates new memory layer for radial bands with selected layer attributes
    vpoly = QgsVectorLayer(fields_in, "RadialBand", "memory")
    feature = QgsFeature()
    provider = vpoly.dataProvider()
    vpoly.startEditing()

    #get band radiuses from the file
    list=[]
    f=open("radial_band.nodelete", "r")
    with open('radial_band.nodelete') as openfileobject:
        for line in openfileobject:
            value=float(line.strip())
            list.append(value)
    f.close()

    #filling radial band table with selected feature coordinates and desired radiuses
    for f in selected_features:
        geom = f.geometry()
        x=float(geom.asPoint().x())
        y=float(geom.asPoint().y())
        item=0
        for item in range(0, len(list)):
            radius=list[item]
    #draw circle to the selected coordinates and radius one by one
            feature.setGeometry( QgsGeometry.fromPoint(QgsPoint(x,y)).buffer(radius,100))

            #need to add something here to get the values from selected features and add the radius value to the Band field
            provider.addFeatures( [feature] )
            item+=item

    vpoly.commitChanges()
    QgsMapLayerRegistry.instance().addMapLayer(vpoly)


texts: How to load multiple vector layers into QGIS' map canvas using python scripting?

texts:
I am trying to write a standalone script in python to create a QGIS map canvas. The vector layers I am attempting to load are Shapefiles.

Currently, I am able to load each vector layer, add them to the map registry and show them in my map canvas application. The issue I am having is that I need all these vector layers to be shown in the map canvas at the same time. Currently the qgis map canvas just shows the last vector layer that was set to "canvas.setExtent(vector_layer.extent())". I've tried creating a list of all my layers and then putting that into canvas.setExtent(list_of_layers.extent()) but lists don't have the method "extent.()".

I've read the pyQGIS cookbook and I understand that the .extent() method is for a certain vector you would like to see... How can I show all the vector layers at one time (without using Python Console)? I've delved into other Q&As that had similar questions like this one; however, I'm not really finding any solutions that actually address my problem.

Here is my following code:

import modules
from qgis.core import *
from qgis.gui import *
from PyQt4.QtCore import *
from PyQt4 import QtGui
import os, sys
import qgis

def main():

    # supply path to qgis isntall location
    QgsApplication.setPrefixPath("/usr", True)

    # create a reference to the QgsApplication, setting the
    # second arguement to False disables the GUI
    app = QgsApplication([], True)

    # load providers
    app.initQgis()

    # create Qt widget
    canvas = QgsMapCanvas()
    canvas.setCanvasColor(Qt.white)

    # enable this for smooth rendering
    canvas.enableAntiAliasing(True)

    # not updated US6SP10M files from ENC_ROOT
    source_dir = "/home/cassandra/desktop/file_formats/Shapefiles"

    # load vector layers
    for files in os.listdir(source_dir):

        # load only the shapefiles
        if files.endswith(".shp"):
            vlayer = QgsVectorLayer(source_dir + "/" + files, files, "ogr")

            # add layer to the registry
            QgsMapLayerRegistry.instance().addMapLayer(vlayer)

            # set extent to the extent of our layer
            canvas.setExtent(vlayer.extent())

            # set the map canvas layer set
            canvas.setLayerSet([QgsMapCanvasLayer(vlayer)])

    # refresh canvas and show it
    canvas.refresh()
    canvas.show()
    app.exec_()
    app.exitQgis()

main()


texts:
I am trying to add a WMS service to qgis but am having some problems. After adding the url and name of the service and clicking on connect, I am receiving the following error message:-

"Could not understand the response.  The wms provider said:
Download of capabilities failed: Connection timed out"

Does anyone know what the issue might be?


texts: Is there any single command line available which converts all the pcd file in a folder into CSV and create individual CSV file for individual pcd file?


texts: Change coordinate system of image exif data

texts:
I just started using Pix4D mapper and straight away ran into the problem that the flight data is in WGS 84 and the GCPs are in Gauss Krüger 3. The pix4d support suggested to

  Transform the image geotags to the output coordinate system before
  importing them

(see this)

So I wrote a python script to do the job. But after converting WGS84 to GK3 I am unable to write the values into the exif  of a new image.

import os, piexif, pyproj
from PIL import Image

img = Image.open(os.path.join(dirPath,fn))
exif_dict = piexif.load(img.info['exif'])

breite = exif_dict['GPS'][piexif.GPSIFD.GPSLatitude]
lange = exif_dict['GPS'][piexif.GPSIFD.GPSLongitude]

breite = breite[0][0] / breite[0][1] + breite[1][0] / (breite[1][1] * 60) + breite[2][0] / (breite[2][1] * 3600)
lange = lange[0][0] / lange[0][1] + lange[1][0] / (lange[1][1] * 60) + lange[2][0] / (lange[2][1] * 3600)
print(breite) #48.81368778730952
print(lange) #9.954511162420633
x, y = pyproj.transform(wgs84, gk3, lange, breite)
print(x) #3570178.732528623
print(y) #5408908.20172699
exif_dict['GPS'][piexif.GPSIFD.GPSLatitude] = (round(y,6) * 1000000, 1000000)
#exif_dict['GPS'][piexif.GPSIFD.GPSLongitude] = (round(x,6) * 1000000, 1000000)

exif_bytes = piexif.dump(exif_dict) #error here
img.save(os.path.join(outPath,fn), "jpeg", exif=exif_bytes)

I get UnboundLocalError: local variable 'new_value' referenced before assignment  in the dump method.

The GPSInfo tag looks like {0: b'\x02\x03\x00\x00', 1: 'N', 2: ((48, 1), (48, 1), (3449322402, 70000000)), 3: 'E', 4: ((9, 1), (57, 1), (1136812930, 70000000)), 5: b'\x00', 6: (3659, 10)}

I am guessing I have to change the format but how or to what, I do not know.


texts: Converting NZ Topo50 Grid Reference to WGS84 Lat/Long

texts:
newbie here so please pardon me if I'm asking a silly question.

Currently I am required to work with Lat/Long coordinates in WGS84 (EPSG:4326) on QGIS for data analysis.

The input data I have comes in NZ Topo 50 Grid References and I need to transform this into WGS84 Lat/Long coordinates. The tool that I have used to do this successfully previously was the detailed online converter from LINZ found here: http://apps.linz.govt.nz/coordinate-conversion/index.aspx?Advanced=1

However, the online converter only takes in 6 figure grid references (e.g. BA32 123456) whereas I need to convert a more granular 10 figure grid reference (e.g. BA32 123.45 678.91).

Other software from LINZ that I have tried are NZMapConv and Concord. However, NzMapConv only allows up to 8 figure grid references for NZ Topo50 and does not allow bulk conversions whereas Concord does not even contain NZ Topo50 as an option (despite having support for bulk conversions).

As such, are there any tools, software, code that I may have missed that would allow me to easily do bulk conversions from NZ Topo 50 10 figure Grid References into Lat/Long coordinates for WGS84? I'm also open to trying to script my own solution but I was unable to find the algorithm or mathematical formula for the conversion.


texts: Easy way to extract raster cells OUTSIDE a vector (Mask) layer in QGIS?

texts: I have a raster file but am only interested in the cells outside a certain polygon. This polygon is very iregular and scattered across the raster. Is there any easy way to extract only the cells outside of this polygon? A sort of "reverse clip" command? I have zero knowledge of python. I am using QGIS. I need a code-free solution.


texts: QGIS "Fill Ring" tool error on SpatiaLite DB?

texts:
I´m try Fill Ring tool on a SpatiaLite DB and I get a error on the Primary Key.
The Primary Key is integer auto increment.
When I insert the ring and try to save editing I receive error message duplicate key.

Ok, I can insert ring and add a new feature in, but I would like to make most quickly.

This is a problem or I am commiting any error?

Thanks!

version: QGIS Valmiera (Windows x86)

I've tested on Chugiak and I get that the inserted ring is not contained on the feature.


texts: Exporting list of lists of numbers from Google Earth Engine

texts:
I have received an output in GEE code editor in the form of a list of lists of numbers as shown below:

var list1 = ee.List([[1.7,2.4,3.1],[1.3,2.7,3.1],[1.4,2.4,3.2],[1.9,2.8,3.7]]);

I want to export this result in the form of a csv and for this, I tried the following:

Export.table.toDrive(ee.FeatureCollection(list1));

But it throws the error citing this:

Collection, argument 'features': Invalid type. Expected: List<Element>. Actual: List<List<Float>>.

Can you help as I am unable to understand where the problem is?


texts: How to locate a point on google maps/earth using cadastral maps in Kenya

texts:
I'm an architect, so not too well versed with all this datum/transformation talk.

My little problem: i want to be able to enter a position using UTM or deg/sec/min coordinates into google earth and fly there.

the cadastral maps we have come with eastings and northings and a general reference that reads: coordinates referred to Equator and 37* E (37deg E), all in meters

station:A
+X (Easting):116379.70
-Y(Northing) : 6309.10

any help out there??


texts: QGIS corrupted plugins upon reinstallation

texts:
I had various bugs in my GIS softwares (QGIS, ArcGIS) after uninstalling PyCharm from my computer (time limited school license). For a school assignment we had to change the PythonPATH and PythonHOME and after uninstalling the said IDE, some plugins stoped working in QGIS.

I uninstalled the programs, ran CCleaner to remove old registries and just re-installed these:

ArcGIS 10.2
OSGeo4W (including QGIS Lyon 2.12.3)
OpenGeoSuite 4.8 (without QGIS)
Python 2.7
PyScripter

Now the Processing toolbox and other plugins are corrupted and there is this error message:

AttributeError: 'module' object has no attribute 'openssl_md_meth_names'

Also, when the program starts I have other Python errors.

Impossible to load fTools plugin provokes error when calling classFactory()

(translated from french)

Here's the logfile:

Traceback (most recent call last):
              File "C:/OSGEO4~1/apps/qgis/./python\qgis\utils.py", line 306, in startPlugin
                plugins[packageName] = package.classFactory(iface)
              File "C:/OSGEO4~1/apps/qgis/./python/plugins\fTools\__init__.py", line 33, in classFactory
                from .fTools import fToolsPlugin
              File "C:/OSGEO4~1/apps/qgis/./python\qgis\utils.py", line 572, in _import
                mod = _builtin_import(name, globals, locals, fromlist, level)
              File "C:/OSGEO4~1/apps/qgis/./python/plugins\fTools\fTools.py", line 52, in <module>
                import doRandom
              File "C:/OSGEO4~1/apps/qgis/./python\qgis\utils.py", line 572, in _import
                mod = _builtin_import(name, globals, locals, fromlist, level)
              File "C:\OSGEO4~1\apps\qgis\python\plugins\fTools\tools\doRandom.py", line 37, in <module>
                import random
              File "C:/OSGEO4~1/apps/qgis/./python\qgis\utils.py", line 572, in _import
                mod = _builtin_import(name, globals, locals, fromlist, level)
              File "C:\Python27\Lib\random.py", line 49, in <module>
                import hashlib as _hashlib
              File "C:/OSGEO4~1/apps/qgis/./python\qgis\utils.py", line 572, in _import
                mod = _builtin_import(name, globals, locals, fromlist, level)
              File "C:\Python27\Lib\hashlib.py", line 138, in <module>
                _hashlib.openssl_md_meth_names)
            AttributeError: 'module' object has no attribute 'openssl_md_meth_names'


texts: Reducing elevation from DEM using ArcGIS for Desktop?

texts:
I am so curious how would I strip off elevation from the DEM.

I don't know how I would explain it further but I have attached an image to clarify this.

I wanted to eliminate the area shaded with red so that I can have the unshaded area.

The image below is a sample cross section of elevation.

I'm using ArcMap 10.1


texts: What parameters of xarray I should set to upload the NetCDF in proper projection?

texts:
I save an xarray to NetCDF. I have all parameters of the projection and I can write them to xarray.
I tried coord "spacial_ref" and attr 'grid_mapping='spatial_ref' while spatial ref is equal to string with my projection. But the saved file isn't projected and QGis writes "unproper projection".
I load NetCDF as a mesh.
How can I set parameters to get a good file?


texts:
I'm trying for some time to convert a weather radar map (gnomonical projection) to mercator (EPSG_3395) and then layout it onto third-party web map so far unsuccessfully.

The  original image:

The projection can be defined by: +proj=gnom +lat_0=50.008 +lon_0=14.447 +x_0=0 +y_0=0

I can convert pixel coordinates from origin of the projection with ratio 1:1000 by:

cs2cs +proj=gnom +lat_0=50.008 +lon_0=14.447 +to +init=epsg:4326 -f "%.9f"

which gives me those corner coordinates:

    -301500 +217500 top     left:  10.062746833 51.879564163
    +427500 -310500 bottom  right: 20.076498030 47.082819838

My next step is to add metadata to image:

gdal_translate -a_srs "+proj=gnom +lat_0=50.008 +lon_0=14.447" -a_ullr -301500 +217500 729000 528000
            infile.png translated.png

Finally when I try to reproject whole image by:

gdalwarp -overwrite -s_srs "+proj=gnom +lat_0=50.008 +lon_0=14.447" -t_srs EPSG:3395 translated.png warped.png

I get wrong result.
I'm new to GDAL. Where am I wrong and how can I rewrite my commands to make the reprojection go well?


texts: How to change the color of a vector layer in pyQGIS?

texts:
Is it possible to change the color while adding a vector layer? I have a vector layer with water areas and want to load it with a blue color style. The layer is added with:

QgsMapLayerRegistry.instance().addMapLayer(self.vlayer)

Can I somehow change the color after or while loading?


texts: React Google map . Fix the marker in the center

texts:
I want the functionality as Uber/ola app where the marker is fixed in the center evenwhen we drag the map the marker should be fixed in the center. And finally when the drag end we need to get the lcoation of the marker. I am using react js (google-map-react library)

{

import React, { Component } from "react";
import Geocode from "react-geocode";
import { geolocated } from "react-geolocated";
import { Map, InfoWindow, Marker, GoogleApiWrapper } from "google-maps-react";

Geocode.setApiKey("key");

class GoogleMap extends Component {
  state = {
    showingInfoWindow: false,
    activeMarker: {},
    selectedPlace: "",
    markerCoordinate: {
      lat: 0,
      lng: 0
    }
  };
  onMarkerClick = (props, marker, e) => {
    this.setState({
      activeMarker: marker,
      showingInfoWindow: true
    });
  };

  onMapClicked = props => {
    if (this.state.showingInfoWindow) {
      this.setState({
        showingInfoWindow: false,
        activeMarker: null
      });
    }
  };

  onMapReady = () => {
    this.setState({
      markerCoordinate: {
        lat: this.props.coords.latitude,
        lng: this.props.coords.longitude
      }
    });

    this.getAddressfromCoordinates(
      this.state.markerCoordinate.lat,
      this.state.markerCoordinate.lng
    );
  };

  onMarkerDragEnd = coord => {
    const { latLng } = coord;
    const lat = latLng.lat();
    const lng = latLng.lng();
    this.setState({
      ...this.state,
      markerCoordinate: {
        lat: this.props.coords.latitude,
        lng: this.props.coords.longitude
      }
    });
    this.getAddressfromCoordinates(lat, lng);
  };

  getAddressfromCoordinates = (lat, lng) => {
    Geocode.fromLatLng(lat, lng).then(
      response => {
        const address = response.results[0].formatted_address;
        this.setState({ selectedPlace: address });
      },
      error => {
        console.error(error);
      }
    );
  };
  change = (center) => {
    console.log(center);
  };
  render() {
    let body =
loading...
;

    if (this.props.coords) {
      body = (
        <Map
          style={{
            height: "50%",
            width: "50%",
            margin: "0px auto"
          }}
          ref="map"
          google={this.props.google}
          zoom={14}
          onReady={this.onMapReady}
          initialCenter={{
            lat: this.props.coords.latitude,
            lng: this.props.coords.longitude
          }}
          onMouseover={(t,map,e)=>this.change(e)}
          onClick={this.onMapClicked}
        >
          <Marker
            draggable
            animation={this.props.google.maps.Animation.DROP}
            onDragend={(t, map, coord) => this.onMarkerDragEnd(coord)}
            onClick={this.onMarkerClick}
          ></Marker>
          <InfoWindow
            marker={this.state.activeMarker}
            visible={this.state.showingInfoWindow}
            onClose={this.onMapClicked}
          >

{this.state.selectedPlace}

          </InfoWindow>
        </Map>
      );
    }
    return (

        {body}

    );
  }
}

export default GoogleApiWrapper({
  apiKey: "key"
})(
  geolocated({
    positionOptions: {
      enableHighAccuracy: false
    },
    userDecisionTimeout: 5000
  })(GoogleMap)
);
}


texts: API to find polygon for cities and administrative boundaries in India

texts:
I am new to GIS and learning basics by going through various forums. I am looking for a way to get polygon coordinates for any place (city, 2nd/3rd order administrative division) in India. I started with Google Maps but it only returns 2 lat/lng coordinates - NE and SW. I am looking for more points so that I can represent the place in its actual shape on the map.

I came across some sites like Wikimapia, Yahoo Geoplanet, Geonames etc which may serve the purpose but wanted to confirm which one would I should go with. I can write a script to call APIs provided or download one large XML/JSON file which has this info and then store it in my DB.


texts: Whether Sentinel 3 OLCI can be used as a MODIS replacement for 250m NDVI?

texts:
Based on the bandwidths MODIS bands 1 and 2 (620nm - 876nm) will match S3 bands Oa7 - Oa17 (620nm - 865nm)

I am just trying to know if anyone had any experience in using S3 data instead of MODIS for any image processing?


texts:
How can I split large File Geodatabase(FGDB) in Arcgis, It contains many Feature datasets and Feature classes. I want to split it without missing any feature or attribute and after some work I could be able to merge it into one FGDB.
Purpose of splitting is, I have one FGDB and 8 staffs.If I split it into 4 parts it would be helpful to finish job earlier.

Say I have 1 GDB,10 datasets and 200 Feature classes..
GDB 1 should contain 10 datasets and 50 feature classes
GDB 2 should contain 10 datasets and 50 feature classes like that.

Feature datasets should remain same because all 10 datasets contains all type of feature classes.


texts: Where can I download Yahoo's Geoplanet Data?

texts:
I was trying to Download the Yahoo Geoplanet Data, so that my application does not have a dependency on the API provided by Yahoo.

When I go to the download link, I see the following message:

  We are currently making the data non-downloadable while we determine a better way to surface the data as a part of the service.

Since this data was originally shared with a CC BY 3.0 license, there should be someone who is hosting this data, and offers it for download.

Where can I download the actual data?


texts: Openlayers 3.5.0: Problems with GeoJSON in IE11

texts:
I seem to have a problem with loading a geojson-file as source for a layer in IE11 using openlayers 3.5.0. IE11 throws an error without any description. The problem seems to occur when the geojson is set as a vector-source for an imagevector layer.

Both Chrome and Firefox do not have this problem and renders my map as desired.

A little piece of my code:

var groupedData = groupByData(data, level, json.visualConfig.measure);
var minValue = getMinValue(groupedData, 'total');
var maxValue = getMaxValue(groupedData, 'total');
var scaleStepSize = 0;
var scaleSize = json.visualConfig.scaleSize;
var newColors = [];
if(maxValue < scaleSize) {
    scaleSize = maxValue;
    for (var i = colors.length - scaleSize; i < colors.length; i++) {
        newColors.push(colors[i]);
    }
} else {
    newColors = colors;
}
var scaleStepSize = Math.ceil((maxValue - minValue) / scaleSize);

// Reset Extent
featuresExtent = ol.extent.createEmpty();

// Layer
var layer = new ol.layer.Image({
    source: new ol.source.ImageVector({
        source: new ol.source.Vector({
            url: file,
            format: new ol.format.GeoJSON({defaultDataProjection :'EPSG:4326', projection: 'EPSG:3857'})
        }),
        style: (function () {
                return function (feature, resolution) {
                    formatScale(minValue, scaleStepSize, newColors, json.visualConfig.measure + ' per ' + name);
                    var color = '#ffffff';
                    groupedData.forEach(function(vData){
                            if(vData.dimensionValue == feature.get(reference) && vData.dimensionValue != null) {
                                    for(var i = 0; i < newColors.length; i++) {
                                        var startValue = minValue + (i * scaleStepSize);
                                        var endValue = minValue + ((i + 1) * scaleStepSize);
                                        if(vData.total >= startValue && vData.total <= endValue) {
                                            color = newColors[i];
                                        }
                                    }
                                    feature.set('measure', measure);
                                    feature.set('measureValue', vData.total);
                                    ol.extent.extend(featuresExtent, feature.getGeometry().getExtent());
                            }
                    });
                    var newStyle = [new ol.style.Style({
                        fill : new ol.style.Fill({
                            color : color
                        }),
                        stroke : new ol.style.Stroke({
                            color : '#C3C3C3',
                            width : 1
                        })
                    })];
                    return newStyle;
                }
            })()
    }),
    minResolution: minZoom,
    maxResolution: ((layers.length == 1) ? 20000 : maxZoom), // eerste laag mag altijd naar maxResolutie = 20000
    opacity: 0.8,
    level: level,
    reference: reference,
    name: name
});


texts:
I have and SLD with following rule

<se:Rule>
      <se:Name>Primary Pole</se:Name>
      <se:Description>
        <se:Title>Primary Pole</se:Title>
      </se:Description>
      <ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">
        <ogc:PropertyIsEqualTo>
          <ogc:PropertyName>type</ogc:PropertyName>
          <ogc:Literal>type1</ogc:Literal>
        </ogc:PropertyIsEqualTo>
      </ogc:Filter>
      <se:PointSymbolizer>
        <se:Graphic>
          <se:ExternalGraphic>
            <se:OnlineResource xlink:type="simple" xlink:href="svgs/type1.svg"/>
            <se:Format>image/svg+xml</se:Format>
          </se:ExternalGraphic>
          <se:Size>16</se:Size>
        </se:Graphic>
      </se:PointSymbolizer>
    </se:Rule>

and for creating the base64 representation of this rule i am using below code

 // assuming we  will provide only point
 BufferedImage image = null;
 image = Glyph.point(rule);
 ByteArrayOutputStream bos = new ByteArrayOutputStream();
 ImageIO.write(image, "png", bos);
 Base64.Encoder encoder = Base64.getEncoder();
 String imageString = encoder.encodeToString(bos.toByteArray());

but issue is the when i show it on client when height and width is 16 image is displayed fine but when i increase the height and width  the  generated base64 image is blured (pixelated)  what i think its generated the image w.r.t the size in rule is there any option to ignore it and provide my own size? Or is there a better way to generate the legend.


texts:
I have successfully created a query for a feature layer. However, if I run the query while zoomed in, it only filters those graphics within view. When I zoom out, the areas that were not originally in view are not affected by the query.

Here is the site in question:
Public Meeting Map

Here is the code


texts: Creating symbol from EMF acts like picture and blocks out surrounding vector data in ArcMap?

texts:
I have a lot of CAD blocks that I want to convert to symbol markers to use in ArcGIS 10.  Everything I read says the emf is converted to vector data, but that does not seem to be the case.

Procedure: open Style Manager>Open the Style I created for the blocks>Click on Marker Symbol folder>Right click for new marker symbol>The only option that allows emf is the Picture Marker option. When I bring the emf picture in it's the block on what is effectively a piece of paper and when I use it as my point symbol, the surrounding area is blanked out.  Changing the tier level in the TOC is not a feasible option, and I think that is only a bandaid and not a solution.

I think what I really need is a character marker but I don't see an option for creating them from an emf file.

What am I doing wrong?


texts: define a precision accuracy QGIS project

texts:
How to define a precision accuracy of decimal places in  QGIS project or where I use digitizing to add new features ?

for example I have a point or polygon and true coordinates of vertex have to many precision like this 14521.124545421 , 454411.4545441414 I want only precision 3 of digits.I say the true coordinates to take some point where I add in QGIS I don't mean in attribute table using X or Y there I can easy to define precision.
first I think to change precision is project - project properties - general - precision but that is only for display....that is easy work in ARCGIS and I want to QGIS.


texts:
I am trying to remove the legend associated with the current CartoDB generated map and replace it with a legend from a map generated when the user selects a radio button. My code is as follows:

function showMap(radioButton, currentLayer, currentMap, mapUrl) {

        if (radioButton.checked == true) {
            if (currentMap.viz != null) {
                var vis = currentMap.viz;
                vis.legends.getLegendByIndex(0).remove();
            }

            cartodb.createLayer(currentMap, mapUrl)
            .addTo(currentMap)

            .on('done', function(layer) {
                //do stuff
                currentLayer=layer

                .on('featureOver', function(e, latlng, pos, data) {
                    console.log(e, latlng, pos, data);
                })
            })

            .on('error', function(err) {
                alert("some error occurred: " + err);
            });
        } <!--end function-->
        else {
            currentLayer.hide();
        }
    }


texts: How to add a color table to a one band tiff using GDAL?

texts:
I have 1-band GeoTIFF file. The color interpretation for band is gray. I want a color table in order to transform my Col.tiff to Col.jpeg or Col.png.

I tried putting a color table in this tiff via a .vrt:

gdal_translate -of VRT  Col.tiff Col.vrt

then adding in the .vrt:

<ColorInterp>Palette</ColorInterp>
<ColorTable>
   <Entry c1="0" c2="0" c3="255" c4="240"/>
   <Entry c1="0" c2="0" c3="254" c4="240"/>
   <Entry c1="0" c2="0" c3="253" c4="240"/>
   <Entry c1="0" c2="0" c3="252" c4="240"/>
</ColorTable>

And translate again :

gdal_translate Col.vrt Col.tiff

But it answered:

Warning 1: Unable to export color table to GeoTIFF file.  Color tables can only be written to 1 band or 2 bands Byte or UInt16 GeoTIFF files.

How can I do it?

For information, there is gdalinfo for the Band:

Band 1 Block=55x18 Type=Float64, ColorInterp=Gray
  Description = Something
  Metadata:
    GRIB_COMMENT=Something[p]
    GRIB_ELEMENT=TMP
    GRIB_FORECAST_SECONDS=0 sec
    GRIB_REF_TIME=1401926400 sec UTC
    GRIB_SHORT_NAME=SMTHG
    GRIB_UNIT=[p]
    GRIB_VALID_TIME=1401926400 sec UTC


texts: Selecting and deselecting markers in leaflet under R shiny?

texts:
I would like to write a script that selects a marker on the map whenever I click on it and then deselect it when I click/select another marker (that is, so that no more that one marker is selected on the map at any time). I managed to get the selection bit working for the multiple "current" marker, but not the deselection of the last clicked marker, partly because I am confused with the time ordering/reactivity of objects in shiny environment. I attempted to use leafletProxy() to substitute default icons with selected icons to mimic selection, but that only took me this far.

library(shiny)
library(leaflet)

ui <- fluidPage(
   fluidRow(leafletOutput("clinics.map"))
)

server <- function(input, output, session) {
     icons <- makeIcon(iconUrl = "clinic_icon.png", iconWidth = 27, iconHeight = 30)
     icons.sel <- makeIcon(iconUrl = "clinic_icon_selected.png", iconWidth = 27, iconHeight = 30)

   output$clinics.map <- renderLeaflet({
      leaflet(data = clinics) %>%
      addProviderTiles("CartoDB.Positron") %>%
      setView(lat = 37.900000, lng = -121.900000, zoom = 9) %>%
      addMarkers(lng = ~lon, lat = ~lat, icon = icons, layerId = ~ID)
   })

   previous.markers <- reactiveValues(markerID = c(), marker.lat = c(), marker.lon = c())

   observeEvent(input$clinics.map_marker_click, {
      marker.click <- input$clinics.map_marker_click
      previous.markers$markerID <- c(previous.markers$markerID, marker.click$id)
      previous.markers$marker.lat <- c(previous.markers$marker.lat, marker.click$lat)
      previous.markers$marker.lon <- c(previous.markers$marker.lon, marker.click$lng)

      last.marker.id <- tail(previous.markers$markerID, 1)
      last.marker.lat <- tail(previous.markers$lat, 1)
      last.marker.lon <- tail(previous.markers$lon, 1)

      leafletProxy("clinics.map", session) %>%
      removeMarker(last.marker.id) %>%
      #addMarkers(last.marker.id,      # This bit doesn't work b/c there's no last marker initially
      #           lng = last.marker.lon,
      #           lat = last.marker.lat,
      #           icon = icons) %>%
      addMarkers(marker.click$id,
                 lng = marker.click$lng,
                 lat = marker.click$lat,
                 icon = icons.sel)
   })

shinyApp(ui, server)


texts: Clip data frame using Outline of Selected graphic(s) in ArcGIS Desktop?

texts:
I did a subtraction of a data frame in ArcGIS using the following steps.
I prepared a layout in ArcGIS. There is a lot of empty space at the top of layout page. So I want to remove this space. I used a graphic operation to remove it.

Created a graphic rectangle that the same size as the Data Frame
Created another rectangle covering empty space
Selected rectangle 2, then rectangle 1, and chose--> Graphic Opperations --> Subtract.
In the last step, I wanted to clip the subtracted image by "Outline of Selected graphic". But 'Outline of selected graphic' is not activated. I can't click it. I'm wondered why this is not activated?


texts: scaling of results: the comparing of osm-files for a few days - a statistical approach

texts:
i am pretty new to Geographic-information Systems.

i am interested in gathering information on schools Some days before i have heard that it is interesting to have a closer look at the taginfo-site (see below):

What does this mean https://taginfo.openstreetmap.org/tags/amenity=school

here we have two requests on the above site:

a. yesterday: 2014-06-23 23:58 UTC

Type count Objekts all

**all:**  592672 0.02% -

[Node] Node 325997 - 0.36%
[Way] Way 260585 - 0.11%
[Relation] Relation 6090 - 0.23%

b. some (guess two or three) days before: we have had the following values:

countl Objekte all: 592 443 0.02%

[Node] Node 325 974 0.36%
[Way] Way 260 384 0.11%
[Relation] Relation 6 085 0.23%

see the details:

today all:  592672 0.02% -
some days before: countl Objekte all: 592 443 0.02%

you see: that differs: we have a differende of 229 x (entries, records or some what)

Does this mean that there is a total growth of records of schools or that some records have been changed!?


texts: Creating groups of points from lat/long pairs using R?

texts:
I have a database which contains Lat/Long pairs to identify the location of points of interest.  I would like to group the points of interest into groups of 10.  The group should be geographically local and contain exactly 10 points.  Each group should be of a minimum area.

I've looked at various implementations in R but none of them (that I can see) allow you to specify a definite cluster size.

I previously asked Grouping map points into fixed cluster sizes? but I don't think I was exact enough in my question to get a good answer.

Geographically local - I think I mean that groups shouldn't significantly overlap.  In my application (allocating people to groups for monitoring purposes) it would be ideal if each group was as small as possible in physical area.
Minimum area - again, trying to keep the group area to a minimum.  I suppose this could be quantified as keeping the each group's area below a specified threshold (to avoid dozens of small groups and one large one).


texts: Linking Feature Classes between Two Different Map Frames so that Definition Queries in one Frame Apply to Another Frame

texts:
In ArcGIS Pro 2.4.1, would like to determine whether it is possible to link two identical feature classes, where feature class 1-1 is in map frame 1-1 and feature class 1-2 is in map frame 1-2. The ultimate end goal would be that a definition query placed in feature class 1-1 would automatically apply to feature class 1-2, and vice versa.

I have attempted to simply copy paste and, as expected, this does not produce the desired result. I have also tried to find a setting within the feature class properties but nothing appears to relate to a function like this.


texts: R apply function in turn to each subset of a spatialPolygonsDataFrame

texts:
I'm trying to get the area of all the polygons in a SpatialPolygonsDataFrame: parcProj.

library("rgeos")
gArea(parcProj[2,])

works fine but:

apply(parcProj,1,gArea)

or

sapply(parcProj@polygons,gArea)

don't.

Any standard way to do an apply besides the obvious:

for(i in 1:length(parcProj)) areas <- c(areas,gArea(parcProj[i,]))


texts:
I'm trying to parse WFS reponse using OpenLayers.Format.GML parser. However, an exception thrown in my try-catch block is TypeError: a is undefined.

protocol = new OpenLayers.Protocol.WFS({
    version: "1.1.0",
    url: WFS_HOST,
    featureType: "pl",
    featureNS: "ws",
    srsName: "EPSG:900913"
});

wfsLayer = new OpenLayers.Layer.Vector("PE Vector", {
    strategies: [new OpenLayers.Strategy.BBOX()],
    protocol: protocol
});

var _Callback = function(resp) {
    console.log(resp);
    console.log(resp.error);

    try {
        var gmlParser = new OpenLayers.Format.GML();
        gmlParser.extractAttributes = true;
        var features = gmlParser.read(resp.responseText);

        if(features){
            for(var i = 0; i < features.length; i++) {
                var feature = features[i];
                                    // More operations
            }
        }
    } catch(e) {
        alert("Error: " + e);
    }
};

response = protocol.read({
    maxFeatures: 100,
    callback: _Callback
});

My environment: GeoServer 2.1.3, OpenLayers 2.11, PostGIS 1.5.


texts: Error to obtain coordinate reference system argument from a raster image in R

texts:
I got an error when i wanted to obtain the CRS from stacklayer. I don´t know what could it be. I need some advice or help to obtain it.
 This is my code.

#LOAD RASTER LIBRARY
>library(rgdal)
>library(raster)
>ndvi<-raster('ndvi2005.tif')
>ndvi2<-raster('ndvi2006.tif')

#RASTER STACK
>stk <- stack(ndvi1, ndvi2)
#DISPLAY STK
>stk
class       : RasterStack
dimensions  : 1160, 1573, 1824680, 2  (nrow, ncol, ncell, nlayers)
resolution  : 0.008627311, 0.008627319  (x, y)
extent      : -74.49244, -60.92168, -20.00769, -10  (xmin, xmax, ymin,    ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
names       : ndvi2005, ndvi2006
min values  :       NA,       NA
max values  :       NA,       NA

then I changed the NA values

>stk[stk==0]<- NA
>stk
class       : RasterBrick
dimensions  : 1160, 1573, 1824680, 2  (nrow, ncol, ncell, nlayers)
resolution  : 0.008627311, 0.008627319  (x, y)
extent      : -74.49244, -60.92168, -20.00769, -10  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names       : ndvi2005, ndvi2006
min values  :     -0.3,     -0.3
max values  :   0.9990,   0.9993

then when i wanted to obtain de coordinate reference system with CRS() i got the next message

>cord<-CRS(stk)
>Error in !is.na(projargs) && !nzchar(projargs) :
 invalid 'x' type in 'x && y'

I know that my projection is +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 but i have in total 100 stacklayers in differents proyections that i need to analize each one, but i dont know what is my mistake or is something wrong with my code because then i want to use the output of CRS(stk)in the function raster(...., CRS=cord) as an input.


texts: Detecting "load" or "done" event after setting new CartoCSS using set() in cartodb.js

texts:
I'm loading a map data tilelayer onto a Leaflet map using cartodb.

Upon user input, the map reconfigures its CartoCSS to display data differently. Example below.

layer.on({
    click: function(e){
        dataLayer.getSubLayer(0).set({
            cartocss: "#export { polygon-opacity: 0; }"
        }).on("load", function(){ console.log("load"); });
     }
});

I'd like to be able to detect when the reconfigured layer has loaded -- there's sometimes a considerable lag -- with a "done" or "load" event. But neither works as drafted above, and the documentation doesn't give any hints as to any special promise/event for .set().


texts:
I am trying to populate a grid with the number of features in each grid cell. I've gotten this script to work with a subset of the data, but when I try to run it on the full (52,000+) cell feature class, it freezes when doing the Field Calculation and won't progress past it. I want to run it outside of Arc so it processes faster than the 20 minutes it took to run inside Arc.

I'm using Identity to split these features by the grid cells, so I get an accurate count, then spatial joining that to the Grid, then joining the FID from the spatial join to the original grid and field calculating the Join Count over.  This is done for several feature classes, so I need it to calculate into the original Grid FC.

Additionally can arcpy do feature layer outputs like running the script inside Arc?

Whenever I have my input as "Input" I get an error 000210 cannot create output.

arcpy.MakeFeatureLayer_management(Grid, "Grid_Lyr")
Feature_Split = "C:\Users\Testing.gdb\Feature_Split"
Feature_Join = "C:\Users\Testing.gdb\Feature_Join"
print "Counting Features"
arcpy.Identity_analysis(Ways, Grid, Feature_Split)
print "Identity Analysis Done"
arcpy.SpatialJoin_analysis("Grid_Lyr", Feature_Split, Feature_Join , "JOIN_ONE_TO_ONE", "KEEP_ALL", "", "CONTAINS")
print "Spatial Join Done"
arcpy.AddJoin_management("Grid_Lyr", "OID", Feature_Join, "TARGET_FID", "KEEP_ALL")
print "Join Added"
arcpy.CalculateField_management("Grid_Lyr", "Feature_Count", "!Feature_Join.Join_Count!", "PYTHON_9.3", "")
arcpy.RemoveJoin_management ("Grid_Lyr")
print "Fields Calculated"


texts:
Is there a function to crop sf map object, similar to maptools::pruneMap(lines, xlim= c(4, 10), ylim= c(10, 15)) used for SpatialPolygon or SpatialLine?

I am considering st_intersection() but there can be proper way.


texts: I need to identify the area (m2) of each cell in a raster layer. I do know how to find the cell size under metadata (which in my case is [0.00980834,-0.00748333]); but how do I convert this into m2 or km2?


texts: Creating multiple LAS datasets from LiDAR files in various directories using model builder?

texts:
I have 60 folders each containing approximately 20 LiDAR files in .las format. I want to create one LAS dataset from .las files from each folder, and each LAS dataset is named by the folder name. I thought it would be straightforward to use iterator in model builder (ArcGIS 10.1) but unfortunately I couldn't find the suitable iterator, as it doesn't seem to iterate by folders. I tried "iterate files" and "iterate multivalue" but, struggle to produce individual LAS dataset by the folder name.

Could anyone give me some ideas or tips to complete my task?


texts: QGIS : how to get selected feature "state" to filter data in QGIS Composer

texts:
In a QGIS layer, I have selected some entities.

Is it possible to show these in a particular way using style definition using a "state" like $currentfeature or @atlas_featureid ?

The goal is to gather automatically only selected values in the composer map and attribute table and set composer map extents to the selection extent.

This would enable the user to simply select and print its data without caring much about printed map extent and printed attribute table synchronisation.


texts: Drawing lines between points recorded through GPS device in ArcGIS

texts:
I have collected data (GPS Coordinates) for an 11 KV Distribution line (pole to pole GPS coordinates) through a GPS device. The distribution line consists of different kind of conductors (i.e. ACSR Dog conductor, ACSR Rabbit Conductor between different points) and electrical transformer of difference rating (i.e. 50 KVA, 100 KVA, 200 KVA). The GPS points have unique IDs. I have also Data on paper which specifies the location of transformers, S/E, kind of conductor between different points.
My problem is how to draw a line between these GPS points with different line segments specifying kind of conductor (ACSR DOG, ACSR RAB between different points) from the attribute table. Also I want to give specific symbols for transformers located at different GPS points.

i am using ArcMap. i have tried points to line which simply connects all the points. i need segments with different properties. I am new to ArcGis. I need the concept to do the work i have mentioned in question


texts:
I'm looking to make a python add in tool for ArcMap to draw a polygon.  Eventually I want the tool to also calculate the area and a bunch of other things but for now I would be happy just drawing an irregular polygon by clicking the mouse for each vertex and double clicking to close.  Anyone have a template of what that would look like as I'm horrible at coding?

The genericode produced by the add in wizzard gives the following:

class ToolClass3(object):
"""Implementation for Test Addin_addin.tool (Tool)"""
def __init__(self):
    self.enabled = True
    self.shape = "NONE" # Can set to "Line", "Circle" or "Rectangle" for interactive shape drawing and to activate the onLine/Polygon/Circle event sinks.
def onMouseDown(self, x, y, button, shift):
    pass
def onMouseDownMap(self, x, y, button, shift):
    pass
def onMouseUp(self, x, y, button, shift):
    pass
def onMouseUpMap(self, x, y, button, shift):
    pass
def onMouseMove(self, x, y, button, shift):
    pass
def onMouseMoveMap(self, x, y, button, shift):
    pass
def onDblClick(self):
    pass
def onKeyDown(self, keycode, shift):
    pass
def onKeyUp(self, keycode, shift):
    pass
def deactivate(self):
    pass
def onCircle(self, circle_geometry):
    pass
def onLine(self, line_geometry):
    pass
def onRectangle(self, rectangle_geometry):
    pass


texts: OpenLayers 3 with geoJSON - struggling with projections

texts:
I am trying to implement support for multiple projections in qgis2web. I believe that the correct approach would be to export the QGIS layers to EPSG:4326 geoJSON. I've implemented that.

However, I cannot get the points to display on the EPSG:3857 OL3 map. I've looked through the OL3 API, and I've only found that under ol.format.GeoJSON, the option defaultDataProjection exists. It's not clear to me whether I should set this to the JSON's projection 4326 or the map's 3857. I've tried both, but the points still do not appear.

A sample qgis2web OL3 map is here: https://www.dropbox.com/s/nyjvn8s0w14a04j/ol3map.zip?dl=0. Can anyone help me figure this out? What am I missing?

UPDATE: I might have figured this out. In the dropboxed code /layers/layers.js line 5, adding {dataProjection: 'EPSG:4326', featureProjection: 'EPSG:3857'} as the second argument to readFeatures might have done the job. Will confirm or deny after testing.


texts: clipping an area out of multiple satellite images in QGIS

texts: How can I clip the exact same area out of  multiple satellite images ( multiple bands)? I am working on a remote sensing project and only need a certain island but the images show much more.


texts: PostGIS functions do not exist even when public is in the search path

texts:
I have a postgres database with the postgis extension enabled in the public schema.  We have put the public schema in the search path, but still can only get the postgis function (i.e. st_x() ) if we explicitly refer to the public schema, but we get:

  select st_x(obspoint_geom) from sc2.observation_points;
                ^ HINT:  No function matches the given name and argument types. You might need to add explicit type casts.
  ********** Error **********

  ERROR: function st_x(public.geometry) does not exist SQL state: 42883
  Hint: No function matches the given name and argument types. You might
  need to add explicit type casts. Character: 8

Any idea what is happening?


texts: Is there a way to see all of the different values in each field?

texts: I'm using QGIS.  I have an attribute table that lists different land classifications under the "Land Type" column (agricultural, industrial, etc.).  I don't know what all the land types are and I would just like to find a list of them instead of scrolling down through thousands of entries. Is there a way to do this?  The metadata table is wrong so I have no documentation to go from.


texts: Showing "layer name" as text in ArcMap layout view?

texts:
I am trying to show the layer names of all visible layers in ArcGis 10.1
when creating a map in layout view.

First I tried to use "insert legend" and show there only the layername, but I can't disable the decription/symbology display.
That did not work, so I think it is better to use "insert dynamic text" instead.

That is easy for the Document name Document Name: <dyn type="document" property="name"/>
but I don't know, what to insert when I want the layernames.

What is the correct type and property?


texts: Using postgis to generate building shades

texts:
I am currently working on a solar potential estimation tool for existing buildings. The idea is to use buildings (=Polygon) shape, their height, and create the resulting shade at a given hour. I am only going to do the test for a few position of the sun. I know about the grass function r.sun.mask but it is overpowered for what I am trying to do.

For now I only want to get the shade when the sun is south, at 18° in the sky. The building should therefore cast a shadow about three times its height.

What I am trying to get:

I have been looking for a while but I did not find any tool in postgis to do this. I was thinking about buffering and then cutting the polygon but could not make it work. Do you know any tool or have any idea how to do it?


texts:
I'm using OpenLayers filter module for making WFS getFeature requests. Already used and, EqualTo and IsBetween, and all of them are working fine (Example below):

const featureRequest = new WFS().writeGetFeature({
  srsName: 'EPSG:25830',
  featurePrefix: 'CP',
  featureTypes: ['CadastralParcel'],
  outputFormat: 'application/json',
  propertyNames: properties,
  filter: equalToFilter('id', id, false)
});

I found that I need to add a filter so that the features returned must have a property value IN a range of values.

Example:

|--------|----------|---------|
|   id   |   geom   |  year   |
|--------|----------|---------|

Need a filter that gets the features with the year field IN 2000, 2010 and 2020

I have not specifically found the filter in the openlayers filter list.

  Does this filter exist? If not, how can I implement this using the
  existing filters?


texts: Carrying information from one shapefile to another based on location using ArcMap

texts:
I am a R user but I couldn't do this on R and decided it might be easier on ArcMap but I don't know how to use it that well.

I have two shapefiles for same county, one of them has 360 features but attached to a long distance of a road, other one 7345 features but short roads. Example; let's say there is a major artery speed limit range 30-55 in 10 miles. So they separated this road 25 pieces in one shapefile, other one has the volume of road for 10 miles there is only one feature.

Line1
class       : SpatialLinesDataFrame
features    : 7162
extent      : -84.70612, -83.98041, 30.27399, 30.68428  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
variables   : 52

Line2
class       : SpatialLinesDataFrame
features    : 376
extent      : 143516.4, 213981, 3353367, 3399153  (xmin, xmax, ymin, ymax)
crs         : +proj=utm +zone=17 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
variables   : 22

LineMerged
class       : SpatialLinesDataFrame
features    : 7538
extent      : 143511.3, 214019.4, 3353367, 3399153  (xmin, xmax, ymin, ymax)
crs         : +proj=utm +zone=17 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
variables   : 75

I want to combine these two files, they have different proj4.

I merged them on ArcMap, I have 7538 features and uses Line2. but volume is only available 376 location when I read it on R.

Is there a way carry information on Line2 to Line1 since it has more features and speed changes on a lot of locations but volume is same. if i can add volume in that location, would be amazing.

Getting information from one shapefile to another using QGIS is not helping me.


texts: QField hide / unhide features which are not created on the current date

texts: During field surveys I use QField to enter sightings of animals. Each sighting is stored with a date. When I visit the same site the next time I would like to be able to turn on / off the display of previous sightings so I only see the sighting of the current date. I tried applying a filter using the current date but did not succeed. How can I realize this?


texts: Download an Image Collection with a button using GeoJSON

texts:
I'm working on a GEE App where users can query for NDVI data of any country or given point, with an specific temporality. I also want the user to be able to download the result map.
I know this can be done with Export.image.toDrive() but this involves clicking on the Task tab, which isn't visible in a GEE App.
The closest I've been to solve my problem it's through a script I found for another question (I lost the link but I'll link the script with the license below). In this script, they use GeoJSON to encapsule the map view, but in the example they only use an Image from a Landsat8 collection, and I'm working with a lot more data.
How can I solve my problem?
Link to the example: https://code.earthengine.google.com/83af250cb4f68d038417cde28f538476
Link to my project: https://code.earthengine.google.com/20d6255337c2cba28112bb63e31a2888


texts:
I want to present data from a csv-file in a fact box next to a map. The map is created using QGIS (2.12.2) print composer. The data in the csv-file is updated regularly and I want to be able to link to the contents of it
by using an expression so that the data within can be presented as part of a text or in a table using a html rendering.

The file looks more or less like this:

zone,area
2,901403
1,109706
n,5061049
V,587821

I want to do this for a simple map. How can it be done?


texts: CartoDB static maps API and named placeholders

texts:
Is it possible to use static maps API using named maps placeholders?

I've a private dataset and I need to enhance some features at the output.

In more details, I've a layer of land uses and I need to create a PDF of an specific plot of my layer. I'd like to get an image of the layer with the selected plot enhanced over the others.

If it's not supported I'll try the hack using phantomjs + cartodbjs


texts: Connection to ArcGIS Server via QGIS fails

texts:
I am trying to load this map that comes from an ESRI Server in QGIS.
I tried two different ways:

DATA SOURCE MANAGER / PROTOCOL HTTPS or FILE and then adding the URL. When I do this it fails and appears this message. (INVALID DATA SOURCE
DATA SOURCE MANAGER / ARCGIS Map SERVER. I add the server connection using the URL. It loads the layer but appears a warning message next to it saying.

Network error: Error transferring https://egis.fire.ca.gov/arcgis/rest/services/FRAP/ownership/FeatureServer/export?bbox=-13849232.281800,3833711.951200,-12704981.773600,5162403.249100&size=1319,1531&format=jpg&layers=show:0&transparent=true&f=image - server replied: Bad Request
https://www.arcgis.com/home/item.html?id=f73858e200634ca888b19ca8c78e3aed .Maybe I am doing it the wrong way, but everytime I load it, I always get the same error.
Does anyone have an idea how to solve it?


texts: How to do unit conversion on a Terrain dataset?

texts: I was wondering if there is a way to change the elevation units of a Terrain dataset (from example from meter to feet)?


texts:
I wish to loop function that clips a shapefile by a feature class. I want to make the clip the same as the input name. The names of the files are WVHGT_201403281200_f180_201404050000.shp and I just want f180 part as the output feature class. The clip feature will remain consistent throughout.

I have the syntax down from here. I am just not sure how the write the for loop to do this.


texts: Import a shapefile to postgis with ogr2ogr gives: Unable to open datasource

texts:
I want to use ogr2ogr in order to import a shapefile in a postgis database.
I have successfully installed ogr2ogr and I run from the pgsql the following command:

ogr2ogr -f "PostgreSQL" PG:"host=localhost user=user_1 password=***** dbname=imports" world_boundaries.shp

What I get back is an error message:

Unable to open datasource `world_boundaries.shp' with the following drivers: --a list of drivers follows (ESRI Shape File etc.)

I have also tried to define the full path of the shapefile but I got the same message.

Also I tried to run the:

ogrinfo world_boundaries.shp

Same thing.

After fixing the issues with the permissions of the file I get the following error:

  ERROR 1: AddGeometryColumn failed for layer world_boundaries, layer creation has failed.
  ERROR 1: Terminating translation prematurely after failed
  translation of layer world_boundaries (use -skipfailures to skip errors)

I also tried to import it through the GUI shp2pgsql and I get the following error:

ALTER TABLE "public".""
Failed in pgui_exec(): ERROR:  permission denied for relation spatial_ref_sys
 CONTEXT:  SQL statement "SELECT SRID         FROM spatial_ref_sys WHERE SRID = new_srid"
PL/pgSQL function addgeometrycolumn(character varying,character   varying,character varying,character varying,integer,character varying,integer,boolean) line 50 at SQL statement
 SQL statement "SELECT AddGeometryColumn('',$1,$2,$3,$4,$5,$6,$7)"
 PL/pgSQL function addgeometrycolumn(character varying,character  varying,character varying,integer,character varying,integer,boolean) line 5 at     SQL statement

Shapefile import failed.

The issue this time was that this database user didn't have sufficient permissions. This fixed it:

GRANT ALL ON TABLE spatial_ref_sys TO my_user_name;

Next error message is:

  Warning 1: Geometry to be inserted is of type 3D Multi Polygon, whereas the layer geometry type is 3D Polygon.
 Insertion is likely to fail
 ERROR 1: INSERT command for new feature failed.
 ERROR:  Geometry type (MultiPolygon) does not match column type (Polygon)

So it seems I need to use the parameter: -nlt MULTIPOLYGON
But when I do so I get another error, which doesn't make any sense to me:

ERROR 1: PostgreSQL driver doesn't currently support database creation.
Please create database with the `createdb' command.
PostgreSQL driver failed to create PG:host=localhost user=my_user_name dbname=my_database password=password -nlt POLYGON

But it gets loaded using the shp2pgsql GUI.

Comment of @elrobis enabled this to finally work. Data loaded in db correctly!

ogr2ogr -f "PostgreSQL" PG:"host=localhost user=geonode dbname=geonode_imports password=geonode" -nlt GEOMETRY wld_bnd_adm0_gaul_2015.shp


texts: Covering polygon with fixed size squares?

texts:
Can you please advice me of an algorithm (which I can implement in a Java) to cover some polygon (not exact cover) with squares of certain size (aXa squares)?

The cover doesn't have to be exact - that is, the cover can be more than the original polygon but the squares can't be intersecting and every square has to have at least one point of intersection with the original polygon.


texts: Uniformly distribute spatial points over a given extent R

texts: Within a shapefile I am given an integer attribute that indicates the number of points that will need to be evenly (uniformly) distributed over the associated polygons extent. These spatial points will need to be stored in a separate shapefile.


texts: How to set appropriate spatial index for a feature class in Oracle?

texts:
I wanted to ask is there any known method how to set spatial index for my feature class in Oracle. I am working with ESRI 10.1 and I have various feature classes covering the whole area of my country (Poland). Those feature classes represent forest, rivers, roads, cities etc. They are of course represented by different geometires so I have polygons, lines, points and annotations.

After loading all data into Oracle I wanted to recalcute spatial index to improve the speed in which the features are drawn on the map. Of course I could use the recalculate function in ArcCatalog for each feature class but I wanted to do it on my own. Maybe it will give better results?

I have found interesting solution at ESRI website, however I would like to ask how did other user deal with this problem:

http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/Guidelines_for_choosing_a_spatial_index_grid_size/006z0000002p000000/


texts: OpenLayers, symbol every 100 pixel on line symbolizer

texts:
First: the above banner links this to How to use a sld file with Openlayers using OSM data? , but look update2.

I would like to add rail-road hatching to my OpenLayers map with a Symbolizer. Example how that is done in GeoServer is at [1] and there I could do that like this:

   <FeatureTypeStyle>
     <Rule>
       <LineSymbolizer>
         <Stroke>
           <CssParameter name="stroke">#333333</CssParameter>
           <CssParameter name="stroke-width">3</CssParameter>
         </Stroke>
       </LineSymbolizer>
       <LineSymbolizer>
         <Stroke>
           <GraphicStroke>
             <Graphic>
               <Mark>
                 <WellKnownName>shape://vertline</WellKnownName>
                 <Stroke>
                   <CssParameter name="stroke">#333333</CssParameter>
                   <CssParameter name="stroke-width">1</CssParameter>
                 </Stroke>
               </Mark>
               <Size>12</Size>
             </Graphic>
           </GraphicStroke>
         </Stroke>
       </LineSymbolizer>
     </Rule>
   </FeatureTypeStyle>

(This is only example directly from link, not exactly making the desired thing, but the shower of idea only)

Picture about what I get with that example code below.

Now, after that on screen, I can go further and get what I want, but the issue is the lines are a pure line VectorFeature coming to my code and I have to do the same thing with OpenLayers.

Question:

How I could come up with vertical lines, like done in GeoServer, now with OpenLayers Symbolizer? How to substitute lines with Z and get an electricity line styling?

My ideas so far:

I accomplished a two line road with black and white styles combined above each other and thought, that having a dash that is very long but has a little cap with white and very broad width and then equal black line with same width, there should come a vertical black line every certain period. My issue is to substitute though the lines with Z symbols and to that I run out of ideas. Playing with lines comes to end here.

My issue is to know how to get a symbol with certain cap to line. Then is everything cool.

Update 1:

Looking up for the docs [2] is not so promising, but maybe a custom style still is possible.

Update 2:

The question linked to here [3] answers only partially the addressed question. At least shape symbol on given code won't work with commands told in link [2].

Changes to OL renderer is needed, so that when parsing sld to svg it "handles" the new stuff properly. How to do that? Well, maybe I post answer later, if I don't get any hints before I invent it myself.

Links together:

[1] http://docs.geoserver.org/stable/en/user/styling/sld-cookbook/lines.html#railroad-hatching

[2] http://docs.openlayers.org/library/feature_styling.html

[3] How to use a sld file with Openlayers using OSM data?


texts:
So apparently the MAC install needs various stuff besides QGIS itself. The QGIS website sends you to an outside source for GDAL and for even the QGIS download itself. So I got a disk image from that source that has some stuff in it although there are no instructions. There is no "installer" that will boot in what was delivered and there are readme's that say you may need different python.

How do I parse how to install this and what additional resources I need?


texts: Reversing "combine" function on raster using R?

texts:
I have a raster that I think was "combined" using an ArcGIS function.

I am wondering if there is a way to extract just one raster layer from the combined dataset in order to use it separately, using R.

If you want to see the data it is available here

It takes this form in R:

class       : RasterLayer

dimensions  : 10995, 10349, 113787255  (nrow, ncol, ncell)

resolution  : 100, 100  (x, y)

extent      : -446880, 588020, -612570, 486930  (xmin, xmax, ymin, ymax)

coord. ref. : +proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD27 +units=m +no_defs +ellps=clrk66 +nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat

data source : G:\Faculty\Mann\Share\Fire\FRAP_FireThreat\FireThreat.tif

names       : FireThreat

values      : 0, 255  (min, max)

attributes
     fields : ID Rowid   COUNT FUEL_RANK FROTCLASS URBAN THREAT THREAT2PEOPLE
        min :  1     0       1        -1         0     0     -1            -1
        max : 97    96 5360409         3         3     1      4           424