texts: How to Create shapefile using Geometry field

texts:
I am working on ArcObjects. I want to create shapefile. I am very close to it but need some better alternative as I feel my solution may fall short in time when there are thousand of records.

I have a geometry field in a database. From that I get the string which is Linstring, Multilinestring or Point

Eg.

  "LINESTRING (-98.71702604799998 37.746413141000062,
  -98.717210399999942 37.738209632000064, -98.717347586999949 37.732091870000033)"

From this I manually extract the coordinates like

string d = st.Substring(12); // to get rid of LINESTRING label
d = d.Substring(0, d.Length - 1); // to get points as x1 y1,x2 y2,x3 y3
arr = d.Split(','); // now to get array [x1 y1, x2 y2, x3 y3]

Here is my code:

private IPolyline CreatePolyline(DataSet ds)
{
    ESRI.ArcGIS.Geometry.IPolyline pSegPoly1 = new ESRI.ArcGIS.Geometry.PolylineClass();

    pSegPoly1.SpatialReference = createSpatialReference();
    ESRI.ArcGIS.Geometry.IGeometryCollection pGeoColl12 = pSegPoly1 as ESRI.ArcGIS.Geometry.IGeometryCollection;

    foreach (DataRow row in ds.Tables[0].Rows)
    {
        string st = row["ShapeGeoWKT"].ToString();
        string[] arr = null;
        if (st.StartsWith("MULTI"))
        {
            string d = st.Substring(18);
            d = d.Substring(0, d.Length - 2);
            arr = d.Split(',');
        }
        else if (st.StartsWith("LINE"))
        {
            string d = st.Substring(12);
            d = d.Substring(0, d.Length - 1);
            arr = d.Split(',');
        }
        ESRI.ArcGIS.Geometry.IPoint[] pntArrayTemp = new ESRI.ArcGIS.Geometry.IPoint[arr.Length];
        int k = 0;
        ESRI.ArcGIS.Geometry.ILine[] linearray = new ESRI.ArcGIS.Geometry.LineClass[pntArrayTemp.Length];
        ESRI.ArcGIS.Geometry.ISegmentCollection[] paths = new ESRI.ArcGIS.Geometry.PathClass[pntArrayTemp.Length];
        foreach (string st1 in arr)
        {
            linearray[k] = new ESRI.ArcGIS.Geometry.LineClass();
            paths[k] = new ESRI.ArcGIS.Geometry.PathClass();
            pntArrayTemp[k] = new ESRI.ArcGIS.Geometry.PointClass();
            string[] temp1 = st1.Trim().Split(' ');
            pntArrayTemp[k].X = Convert.ToDouble(temp1[0].Replace(")", "").Replace("(", "").Trim());
            pntArrayTemp[k].Y = Convert.ToDouble(temp1[1].Replace(")", "").Replace("(", "").Trim());
            object obj1 = Type.Missing;
            paths[k].AddSegment((ISegment)linearray[k], ref obj1, ref obj1);
            k++;

        }
        for (int j = 0; j < paths.Length - 1; j++)
        {
            object obj12 = Type.Missing;
            pGeoColl12.AddGeometry((ESRI.ArcGIS.Geometry.IGeometry)paths[j], ref obj12, ref obj12);
            linearray[j].PutCoords(pntArrayTemp[j], pntArrayTemp[j + 1]);
        }

    }
    pGeoColl12.GeometriesChanged();
    return pGeoColl12 as IPolyline;
}

From this I get an array of the co-ordinates which I use to create shape for line. It is a very tedious procedure for many records (100k or more).

Is their any simple easy method in ArcObjects that I can use to get the coordinates from LINESTRING on the fly or I can use the LINESTRING to draw lines?


texts:
I want to generate tiles in custom size but can't find a tiler that can do it.
In gdal2tiles I cant find the size option.
TileMapper costs too much.

Is there any opensource utils that can generate tiles in custom sizes?


texts: How to loop the kernel for covering all pixels from a raster in 5x5 convolution window by using OpenCV?

texts:
Here, I have python code which is helping me in selecting 5 by 5 window from upper left corner of a raster (.tif). I want a process to automate through entire image in 5 by 5 convolution window.

import cv2
import numpy as np
from matplotlib import pyplot as plt
img = cv2.imread('class.tif')
print(img)
kernel = np.ones((5,5),np.float32)
print(kernel)

How to loop the kernel which will cover all pixels from a raster in 5 by 5 convolution window for getting all pixel values in this convolution window? I don't want to generate another raster. I want location of abnormal values within each window after statistical tests. How is it possible?


texts:
I'm trying to connect my QGIS to a WCS service using OAuth2 authentication (see curl example below to understand the authentication workflow).

the problem

The option "add WCS layer" seemed to be the obvious right choice, but the OAuth2 authentication makes it difficult. I would need to

get an access token
add this token to the http header along with the actual wcs request

Now. When I try to connect using "add wcs layer", I can choose to enter a username and password for basic auth, but I just cannot find any possibility to use the OAuth2 authentication protocol.

curl example

Using curl, the whole thing works as follows:

curl -i -H "Content-Type: application/x-www-form-urlencoded" -X POST -d "grant_type=client_credentials&scope=<myTestScope>&client_id=<myClientId>&client_secret=<myClientSecret>" <https-urlToGenerateToken>

I get as a repsonse my token in json format:

{
  "token_type":"Bearer",
  "access_token": <myToken>,
  "expires_in":3600,
  "scope": <myTestScope>
}

This, I can then use to call the service:

curl -i -H "Authorization: Bearer <myToken>" <wcsURL>/getCapabilities

the question

Is there a standard way to use this OAuth2 authentication workflow in qgis when connecting to a service, and I just couldn't find it? If there is no standard way, can I somehow "hook in" and add custom headers?

Edit - Tried GDAL driver:

As suggested in the comments section, I tried to use the GDAL driver. I didn't quite figure out how to use it from within QGIS (but I assume it's possible). Instead, I tried it from the command line (creating the respective xml-file and calling gdalinfo file.xml).

Could connect to the test server (no authentication)
Could connect to the server when using Basic Auth, by adding the <UserPwd> tag as documented in the GDAL Driver doc.
The problem with OAuth2 is still not solved. Setting <HttpAuth> to ANY and stuff doesn't work. Also, this makes total sense, as I would really need to make 2 different calls (1. token, 2. wcs-request) to two different url's, but the GDAL Driver accepts only one ServiceURL.


texts: Can GDAL report if a .tiff is tiled or not?

texts:
A certain piece of software that I use doesn't support .tiff DEM's that are tiled. Is there a way with a GDAL utility to determine if a .tiff DEM is tiled or not?

gdalinfo dem.tif

will report the compression type and the interleave type (band/pixel), but as far as I can tell, the interleave type doesn't have anything to do with tiling.

you can translate the tiff DEM with a creation option flag to make an untiled DEM, but that requires a translation and not just a report.

gdal_translate -co TILED=NO dem.tif dem_notile.tif


texts: Changing the time to GMT+05:30 in twitter tool in cartodb

texts: I am much more interested in exploring geotagged twitter maps using cartodb. By default, We will be giving a search term as well as the From/To (Date with time) usually it will be 30days. In this From/To, time will be in GMT+0 and my time is in GMT+05:30. I would like to know how to change the default GMT+0 settings in this twitter tool?


texts: Keeping attribute table open during QGIS close/re-open?

texts: When saving and closing QGIS with several open attribute tables, QGIS does not open these when opening the project next time. This is very annoying. Is there a workaround for this?


texts: Are WGS84 LatLon coordinates and EGM96 LatLon coordinates the same and difference for [LatLonEle] is elevation only?

texts:
Does WGS84 LatLon coordinates and EGM96 LatLon coordinates are the same?

Is difference for [Lat;Lon;Ele] is in elevation only?

If I have, for example, a [Lat;Lon] point in WGS84, I can take elevation from NASADEM_HGT.001 dataset (which is in EGM96) without any conversion of Lat\Lon and this will be accurate elevation for this point in EGM96. By adding calculated undulation for [Lat;Lon], I can find WGS84 Elevation.

Am I correct ?

How can I 'proof' that to myself? Any docs?


texts: Why can't I see any string fields in graduated renderer?

texts:
After the join I can see all my data in the attributes of the layer but when I want to select statistics in the graduated style to categorize, I don't see them. I see that the real type variables and not those of type string. Are there any parameters to be changed?

original:

après avoir fait la jointure je vois bien l'ensemble de mes données dans les attributs de la carte mais quand je veux sélectionner les statistiques dans l'onglet style pour catégoriser, je ne les vois pas. Je ne vois que les variables de type real et pas celles de type string. Y a t-il un paramètre à modifier?


texts: Creating heat map with fixed max color value in QGIS?

texts:
I am trying to create a heat map using raster data in QGIS 3.0. I have made them previously in ArcGIS pretty easily, but I am having difficulties in QGIS. The difficulty I am having is setting a fixed max color value (in this case, it's $400M). The procedure I take is as follows:

Load the raster set into QGIS
Convert to a heat map via Symbology (no changes made to default options)
Open Processing Toolbox, Interpolation, Heatmap (Kernel Density Estimation)
Save as .tif file.
Load in .tif file, go to Style
Change max value to the value I am interested in ($400M).

In the end, I just get a big circle looking like this.

In a nutshell, I want something like this (the max value is about 200M)

But instead, the best I can do is this, before I export out the heatmap as a .tif (after step 2) (where 200M is the max color, instead of 400M)

Any ideas here?


texts: Definitive list of Proj4 supported coordinate systems

texts:
I need to determine the minimum set of coordinate systems supported by all tiers of my project's stack. I have a list of EPSG codes supported by Oracle and GeoServer but I'm having trouble finding a list of systems supported by Proj4 (I intend to use the latest release of Proj4JS with OpenLayers).

The Proj4 site provides this page that "is intended to eventually contain information on all the projections supported by PROJ.4" but currently doesn't.

Does anyone know where I might find this information?


texts: LAStools (3D Analyst) scale factor application

texts: I'm trying to use the las2las (transform) module from the LASTools toolset to apply a scale factor to an LAS dataset. The problem is that the client wants us to apply the scale factor from a known XY coordinate. Does anyone know of a good workflow to incorporate an origin point in LAStools?


texts:
At the moment I am clipping osm.pbf files. From the osm file, I only want to extract the buildings. The outcome should be shapefile with all the buildings within a nuts2 area.

The following code works:

pbf = 'A://CORINE//OSM//%s-latest.osm.pbf' % (country)
nuts_wgs = "A:\\CORINE\\%s\\NUTS2_SHAPE\\%s_wgs.shp" % (country,nuts)
wgs = 'A://CORINE//%s//NUTS2_BUILDING//%s_buildings_wgs.shp' % (country,nuts)

os.system('ogr2ogr -progress -f "ESRI shapefile" ' + wgs + ' ' + pbf + ' -clipsrc ' + nuts_wgs + ' -sql "select building from multipolygons where building is not null" -nlt POLYGON -skipfailures')

But this code is relative slow. I have found some documentation about clipsrcsql and clipsrclayer. If I add one of those to the code, it will be much faster (approx. 6 times faster), but a clip to the nuts2 file does not occur anymore (it simply gives all the buildings of the country as an output).

Does anyone have an idea either on:

1) how clipsrcsql or clipsrclayer work (documentation is very limited)?

2) how to improve the speed of the calculation in general. Especially for areas with many buildings (e.g. Netherlands, France, Germany), a speed improvement is required?


texts: How to use wildcards for characters inside the name

texts:
I have shapefiles in that form

RCs4s3s2c10_S2C10
RCs4s3s2c20_S2C20
RCs4s3s2c30_S2C30
RCs4s3s2c40_S2C40
RCs4s3s2c50_S2C50
RCs4s3s2c60_S2C60
RCs4s3s2cb_S2CB
RCs4s3s2c10_1_S2C10
RCs4s3s2c20_1_S2C20
RCs4s3s2c30_1_S2C30
RCs4s3s2c40_1_S2C40
RCs4s3s2c50_1_S2C50
RCs4s3s2c60_1_S2C60
RCs4s3s2cb_1_S2CB
RCs4s3s2c10_2_S2C10
RCs4s3s2c10_2_S2C10
RCs4s3s2c20_2_S2C20
RCs4s3s2c30_2_S2C30
RCs4s3s2c40_2_S2C40
RCs4s3s2c50_2_S2C50
RCs4s3s2c60_2_S2C60
RCs4s3s2c10_2_S2C10
RCs4s3s2cb_2_S2CB

How can i use a second wildcard parameter in my code for the first group of fc which don't have a number between _ cause with the following code

featureclasses = arcpy.ListFeatureClasses("RCs*", "Polyline", "")
    for fc in featureclasses:
          ....

i get everything but i want only the first group of feature classes


texts: How to Call ArcHydro Tools in Stand Alone Python Script

texts:
I am currently using ArcGIS ArcHydro tools to perform 'Catchment Polygon Processing'.

However, it is proving time consuming to perform this for multiple rasters.
How might a write a python script to import ArcHydro tools?


texts: How to generate a report in Arcmap 10.2 which counts the number of entries in an attribute table of a certain value

texts:
I am generating an attribute table which may have many thousands of entries. There is a column (one of perhaps 20 or so) for 'vegetation type'. I need to generate a report (or some other easily exportable format) which automatically counts and lists how many times each vegetation type is included in the table. There's far too many entries to do manually and I will need to run the process a number of times so I'm hoping there's a way to automate the process.

For example the attribute table may look like
FID  COLOR
0     Red
1     Blue
2     Green
3     blue
4     blue
5     green
6     blue
7     green
8     red
9     blue
10    yellow

The report would then say something like
Blue: 5
Red: 2
Green: 3
Yellow: 1
Total: 11


texts: Calculating percentage with comparing two string fields

texts:
i have data like below

i need output like below


texts: Seeking equation/algorithm for angle between straight line and grid north?

texts:
I have a cylindrical equidistant projection map containing climate data(equal angle and equal distances). On this projection I have a straight line connecting two latitude and longitude points.

How do I calculate the angle between this line and the line of longitude of this grid ?

I am requesting the equation/algorithm for this problem because I do not know whether the cosine of the angle between the line vector and the longitude line is enough or whether the projection of the map has a role in this as well.


texts: Name attribute in a new field as a function of attribute in another field

texts:
I have table with a field containing attribute names in another language (oldname1, oldname2,...). I want to create a new field in which I can add (with field calculator) the translated name of each attribute in my language (newname1,newname1,...).

I thought the cursor would be the best way to realize this.
The original field is called fieldA and the new field fieldB

I tried with a cursor but my script isn't working. The error is given for line "if row in 'fieldA'= 'oldname1'

myTable = ":\\path to table"

cursor = arcpy.da.SearchCursor(myTable, ['fieldA', 'fieldB'])
for row in cursor:
    if row in 'fieldA'= "oldname"
    arcpy.CalculateField_management(myTable, 'fieldB',"newname1")
    elif...
    .....

I tried to add cursors for each field as variable before the loop:

fieldA= arcpy.da.SearchCursor(myTable, ['fieldA'])
fieldB= arcpy.da.SearchCursor(myTable, ['fieldB'])

But it is not helping. I also try to write ['fieldA'] instead of 'fieldA'

Can I create variables for each field and use it in the cursor? or do you have an idea of how I could fix this?


texts:
I'm following the http://workshop.pgrouting.org/ but my routing result not showed up at all, like the picture below:

my firebug also didn't give me an error, I'm using geoserver 2.4.2, postgresql 9.2 x64, postgis 2.1.1, and pgrouting 2.0.
Here are the code that originally came from the pgrouting latest workshop:

<!DOCTYPE html>

    ol3 pgRouting client

      #ol-map {
        width: 100%;
        height: 400px;
      }

      start
      final

<button id="clear">clear</button>

<script type="text/javascript">

  var map = new ol.Map({
    target: 'ol-map',
    renderer: ol.RendererHint.CANVAS,
    layers: [
      new ol.layer.TileLayer({
        source: new ol.source.OSM()
      })
    ],
    view: new ol.View2D({
      center: [-127998, 6974591],
      zoom: 10
    })
  });

  var params = {
    LAYERS: 'pgrouting:pgrouting',
    FORMAT: 'image/png'
  };

  var startPoint = new ol.Overlay({
    map: map,
    element: document.getElementById('start-point')
  });
  var finalPoint = new ol.Overlay({
    map: map,
    element: document.getElementById('final-point')
  });

  var transform = ol.proj.getTransform('EPSG:3857', 'EPSG:4326');

  map.on('click', function(event) {
    var coordinate = event.getCoordinate();
    if (startPoint.getPosition() == undefined) {
      // first click
      startPoint.setPosition(coordinate);
    } else if (finalPoint.getPosition() == undefined) {
      // second click
      finalPoint.setPosition(coordinate);

      // transform the coordinates from the map projection (EPSG:3857)
      // into the server projection (EPSG:4326)
      var startCoord = transform(startPoint.getPosition());
      var finalCoord = transform(finalPoint.getPosition());
      var viewparams = [
        'x1:' + startCoord[0], 'y1:' + startCoord[1],
        'x2:' + finalCoord[0], 'y2:' + finalCoord[1]
      ];
      params.viewparams = viewparams.join(';');

      // we now have the two points, create the result layer and add it to the map
      result = new ol.layer.ImageLayer({
        source: new ol.source.SingleImageWMS({
          url: 'http://localhost:9090/geoserver/pgrouting/wms',
          params: params
        })
      });
      map.addLayer(result);
    }
  });

  document.getElementById('clear').addEventListener('click', function(event) {
    // hide the overlays
    startPoint.setPosition(undefined);
    finalPoint.setPosition(undefined);

    // remove the result layer
    map.removeLayer(result);
  });

</script>

am i missing something? I've been followed the workshop and check it twice.
Please help me I've been through so much error whether it using the older tutorial sample or this is the latest tutorial.

best regards


texts: Python library for geometrical calculations using with GRASS GIS?

texts: I have to do some geometrical calculations in 2D(point, lines, line segments and intersection points) in Python, and then use it in GRASS GIS. What is the best library to use for it? I'm trying with SymPy, but I'm not sure is that good choice.


texts: QGIS print layout - map refresh not working

texts:
I would like to synchronize my QGIS map with the printing layout map.
According to the tutorial here:
https://docs.qgis.org/3.4/en/docs/training_manual/map_composer/map_composer.html
I can refresh my map by clicking the "Refresh" symbol.
Unfortunately, despite the command "Refreshing map" appearing at the very bottom when clicking the icon, the map remains exactly the same as it was embedded initially.

Is it possible to sync the QGIS map canvas with the map provided in the print layout?


texts: Trouble with accesing API code: "ReferenceError: esri is not defined"

texts:
I am developing a web map application in the ArcGIS API for Javascript.

I two have ArcGISDynamicMapServiceLayer layers that I am trying to run a simple query task on. The code I am using works in other applications doing the same task with different data.

My problem is, when I load the app (tried both Firefox and Chrome) in a browser, the map itself won't load (indefinite load screen). Console in both browsers shows this error:

  ReferenceError: esri is not defined.

Comment #13 on this Esri forum post alludes to the DOM issue, but I need help understanding how to handle the issue. Also, it's possible that the DOM thing is not the issue at all.

So, does anyone have experience with this error? How have you worked around it?

Thanks so much for your input.

UPDATE:
This is the code that is actually throwing the error:

  var scalebar = new esri.dijit.Scalebar({
  map: map,
  attachTo:"bottom-left"
  });

Note that I have included the dojo require statement for that class:

  dojo.require("esri.dijit.Scalebar");

And if I simply comment out that function, the same undefined error comes with the next esri object class referenced (which also has its appropriate dojo require statement in the doc).


texts: Versions of Python & QGIS for new user to both to install?

texts:
I have used GIS (very basic functions), but now want to be able to customize to my needs.

My initial research leads me to believe that Python + QGIS is the way to go.

I understand that there are several plugins that will be useful along the way, but do not have enough knowledge to discern which ones and if they are all compatible with Python V3.

The install will be on Windows 7 machine.

Since I am learning Python and QGIS for the first time I want to use the newest versions possible-provided they (and their respective plug-ins) are compatible.

What versions of Python and QGIS should I, as a new user to both, install?


texts: Create categorized rendered from styles in QGIS 3.4.3

texts:
How does this tool work? I've been trying to use it but it needs xml file and I dont know how to create one and running the tool without inputting a style database returns a single colorised map
Is the option of saving a style in postgresql no longer available in QGIS 3.4.3?
Tool under Cartography toolset
Image showing: no option to save style as style db visible


texts:
I retrieve a specific feature from my WFS Layer with a GetFeature request:

request=GetFeature&cql_filter=myattr1=x&

Now, I want to filter for two attributes, thus, I make:

request=GetFeature&cql_filter=myattr1=x;myattr2=y&

However, what I get now is:

<ows:ExceptionText>
java.lang.IndexOutOfBoundsException: Index: 0, Size: 0
Index: 0, Size: 0
</ows:ExceptionText>

There must be a syntax error. Where?


texts: QGIS plugins for creating circular route on footpath network

texts:
I'm looking for a way to create random circular walks of a specified distance from a given point. The end result will be similar to https://routeshuffle.com/ and http://www.routeloops.com/, but along a right of way network with cost surfaces tailored to walking.

I'm able to create a path between two points easily using least cost path, and could create random points to imitate a circular route; however I'm hoping there's a tool out there that will enable me to generate this random circular path without requiring these intermediate points.


texts: TileMill Reference Layer Plugin doesn't work w/ non-MapBox TileJSON

texts:
Has anyone been able to get TileMill's reference layer plugin to work w/ a TileJSON definition?

In its description, it says it can work either with a mapbox map id or a tileJSON url.  I've tried it out successfully with map ids, but when I give it a url to my tileJSON here, nothing loads and the console outputs:

[tilemill] Client Error: Script error.
[tilemill]     at :0
[tilemill] Client Error: TypeError: 'undefined' is not an object (evaluating 'a.tiles')
[tilemill]     at http://127.0.0.1:20009/assets/tilemill/js/vendor.js:5164

Anyone had any luck with this?

Update: the tileJSON file must be hosted on a server that supports CORS.  That solved one problem, but still doesn't work w/ the plugin (same error output).


texts: Which STRM resolution Srtm2Osm.exe uses?

texts:
The piece of program Srtm2Osm uses SRTM elevation model from NASA to produce countour lines in OSM format.

Does it use the SRTM1 at 30m or SRTM3 at 90m?


texts: Upgraded Geoserver - How to "reset" / start over with Geowebcache as it is throwing exceptions

texts:
So I upgraded Geoserver 2.9 running in Jetty on Windows to Geoserver 2.14 running under Tomcat.  I stopped the old, set the GEOSERVER_DATA_DIR parameter in Web.xml in the WEB-INF folder to the location used by Geoserver 2.9.

All is working great except for Geowebcache.  When I select "Tile Layers" under "Tile Caching" in the UI, I get an exception.

org.apache.wicket.WicketRuntimeException: Error attaching this container for rendering: [OddEvenItem [Component id = 4]]
...
Caused by:
Caused by: java.lang.NullPointerException
    at org.geoserver.gwc.layer.StyleParameterFilter.setLayer(StyleParameterFilter.java:188)
    at org.geoserver.gwc.layer.TileLayerInfoUtil.checkAutomaticStyles(TileLayerInfoUtil.java:162)
    at org.geoserver.gwc.layer.TileLayerInfoUtil.checkAutomaticStyles(TileLayerInfoUtil.java:140)
    at org.geoserver.gwc.layer.GeoServerTileLayer.getPublishedInfo(GeoServerTileLayer.java:345)
    at org.geoserver.gwc.layer.GeoServerTileLayer.getLayerInfo(GeoServerTileLayer.java:326)
    at org.geoserver.gwc.web.GWCIconFactory.getSpecificLayerIcon(GWCIconFactory.java:54)
    at org.geoserver.gwc.web.layer.CachedLayerProvider$1.getPropertyValue(CachedLayerProvider.java:43)
    at org.geoserver.gwc.web.layer.CachedLayerProvider$1.getPropertyValue(CachedLayerProvider.java:37)
    at org.geoserver.gwc.web.layer.CachedLayersPage$1.getComponentForProperty(CachedLayersPage.java:86)
    at org.geoserver.web.wicket.GeoServerTablePanel$5.populateItem(GeoServerTablePanel.java:238)
    at org.apache.wicket.markup.html.list.ListView.onPopulate(ListView.java:523)
    at org.apache.wicket.markup.repeater.AbstractRepeater.onBeforeRender(AbstractRepeater.java:124)
    at org.apache.wicket.Component.internalBeforeRender(Component.java:950)
    at org.apache.wicket.Component.beforeRender(Component.java:1018)
    at org.apache.wicket.MarkupContainer.onBeforeRenderChildren(MarkupContainer.java:1836)
    ... 144 more

If anyone knows where to look to diagnose this, that would be great.

There is not much in the cache at this point in time, so I'd be happy with instructions on how to "start over" with the caching of layers.  Is this a straightforward thing to do?

Can I "clear" the GWC folder somehow, restart Geoserver, and manually turn caching back on for each layer that should be cached?


texts: Find closest polygon from point and get its attributes?

texts:
Given a set of points and a set of polygons, I am trying to find for each point what the closest polygon is and bring over the attributes from that polygon to the point (together with the actual distance). Ideally, I would end up with the original set of points with 2 more fields: 'distance' and 'polygon_id'.

I'm trying to do so with the following query (using ST_DWithin to limit the calculations)

CREATE TABLE test AS
SELECT t1.point_id, t1.geom, t2.polyg_id, (ST_Distance(t1.geom, t2.geom)) as distance
FROM points t1, polygs t2
WHERE ST_Dwithin(t1.geom, t2.geom, 100)
ORDER BY (ST_Distance(t1.geom, t2.geom)) ASC
LIMIT 1

This query is just returning one record though. The 'LIMIT 1' is acting on the whole result. How can I make that 'LIMIT 1' act on the results for each segment individually?


texts: Rotate feature around defined fixed point using ArcGIS 10.4

texts: I have a polygon feature and want to rotate it around defined (fixed point). There is a Rotate tool in Editor toolbar but I can't define fixed rotation point. I am using ArcGIS 10.4 and want to know how to do that using this software.


texts: Assessing accuracy of image classification using ArcGIS Desktop?

texts:
I have classified images in ArcGIS Desktop with the support vector machine classifier, and am now struggling with the accuracy assessment. I have only been able to find tutorials on how to do it manually by creating points manually and calculating the kappa index with excel... but I cannot find any detailed steps on how to 'create accuracy assessment points' or 'compute confusion matrix' through ArcGIS, (under spatial analyst tools).

Any pointers?


texts: How to create a locator through Open Street Map Data

texts: I need to make a locator for UK addresses and I know Open Street Map has very good data for the UK, I was wondering how I can convert the data for house addresses, which is in xml form, into a locator in Arcmap.


texts:
We are currently testing out OGC GeoServer implementing WFS 3.0 and very exited about it.
We usually use OpenLayers for displaying the data.
Very often we need further data visualizaions - basic diagrams (bar, line, pie) maybe a KPI.
Of course we can always request the wfs-data, convert it if necessary and feed it to a chart library.
But I am curious, if there are any frameworks that might help in querying, preparing and/or visualising wfs data. I stumbled upon this github project which enables you to automate the described process to a certain extend with configurations like this:
 title: 'Total Population',
 server: 'https://geoserver.siasar.org/geoserver',
 layer: 'siasar:communities',
 property: 'population',
 function: 'Sum'

I wonder if there are similar projects doing this. Maybe Server extensions that enable more detailed queries of and arithmetic operatios on the WFS.


texts: How to choose pixel value for NoData when clipping satellite imagery?

texts:
I downloaded a Sentinel 2A image because I want to obtain different index for land usage. I need to clip every band of the satellite imagery using the Clip Multi Layers plug in. The plug-in defaults to a value of -999 for NoData cells.

What is the NoData pixel value I should use to clip the imagery so as not to interfere with pixel data.


texts: opening a georefrenced image in IDL in ENVI

texts:
I am working with idl in envi,I want to use ENVI_CONVERT_FILE_COORDINATES in case of finding pixel value of specified point,for this porpuse I have lat and long data of my point but I need to open image in IDL as a georefrenced image
I noticed that image that is opened using ENVI_OPEN_FILE is not georefrenced,But if I open the same image using ENVI GUI it is georefrenced. So I think the problem is with opening image with its georefrenced data.

How do I open an georeferenced image in IDL?


texts:
Suppose I have a dataset where each record has an x,y and r (the centroid and radius of a circle).  How can I create a polygon feature class containing parametric circles?

Sorry to those who have a great open source solution, but I'm stuck using python in ArcGIS 10.


texts: Conserving Labels during Export to CAD in ArcGIS for Desktop?

texts:
I can export files to cad (i use ArcGIS 10.2), but the converted file does not involve labels,.

Then i tried to use another option, right click to layer name in Table of Contents then click Convert Labels to Annotation then export to cad but still same problem.

The result is a useless vector data without labels, so how can i export a layer to cad with labels?


texts: FME batch processing results filename "NULL"

texts:
I have a FME workspace and a workspace runner. The workspace reads several sources (PostGIS & Shapefile) but only one is batched. In other words, all shapefiles in one folders are all analysed against files from other sources, which remain the same on every run.

My problem is that the writer has a parameter "fme_basename" as the shapefile name value, but when I run the workspace runner, only the first shapefile gets the correct name (in this case L42213D) but the second is "NULL". Also all the rest of the batched files result as NULL so the second file is just overwritten until the process ends. This worked in my other workspace which was almost similar, but I can't get it working here...


texts: I'm writing Python code for ArcGIS using the arcpy library and now I'm experiencing a little problem with labels. In fact I'm not able to rotate them in layout via arcpy (I should rotate all labels of 90 degrees). I found a way to display correct field, round the number, but not to position the label correctly.


texts:
Georeferenced raster in the QGIS not filmed the final display by raster footprints.

Everything is in EPSG 32633.
Path in  GdalToools:

Red is original footprint for image.
I need rotate image to red area
Any idea to solve the problem?


texts: v.to.rast command in one shell script with other commands

texts:
I am using the latest GRASS Gis and QGis versions.
In order for a study project, I have to rasterise a vector layer using v.to.rast. This has to be done in my already existing shell script, in which there are various other commands, for instance "mkdir" in order to create folders or ogr2ogr commmands.
Unfortunately, these two commands dont work in my general shell script, only in the grass gis cmd, which opens by opening grass gis.
Is there some command missing, like for opening Grass Gis or creating a new location? If so, how exactly is the command called?
Unfortunately these two commands dont work:
v.import input="E:\Universitaet\WS 2020 - 2021\Seminar GIS\Gruppenarbeit\Dateien\\Vegetation\Irrelevant_data_deleted\vegetation_two_columns.dbf"
v.to.rast input="vegetation_two_columns.dbf" output="E:\Universitaet\WS 2020 - 2021\Seminar GIS\Gruppenarbeit\Dateien\Wildfireanalysis\Vegetation\Rasterised\Vegetation_Santa-Cruz_rasterised"
Thanks in advance!
Unfortunately, the GDAL command Rasterise (vector ro raster) is no option either, since it doesnt allow string variables.


texts: Extract all the polygon coordinates from a SpatialPolygonsDataframe?

texts:
I want to do something relatively simply but it seems to be very difficult to do, though I am very new to using R for spatial data so apologies if I've missed something obvious.

I have a list containing 161 shapefiles of bird distributions in the world, they are each in the format of a SpatialPolygonsDataFrame. I also have world climate data in a raster stack format for 8 variables. Both sets of files are in the same projection (lat and long WGS84). I want to create new columns for the climate data with 1's where the species polygon overlaps with the climate data and 0 where it doesn't. I tried various functions to do this, quite a few initial options didn't work, the best option appeared to be:

pnt.in.poly(XY coordinates of climate data points, XY coordinates of species polygon)

I changed my raster of climate data into points to get the XY coordinates. I unified the polygons for each species shapefile and then used:

species@polygons[[1]]@Polygons[[1]]@cords

to get the coordinates for my species polygons as a matrix. However, I discovered that:

species@polygons[[1]]@Polygons[[1]]@coords

wasn't actually getting all the coordinates of the polygons only some of them.

I tried using fortify from the ggplot package, and it gives me a dataframe with coordinates but some of them are for holes in the polygons. I ran that with the pnt.in.poly function excluding the holes but it took an hour to run and when I looked at the output maps after it matched for some of the range but not all of it- there were holes in the distribution where there shouldn't be holes.

I also tried this method that I found somewhere else:

lonlat <- lapply(slot(species, "polygons"), function(x) lapply(slot(x,"Polygons"), function(y) slot(y, "coords")))

which gives me a list with 436 items- but I'm still not sure how to get the coordinates into a matrix form that will work in pnt.in.poly.

I'm sure this should be easier to do! Does anyone know a better way?


texts: Measuring the sides of a polygon one by one

texts:
At the moment I am working on an OpenLayers project or rather I am learning on a project.
I stumbled upon a feature:
When you click on a polygon on the OpenLayers map and that plygon is selected and marked, I would like you to bring me from PostGis or calculate openlayers the size of the sides of the polygon one by one and pass them on that line in "m "..

geometry coordinates:
flatCoordinates: Array (12)
0: 3024661.117360038
1: 5657480.092961456
2: 3024466.963703795
3: 5656972.883250752
4: 3024488.1992827402
5: 5657122.118405258
6: 3024485.744880199
7: 5657251.372586279
8: 3024551.9375920547
9: 5657592.494035185
10: 3024661.117360038
11: 5657480.092961456
length: 12

I tried something after a long search on this forum and I found this:
feature.getGeometry (). getArea ()

but it only returns a number to me, not an array with all the sides of points A-B, B-C, C-D and so on.
Here I would ask for your help ...


texts: Cut a linear feature according to the directions or azimuths

texts:
I have a shapefile linear that contains a single items as in the image and I want to divide it into several section respecting the azimuth and the line direction.

i have this code :

    import os
    import fiona
    from shapely.geometry import mapping
    from shapely.geometry import Point,LineString
    import math

    os.chdir('SHP')

    def azimut(point1, point2):
           '''Retourne l'azimuth de la ligne entre 2 points shapely'''
           angle = math.atan2(point2.x - point1.x, point2.y - point1.y)
           return math.degrees(angle)
           print angle
           **if angle &amp;&gt;= 0: else math.degrees(angle) + 180** **Error syntax**

    **def paires(liste):
         '''parcourt une liste par paires'''
         for i in range(1, len(liste)):
             yield liste[i-1], liste[i]**

    with fiona.collection('BOM1.shp', 'r') as entree:
        # copie du schema de la couche et création d'un nouveau champ 'azimut'
        schema = entree.schema.copy()
        schema['properties']['azimut'] = 'int'
        # création d' une nouvelle couche avec le schéma résultant
        with fiona.collection('testligne_azim.shp', 'w', 'ESRI Shapefile', schema) as sortie:
            for ligne in entree:
            # utilisation de la fonction paire() pour extraire les segments de lignes
                for seg_start, seg_end in paires(ligne['geometry']['coordinates']):
                    # création d'une ligne en fonction des points des segments
                    line_start =Point(seg_start)
                    line_end = Point(seg_end)
                    segment = LineString([line_start.coords[0],line_end.coords[0]])
                    # copie des attributs d'entrée et ajout de la valeur résultante
                    # de la fonction azimut()
                    elem = {}
                    elem['properties'] = ligne['properties']
                    elem['properties']['azimut'] = azimut(line_start, line_end)
                    elem['geometry'] = mapping(segment)
                    sortie.write(elem)

but the problem is in :

if angle &amp;&gt;= 0: else math.degrees(angle) + 180

i have syntax error angle &>= 0

and in:

def paires(liste):
         '''parcourt une liste par paires'''
         for i in range(1, len(liste)):
             yield liste[i-1], liste[i]

i have this error position 29: ordinal not in range(128)


texts: Which value should I assign my digitalized features?

texts:
I have a scanned an old document that contains geochemical interpolated values (no way to work with the digital file). This is the original image:

And this is the digitalization I made from it:

Now, given that the original image is an interpolation of geochemical measures taken in field, which values should I give to each features of my polygon?

For example the crescent moon shape in the upper left corner has a maximum value of 100 (in the original file). Is the half of the original value a good idea (so 50 in this case)?

Someone has any suggestion?

I hope I made myself clear.

p.s. software used: QGIS 2.0


texts: How does Environments > Raster analysis > Cell size works in ArcGIS?

texts:
How does Environments > Raster analysis > Cell size works, when creating a raster from a polygon using ArcGIS?

I have tried to handle both "Maximum of Inputs" and "Minimum of Inputs" and nothing seems to change: the tool creates the same raster.


texts: How to use GDAL/python to fill the no data value

texts:
I am using GDAL tool to merge many geotiffs and get the following result:

how to fill the no data and change the no data value to 1?


texts: How to convert raster data into regularly sampled points?

texts: I am wondering if there's another way to convert raster data (ex: DEM) into regularly sampled points. I've used Raster to Point under ArcToolbox previously, but if the raster data is huge, that's a lot of points! Is there a way to specify the number of points or to specify the cell width and height so they're nicely spread out? I am aware of the Create a Fishnet tool, but that doesn't seem to be working for me at the moment. Not sure what I'm doing wrong.


texts:
I am trying to style a .json layer with a numbe of classes based on an attribute.

I and getting the error below which is apparently in ol3.js based on debugging in chrome:

 Uncaught TypeError: Cannot read property 'e' of undefined

I have tried to use ol-debug but then I get errors earlier on

The code I am using for the layer is

var lyr_RoFS = new ol.layer.Vector({
      source: new ol.source.GeoJSON({url: 'layers/RoFS.geojson'}),
      title: "Risk of Flooding",
      style: (function() {
          var defaultStyle = [new ol.style.Style({
             fill: new ol.style.Fill({color: 'navy'}),
             stroke: new ol.style.Stroke( {color: 'black', width: 1})
             })];
          var ruleStyle = [new ol.style.Style({
             fill: new ol.style.Fill({color: 'olive'}),
             stroke: new ol.style.Stroke({color: 'black', width: 1})
             })];
          return function(feature, resolution) {
              if (feature.get('prob_4band') === 'High') {
                 return ruleStyle;
              } else {
                 return defaultStyle;
              }
          };
      })


texts: Google Earth Engine: Calculate milliseconds since epoch from datetime of uploaded shapefile?

texts:
I want to get the milliseconds since epoch of thousands of points, each with different timestamps. Here is what I believe I need to do, but my Google Earth Engine knowledge is limited and I can't find a good way to do this:

Upload a shapefile with a column containing date/time (DONE)

Set the time/date string as a Date column (using .set()) and convert it to a Date object with ee.Date())

Calculate the "time" column using .millis() on the column created in step 2

My best effort so far has been this:
var times = points.set('timestamp', ee.Date(points.select('timestamp')).millis());

where:
(timestamp = the "date/time" column, formatted as dd/MM/YYYY HH:mm:ss)
(points = name of uploaded shapefile table, the FeatureCollection)
The console error says:
FeatureCollection (Error)
Date: Cannot interpret  as a Date.


texts: How to easily identify/classify coniferous forest from one date Landsat image?

texts:
I am pretty new in remote sensing and I am trying to identify/classify coniferous forest cover from single date Landsat scene. According to my preliminary web research I have these possibilities:

convert scene to NDVI values. Using modal values of NDVI histogram, I can separate scene pixels into forested and non-forested area
use modal value of band 2,3 and 5 (B2) to identify "forest peak" and class scene to forest/non-forest (Huang, 2008: Use of a dark object concept and support vector machines to automate forest cover change analysis). Other scenes characteristics (rocks, rivers) have to be removed using Tasseled cap brightness values

Do you know another simple approach to classify forest cover in mountainous area? I dont really want to apply maximum likelihood classification. Maybe it is better to used unsupervised classification?

I am using ERDAS, ArcGIS 10.2 and R


texts:
I'm starting to play with EE and I jumped right into the online web-based IDE for the Earth Engine JavaScript API (https://developers.google.com/earth-engine/playground).

I see there is also a Python API but it seems much more complex to play with.

Also saw somewhre (https://youtu.be/m1ejxSi3l8s?t=418) Node.js might be used.

What are the alternatives?

In the online code editor it is not possible i.e. to edit two files at same time, which might be a limitation overcome by other approaches. I guess setting up Python environment will allow me to do that, may be to use Sublime Text as an editor?


texts: How to solve "Error creating GEOS Coordinate Sequence" in Postgis?

texts:
When I run this query:

TRUNCATE sertaozinho.apps;

INSERT INTO sertaozinho.apps (id,tipo,the_geom)
VALUES (DEFAULT, 'app' ,(
    SELECT st_union(
    CASE
    WHEN area.cobertura = 'lago' THEN ST_Transform(st_buffer(area.the_geom,15), 32723)
    WHEN area.cobertura = 'rio' THEN ST_Transform(st_buffer(area.the_geom,100) , 32723)
    WHEN linhas.name = 'corrego' THEN ST_Transform(st_buffer(linhas.the_geom,30) , 32723)
    WHEN pontos.tipo = 'nascente' THEN ST_Transform(st_buffer(pontos.the_geom,50) , 32723)
    END )
    FROM sertaozinho.area, sertaozinho.linhas, sertaozinho.pontos
    )
    )

I get this output:

ERROR: Error creating GEOS Coordinate Sequence
Estado de SQL:XX000

The query is only a sequence of buffers and then an union. The dataset contains about 1000 polygons and 200 lines. I have already succeeded in running this query, but after some minor edits in the polygons it stopped working. Any clue?
EDIT: Sometimes after this query the server stops and needs to be restarted.


texts: Getting NDVI values into .csv format using Python

texts:
I have written the following code to get the NDVI values. I want this NDVI array in .csv format for data analysis.

import numpy as np
import rasterio
from rasterio.plot import show

b4 = rasterio.open('C:\\Users\\ARSH\\Desktop\\LC08_L1TP_143051_20200120_20200120_01_RT.tar\\LC08_L1TP_143051_20200120_20200120_01_RT_B4.tif')
b5 = rasterio.open('C:\\Users\\ARSH\\Desktop\\LC08_L1TP_143051_20200120_20200120_01_RT.tar\\LC08_L1TP_143051_20200120_20200120_01_RT_B5.tif')
red = b4.read(1).astype('float64')
nir = b5.read(1).astype('float64')
np.seterr(divide='ignore', invalid='ignore')

check = np.logical_or ( red > 0, nir > 0 )
ndvi = np.where ( check,  (nir - red ) / ( nir + red ), -999 )
ndvi

I am getting the below result

I also tried using pyexcel

import pyexcel
pyexcel.save_as(array=ndvi,dest_file_name="C:\\Users\\ARSH\\Desktop\\ndvir.csv",dest_delimiter=':')

but my .csv file is not taking values in a single column. All values are coming in a row.

How can I do this?


texts: Creating new Projected Coordinate System for the World in ArcMap?

texts:
I'm trying to make a completely new projection of the world in ArcMap. Many different projections are available, and I know they can be adapted by changing the parameters.

However, my goal is to create a completely new one (kind of the opposite of Robinson), as shown below. Does anyone know how to achieve this?


texts: Converting point cloud to raster with PDAL

texts:
I would like to convert my test.las file to a DTM file
I tried using this command:

(LIDAR) C:\Users\Dell\Documents\LIDAR>pdal pipeline test.json

But it's giving me as an error:
(PDAL Debug) Debugging...
PDAL: Pipeline: parse error at line 1, column 1: syntax error while parsing value
- invalid literal; last read: '■'

I tried to test with another file.json, but I'm having the same problem.
here is my test.json:
{
    "pipeline": [
        "C:/Users/Dell/Desktop/DATA/lidar.las",
        {
            "filename": "RasterTEST",
            "output_type": "mean",
            "gdaldriver": "GTiff",
            "radius": "10*1.41",
            "type": "writers.gdal",
            "resolution": "10"
        }
    ]
}


texts: QGIS field calculator in the graphic modeler: Is it possible to refer to a field name that contains a substring?

texts:
I am building a model in which I join attributes from a .csv file (let's name it "upload" here) to a layer.
I want to work with some of the new attributes in the next step of the model and I am having the problem that the field names in "upload" are always different except for the first part of the name.
So what I am looking for is a way of doing the following, e.g. in the field calculator:
Take the field of a certain layer that has a name beginning with "substring" and multiply its values by x.
Alternatively, I could imagine a solution where I refer to a field at a certain position (like the last field in a layer).
So far, I only managed to refer to fields of a completely known name.


texts: Describing Geodatabase Content With Python and XML Metadata

texts:
I am trying to create a descriptive database catalog by implementing the great solution described here. This user's solution creates a csv file that identifies a geodatabase's content, including elements from its metadata. In the GetMetadataItems function of this user's script, the abstract and purpose metadata items are accessed by using the elements as "idinfo/descript/abstract" and "idinfo/descript/purpose", as shown below in a portion of user's script.

def GetMetadataItems(dataset):
    """Retrieves certain metadata text elements from the specified dataset"""
    tree = GetMetadataElementTree(dataset)
    originator = GetElementText(tree, "idinfo/citation/citeinfo/origin") # Originator
    pocorg = GetFirstElementText(tree, ("idinfo/ptcontac/cntinfo/cntperp/cntorg", # Point of contact organization (person primary contact)
                                    "idinfo/ptcontac/cntinfo/cntorgp/cntorg"))      # Point of contact organization (organization primary contact)
    abstract = GetElementText(tree, "idinfo/descript/abstract") # Abstract
    purpose = GetElementText(tree, "idinfo/descript/purpose") # Purpose
    searchkeys = ListElementsText(tree, "dataIdInfo/searchKeys/keyword") # Search keywords
    themekeys = ListElementsText(tree, "idinfo/keywords/theme/themekey") # Theme keywords
    del tree
    metadataItems = (originator, pocorg, abstract, purpose, searchkeys, themekeys)
    return metadataItems

When I run this script, however, most of the values for the purpose and abstract columns are blank. I have tried to change these strings to read "ItemDescription/idCitation/idAbs" and "ItemDescription/idCitation/idPurp" as well as "idCitation/idAbs" and "idCitation/idPurp" My guess is that this tool is not working for me because of this portion of the script. How can I make sure the script is able to access these metadata items in my geodatabase?

Perhaps I need to  ensure the definitions in the GetMetadataItems function adhere to the specific metadata format of my geodatabase.


texts: OpenLayers2: show polygon of images in mobile

texts:
I'm using Phonegap 2.6 and OpenLayers 2.12 for my mobile application.

Feature:
in a panel (jQuery Mobile), I'm clicking in a checkbox and georeferenced images are displayed in the map that is already loaded. These images fits with that we call "blocks" in the map.

With Firefox:
usually, I'm working with Firefox at the beggining of developing stage. So the feature is working, the OpenLayers image is well created and displayed on the map.

With the device:
that doesn't work. Instead it draws a polygon around each image.

The method that provokes that is addLayers

Do you have any idea of what happens?
The only difference between the working and not working case, is that with the device my code is encapsulate with phonegap, but I have no clue why that should be a problem.


texts:
We have downloaded OSM data for our region and loaded them into PostGIS using osm2pgsql. So far so good. On top we have GeoServer and we would like to expose the OSM data using WMS. We would as well like to style data like in "standard" OpenStreetMaps (e.g. osm.org) as this cartography is well designed. We have tried to download som standard SLD for OSM data, but there seem to all sorts of problems using these SLDs (XML parsing errors, invalid content, etc.).

Any good ideas how to reuse an OSM SLD to style OSM through WMS?


texts:
I am trying to use Export Web Map Task but it has a limit of 1000 records, how can I use this service or replicate this service to execute my print command.

Any suggestions? I am working on ArcServer 10 not 10.1 and so cannot use the print task and print widget.


texts: "Flight recorder" for automobiles/motorcycles?

texts:
As I sat in traffic court last week for a non-moving violation, many of the cases were clearly he said/she said issues.. particularly with stop sign violations.   I can easily see how this specific problem could be a source of abuse by both police and citizens.

Case in point - I stop at all stop signs.  I usually count to 2 before moving, but because of the location of the white line, It may look like I never really stopped to someone coming along the cross street.

Are there any court-admissible devices that will accurately (and affordably) record vehicle movement, location, and perhaps speed?


texts: ArcGIS Field Calculator reformat Lat Long values

texts:
I have a table which is filled with point values that have been exported from an external source. Unfortunately the Lat Long values are formatted strangely and this is preventing me from simply displaying XY data to add the points to the map.

An example Latitude value is: 15 13 51S

I am wanting to create a field calculation in a new field which can convert that value into a decimal degree value so I can easily display the XY data.

My understanding is that I need to split the value into four parts (the degrees, minutes, seconds and direction values) and then divide the minutes value by 60, divide the seconds value by 3600, add those two values to the degree value and assign a negative value if the direction is S.

I have little experience using field calculator for complex calculations like this.

How can Python can be used to do it?


texts: How does Spatialite handle concurrent access?

texts: I understand that Spatialite is meant to be used by a single user/machine at a time.  However if the data is stored on server, there is no guarantee that someone on another machine will not also open the same database and use it.  Has anyone got any experience with how QGIS and spatialite react to being in a multiuser environment.  In most cases, the layer would be a read only layer.  QGIS allows multiple people to read the same shapefiles, which is a similar situation.  What happens if multiple people have concurrent editing sessions open?


texts:
I am getting very difficulties to add a map canvas widget through QtDesigner. I want a window with some buttons and a space to open the map canvas. How do I specify that? What widget did I choose?

Also I want to open the OpenLayers maps (Google Maps, Bing maps, etc) to overlap in the window with vector files.

Is there any widget which combine these two forms of view?


texts: v.rast.stats GRASS tool missing in QGIS 2.6.0 brighton

texts: I'm still relatively new to the QGIS world, and the GRASS toolbox. I am trying to do zonal statistics on a shapefile I have in which I want to know the max and min elevation from my DEM. I have used the QGIS zonal statistics tool, but it only returns a count, sum and mean elevation value. After doing some research the v.rast.stats GRASS tool seems ideal. It is however missing from my toolbox, and I have no idea how to go about getting it.


texts:
I've written a program that maps out a specific point file of interest. The program produces a map (clipped map attached below) but I can't figure out how to change the symbology to different colors. Is this possible? I've been looking through the arcpy.mapping library on ESRI and I haven't been able to find much on the subject.

if land_result == 1:
 with arcpy.da.SearchCursor(inFC, ['edited_nam', 'Field12']) as cursor:
  for x in cursor:
    print 'The %s %s contains tribal lands within its buffer zone' % (x[0], x[1])

  mxd = arcpy.mapping.MapDocument(r'P:/Python/topo_template.mxd')

  add_substation = arcpy.mapping.Layer(output_tribal_name)
  add_triballand = arcpy.mapping.Layer(triballands)

  df = arcpy.mapping.ListDataFrames(mxd, 'Data_Frame')[0]
  df.extent = add_substation.getExtent ()

  for df in arcpy.mapping.ListDataFrames(mxd):
    df.rotation = 0
    df.scale = 32000

  sub = arcpy.mapping.AddLayer(df, add_substation, 'AUTO_ARRANGE')
  wet = arcpy.mapping.AddLayer(df, add_triballand, 'AUTO_ARRANGE')

  north_arrow = arcpy.mapping.ListLayoutElements(mxd, 'MAPSURROUND_ELEMENT', 'North*')
  narrow = north_arrow[0]
  narrow.elementPositionX = 4.25
  narrow.elementPositionY = 1.25

  map_legend = arcpy.mapping.ListLayoutElements(mxd, 'LEGEND_ELEMENT', 'Map_L*')
  legend = map_legend[0]
  legend.elementPositionX = 0.5
  legend.elementPositionY = 1.33

  map_scale_bar = arcpy.mapping.ListLayoutElements(mxd, 'MAPSURROUND_ELEMENT', 'Scale*')
  scalebar = map_scale_bar[0]
  scalebar.elementPositionX = 4.25
  scalebar.elementPositionY = 0.33
  scalebar.elementWidth     = 3.75
  scalebar.elementHeight    = 0.3

  map_title = arcpy.mapping.ListLayoutElements(mxd, 'TEXT_ELEMENT', 'Map_T*')
  title = map_title[0]
  title.elementPositionX = 4.25
  title.elementPositionY = 10.75
  title.text = output + ' Wetland Buffer'

  map_data_frame = arcpy.mapping.ListLayoutElements(mxd, 'DATAFRAME_ELEMENT', 'Data_F*')
  dataframe = map_data_frame[0]
  dataframe.elementPositionX = 4.25
  dataframe.elementPositionY = 2.75
  dataframe.elementHeight = 7.5
  dataframe.elementWidth  = 7.5

  mpc_image = arcpy.mapping.ListLayoutElements(mxd, 'PICTURE_ELEMENT', 'MPC*')
  mpcimage = mpc_image[0]
  dataframe.elementpositionX = 8.0
  dataframe.elementpositionY = 1.55

  arcpy.mapping.ExportToPDF(mxd, r'P:/Python/' + output + 'topo_tlands.pdf')

  del mxd


texts:
I am trying to help out a friend of mine. She learnt GIS in school in the past using ArcGIS but has recently installed QGIS package on her Mac OSX.

When she tries to render raster into shapefile, she gets error that gdal_polygonize.py is not found.

Does it file come with a specific python package?

I could not find this exact file on the Internet. As a side note, she is able to import the gdal package in her Python.

Can anybody help?


texts: Stacking/Aligning High-Res UAV Mosaics from Different Dates in QGIS?

texts:
I have two mosaics taken in different band sets respectively, on two different flights- of a crop field. The flights were taken a day apart, so the vegetation cover has not changed. The only issue is the shadows, and slight distortions between the mosaics - as they were stitched from UAV imagery.

I tried georeferencing one to the other in QGIS 2.18; however, there is still slight offset which makes the creation of indices problematic.

How can I get a more accurate pixel to pixel alignment using QGIS?


texts: Calculating Viewshed of a region rather than a point

texts:
There are various calculations to find the viewshed of a point and all of them have different runtimes. the fastest is the calculation of line of sight angle and compare with upcoming cells.

But what about when we have a region rather than a point. Imagine a city which we want to compute its view on the surrounding terrain, no matter where the observer or an unlimited number of observers standing inside that city's boundary.

If we want to repeat the conventional viewshed for all region's cells, this would be a drastic computation.

Aside from This, is there any other method available? please provide links if possible.


texts: Arcmap: how to select features with attachments?

texts:
I have an attachment enabled feature class that has photos attached to points. However while any feature in the layer might have a photo, only 10-15% actually do. How can I narrow in on just those?

The only methods I've discovered rely on one by one interrogation (identify popup windows or attribute table while editing), which is not at all interesting to me in this context.


texts: Clip polygon with the extents of contour polylines

texts: I have a polygon shape that represents the boundary of Fennoscandinavia and a polyline shapefile representing two contour intervals for a trait. I would like to create a new polygon representing the portion of the Fennoscandinavia boundary polygon between the two polylines using ArcMap 10. There does not appear to be any tools in the 'data management', 'conversion' or 'Spatial Analyst' toolboxes that work. Any suggestions?


texts: Time Manager in QGIS works but does not display timestamp on map

texts:
I am using QGIS 3.12 to run the Time Manager plugin. The data is formatted correctly in Y-M-D form and Time Manager runs as intended (data values change with time); however, timestamp is not displayed in the bottom right corner despite "display frame start time on map" being clicked on.   I checked "time display options" and selected a color to contrast the map background.  I also tried with different basemaps, and without a basemap, and the timestamp is not displayed in all cases.

After much toggling around, selecting W view in "time display options" makes the timestamp barely visible off the frame. I have tried changing the map resolution and even the projection system, but the timestamp appears almost hidden to a corner. Changing the display option to any other locations, such as NW, E, W, makes it disappear again.


texts: Accuracy problems with 'Feature To Polygon' in ArcGIS 10.7?

texts:
I have a question regarding the accuracy of Feature to Polygon in ArcGIS 10.7. My data is available here.
I have a glacier shapefile and a polyline. I want to split the glacier using the line feature with Feature to Polygon. I get no error messages, however, the results seems a bit off as the cut appears more than 100m away from the actual line. I have repeated this with several other lines, using the GUI as well as arcpy - the error stays the same.
What am I missing?
Does this have something to do with the projection?
Both objects are in the same projected coordinate system (WGS_1984_UTM_Zone_43N).
EDIT:
I'm trying to do this in Python. I have access to the Advanced License, so I've tried CutPolygonByLines, but the code gives me errors. SplitPolygonsWithLines works, but only as a toolbox as the code provided is pretty unreadable.


texts: Is there any restriction on file size while using arcgis tool box or arcpy methods

texts: I am working with ArcToolBox's tools and their python supports. I found there is some restriction on file size of shape files/feature classes etc. It makes me interested to know about is there any restriction with file size on these tools/python methods also? As an example, suppose I am using "Con" method of "Spatial Analyst Tools". Does this method has any size restriction? I read "tool help" document. But I found nothing about file size here.


texts: Importing OPenStreetMap data into ArcSDE using osm2pgsql?

texts:
I am trying to import OpenStreetMap (osm) data into a PostgreSQL database using osm2pgsql, and then I need to convert that data so I can use it in ArcSDE/ArcGIS Desktop 10.2.

Here is what I do:

Using the PGAdmin GUI tool, I create a database call OSM
I then enabled the PostGIS on the database:

CREATE EXTENSION postgis;

Since I am on a Window environment, I created a path variable to the osm2pgsql.exe
Open a command prompt
Load a small sample set of OSM data:

osm2pgsql -l -c -d osm -U postgres -H localhost -S C:\osm2psql\default.style C:\temp\refine.osm

At this point I get the error

  Projection code failed to initialise

Since I encountered that error, I figured I had to install the proj v.4 lib and I added that to the environmental path, but I still get the issue.

How do you resolve this issue?

I need my data in 4326, not in 900913.  I think this is the reason why my data will not show in ArcGIS.

If I load the data without the -l command, which means lat/long, I get an invalid coordinate ID error from ArcGIS Desktop.


texts: Rotating ortho-satellite imagery so that square polygon perpendicular to the viewer

texts:
I'm trying to figure out how to rotate this satellite geotiff so that the rooftop polygon would appear as square as possible to the viewer rather than a diamond shape. See example below.

How to calculate the angle that the polygon and the image should be rotated to appear squarish? And how to rotate the geotiff by that amount (using gdal say)?

I've tried measuring the angles on the polygon with ST_Dump and ST_Angle but that returns nulls. I tried ST_Azimuth which returns values, but I don't understand how to use this information to solve my problem.

Example:

POLYGON((-104.947919 39.530436,-104.948051 39.530497,-104.948127 39.530398,-104.947995 39.530337,-104.947919 39.530436))


texts: Linear regresion at each raster pixel's value to predict future value (in R language)

texts:
I want to predict future value with existing time series raster. For simplicity, I want use linear regression at each raster pixel's value to predict future value
I have successfull run this code. I have read it from:
https://stackoverflow.com/questions/47435206/cant-calculate-pixel-wise-regression-in-r-on-raster-stack-with-fun?rq=1
library(raster)

# Example data
r <- raster(nrow=15, ncol=10)
set.seed(0)

# Now I make 6 raster (1 raster/months), then assign each pixel's value randomly
s <- stack(lapply(1:6, function(i) setValues(r, rnorm(ncell(r), i, 3))))
names(s) <- paste0('Month', c(1,2,3,4,5,6))

# Extract each pixel values
x <- values(s)

# Model with linreg
m <- lm(Month6 ~ ., data=data.frame(x))

# Prediction raster
p <- predict(s, m)

If you run that code, p will be a raster. But, I still confused. How to make raster in 'Month8' based on 6 previous raster?
What I mean is, each pixels has different linreg equations (where X=Month1, ..., Months6). If I input X=Month8, I will have 150 cells of Y for 8th Month that represent in each pixel of raster.
What I have done
# Lets try make a data frame for clear insight in my data
x <- values(s)
DF <- data.frame(x)

# Make X as month, and y is target.
library(data.table)
DF_T <- transpose(DF)
Month <- seq(1,nrow(DF_T))
DF_T <- cbind(Month, DF_T)

# Make prediction for first pixel
V1_lr <- lm(V1 ~ Month, data=DF_T)

# prediction for 8th Months in a pixel
V1_p <- predict(V1_lr, data.frame(Month=8))
V1_p

This is just one pixel. I want the entire raster


texts: How can I calculate the "average" or "median" result between a number of MultiLine Strings?

texts: So I have a number (let's say 12) of multiline strings representing something that changes over time, like forest borders or shorelines. I want to create a single line that is the result of the average or median calculation of those 12 strings. So if I have 1 ML String of a feature per month, I can calculate the average of those MLs for that year.


texts: Calculate distance between centroid and edge of polygon on given azimuth

texts:
I need to find the distance between the centroid and the edge of polygon on given azimuth.
One solution, that comes to my mind is to create a line starting on centroid and make it long enough to cross the edge (so for example make it a fixed 100km or come up with something more elaborate) and then split it by polygon.
That way, first segment gives me desired length.
However, this solution seems to be not very reliable and unnecessarily complex. Is there a simpler, better way to compute such distance?


texts:
I've been trying to add waypoints to Arcmap for some while now. I have a Mac OS X 10.6.8 and have finally found GPSBabel that is compatible with mac. My current problem is to convert the waypoints from the gps to a format which allows me to add the waypoints to Arcmap ( I don't know which output format to choose and convert it to a shapefile or x;y data to add to Arcmap)

Note: DNRgps is NOT compatible with Mac


texts: Adding Properties to a feature in Google Earth Engine

texts: Is there a way to add a property to a feature in the GEE. I know that you can edit existing properties with set(), but haven't found a way to add a property.


texts:
I'm tring to read a raster map from my grassData.

I'm using this code

import os
import grass.script.setup as gsetup
from grass.pygrass import raster

gisbase = os.environ['GISBASE']
gisdbase = os.path.abspath("/home/myuser/grassData")
location = 'mylocation'
mapset = 'mymapset'
gsetup.init(gisbase, gisdbase, location, mapset)

elev=raster.RasterRow('elevation')
elev.open()
for row in elev:
    print row

But I get all values with nan, where is the error?

If I do elev.info() I get this result

elevation@
rows: 4502
cols: 5554
north: 45.0743663458 south: 45.0500262069 nsres:5.40651685671e-06
east:  7.67457811722 west: 7.6445503225 ewres:5.40651687473e-06
range: 233.050003052, 310.049987793
proj: 3

Tnx

Lorenzo


texts: Automatically opening PDF output report from ArcGIS web app?

texts:
We have a custom tool that we built using the R-ArcGIS bridge. It queries a bunch of tables stored in a SQLite database and uses LaTeX to create a PDF report as the output.

Within ArcMap or ArcGIS Pro, the tool runs great and is programmed within the script to automatically open the saved PDF report using this code and it automatically opens fine:

system(paste0('open "', pdf.path, '"'))

We have created a geoprocessing tool from our server so that we can run this as a web app. It works great and saves the PDF file where we want. However, the PDF report that is generated is not being automatically opened as a popup when running it from the web application.

Any ideas for how to open a PDF saved on the server from the custom tool in the web app?


texts: Using gdalwarp with a crs which doesn't have an EPSG number

texts:
I need to project a raster from EPSG:32642 to this crs.

It is a relatively new addition to the EPSG registry and the EPSG number is not yet in gdal (at least in my version: 1.11.3):

steph@steph-W840AU:~$ gdalsrsinfo "EPSG:7694"
ERROR 1: ERROR - failed to load SRS definition from EPSG:7694

I thus wanted to use gdalwarp with the WKT definition:

gdalwarp -overwrite
-s_srs 'PROJCS["WGS 84 / UTM zone 42N",
GEOGCS["WGS 84",DATUM["WGS_1984",
SPHEROID["WGS84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",69],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1,AUTHORITY["EPSG","9001"]],
AUTHORITY["EPSG","32642"]]'
-t_srs 'PROJCRS["Kyrg-06 / zone 3",
BASEGEODCRS["Kyrg-06",    DATUM["Kyrgyzstan Geodetic Datum 2006",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1.0]]]],
CONVERSION["Kyrgyzstan zone 3",
METHOD["Transverse Mercator",ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.01745329252]],
PARAMETER["Longitude of natural origin",74.516666666667,ANGLEUNIT["degree",0.01745329252]],
PARAMETER["Scale factor at natural origin",1,SCALEUNIT["unity",1.0]],
PARAMETER["False easting",3300000,LENGTHUNIT["metre",1.0]],
PARAMETER["False northing",14743.5,LENGTHUNIT["metre",1.0]]],
CS[cartesian,2],
AXIS["easting (E)",east,ORDER[1]],
AXIS["northing (N)",north,ORDER[2]],
LENGTHUNIT["metre",1.0],  ID["EPSG",7694]]'
/home/steph/Temp/python_tmp/input/rgbir.tif
/home/steph/Temp/python_tmp/output/20160816-1117-UAVanalyse/rgbir_proj.tif

Which returns:

ERROR 1: Translating source or target SRS failed:

The WKT definition is copy-paste from here, so I assumed it should be working?

By the way, the definition of 32642 is not the problem: if I try to project to a well-known crs, then all works smoothly:

gdalwarp -overwrite
-s_srs 'PROJCS["WGS 84 / UTM zone 42N",
GEOGCS["WGS 84",DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",69],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1,AUTHORITY["EPSG","9001"]],
AUTHORITY["EPSG","32642"]]'
-t_srs 'EPSG:4326'
/home/steph/Temp/python_tmp/input/dsm.tif
/home/steph/Temp/python_tmp/output/20160816-1127-UAVanalyse/dsm_proj.tif


texts: combine two FOR LOOPS (MakeFeatureLayer_management)

texts:
I want to combine two FOR LOOPS but I have no idea how.
It should give me in each  TYPE Layer "under" Layers, where are new Layers for each MONTH (see picture)

The two loops in itself works.

# MakeFeatureLayer_management for each whale/dolphin type
whale_names = [
"Atlantic Spotted Dolphin",
"Pilot Whale",
"Bottlenose Dolphin",
"Rough-Toothed Dolphin",
"Bryde´s Whale",
"Common Dolphin"]

 for name in whale_names:
     arcpy.MakeFeatureLayer_management("Sichtungen95-14", "test-" + name, "SPECIES = '%s'" % (name))
     # MakeFeatureLayer_management for each month
 for month in range (01,12):
     arcpy.MakeFeatureLayer_management ("Sichtungen95-14", "test%s" % month, "DATE >= date'01.01.1995' AND DATE < date'01.01.2015' AND EXTRACT(MONTH FROM DATE) = %s" % (month))


texts: Modify the same feature multiple times ArcGIS Pro SDK EditOperation

texts:
I am atttempting to modify a feature twice using EditOperation in a loop. I am redacting certain subcurves from a given polyline at different sections of the line. The output polyline should be the original line without the two subcurves redacted by the for loop. But as you can see from the image below, only the redaction from the second loop is performed. Using the undo stack, I was able to see that after the first loop, the polyline is indeed modified correctly but on the second modification,(editop.Modify(poly, oid, redactedPoly);) instead of using the modified feature in the first loop as the source poly, it seems like it's using the original unmodified polyline to perform the modification in loop two. How can I get it to modify the polyline correctly. I also attempted to use chained edit operation but the results were the same.

var poly = MapView.Active.Map.GetLayersAsFlattenedList().First() as FeatureLayer;
using (var TargetCursor = poly.Search(queryFTarget))
                                {
                                    while (TargetCursor.MoveNext())
                                    {
                                        targetFeat = TargetCursor.Current as Feature;
                                        var oid = targetFeat.GetObjectID();
                                        var targetPoly = targetFeat.GetShape() as Polyline;
                                        if (targetPoly == null || targetPoly.IsEmpty) //if the target feature is null or empty
                                        {
                                            file.WriteLine($"Unable to redact '{sourceID}' and '{targetID}'.");
                                            file.WriteLine($"Feature '{targetID}' has no shape!");
                                        }
                                        else
                                        {
                                            var distancesAlongCurve = GetDistancesAlongCurve(targetPoly, ptCollection);//Get the distances along the source line for each demarkation point as ratios
                                            distancesAlongCurve.Sort();

                                            for (var i = 0; i < distancesAlongCurve.Count; i += 2)
                                            {
                                                var GetsubCurve = GeometryEngine.Instance.GetSubCurve(targetPoly, distancesAlongCurve[i], distancesAlongCurve[i + 1], AsRatioOrLength.AsRatio);
                                                var subCurve = GeometryEngine.Instance.SimplifyAsFeature(GetsubCurve) as Polyline;

                                                if (subCurve != null && !subCurve.IsEmpty)
                                                {
                                                    var difference = GeometryEngine.Instance.Difference(targetPoly, subCurve);
                                                    var redactedPoly = GeometryEngine.Instance.SimplifyAsFeature(difference) as Polyline;

                                                    var editop = new EditOperation();
                                                    editop.Name = "Redact Polyline";
                                                    editop.Modify(poly, oid, redactedPoly);
                                                    await editop.ExecuteAsync();
                                                }
                                            }
                                        }
                                    }
                                }


texts: Add and Remove layers in openlayers on-the-fly

texts:
I am using OpenLayers 2.13.1 in my application and I have a couple of layers which need to be added/removed to the map on-the-fly.

Here are the core codes:

var layers=[layer1,layer2,layer3];  // both of them are instances of ol.Layer.WMTS
var map=.....;
removeLayer();
map.events.register("moveend",map,resetLayer);  // when the center of the map changed, re-decide which layer should be visible
.....

function resetLayer(){
    var activated=[];

    var bounds=map.getExtent(),zoom=map.getZoom()
    for(var i=0,len=layers.length;i<len;i++){
        var lay=layers[i];
        if(checkVisible(bounds,zoom,layer1)){
            if(!lay.activated){
                map.addLayer(lay);   // #1
                lay.activated=true;
            }
            activated.push(lay);
        }else{
            map.removeLayer(lay);
            lay.activated=false;
        }
    }

    if(activated.length==0){
        //no layers avaiable,add the first layer
        var lay=layers[0];
            if(!lay.activated){
                map.addLayer(lay);  //#2
                lay.activated=true;
            }
    }
}

Now, I came to the too much recursion problem.

It seems that when the moveend event fired, the resetLayer will be called (expected). Once I re-add/remove the layer to the map (see comments of #1 and #2), the map object will reset the baselayer which will cause the setCenter to be called for the map, and the moveend event will be fired once again.

Any ideas on how to fix this?


texts: Range of easting value for each MGRS zone(UTM zone)

texts:
From MGRS grid reference which also can be found from https://www.maptools.com/tutorials/100km_square_id , the minimum column id for the first zone is 'A' and it's given when easting is smaller than 200km.

Is this mean that it is guaranteed that UTM easting value in a zone is always larger than 100km? Or, should I consider the width of first column is 200km?

Similarly, for the last column 'H' in the first zone, is it guaranteed that any UTM easting value in a zone is smaller than 900km? Or should I just handle any value larger than 800km as 'H' column?


texts: Plotting Single Band & Multi-Band Stacked TIFF File using Rasterio in Python

texts:
I have 2 raster files (one single band and other multi-band stacked (465 bands) tiff file) which I am trying to plot using the following code:

import rasterio as rio
from rasterio.plot import show

## Plotting Single-Band Image
with rio.open(r'C:\Users\user\Desktop\et.tif') as src_plot:
    show(src_plot)

## Plotting Multi-Band Stacked Image
with rio.open(r'D:\Weather data\et_stack.tif') as src_plot:
    show(src_plot.read([1,2,3]))

Using the above code both the images are getting displayed correctly. But the issue which I am facing is that when I am plotting the stacked image using bands 1, 2 & 3 (or any other band combinations), the image is getting displayed but with black background. This is not the case when I am plotting single band image i.e. single band image gets displayed with white background using the above code.

Can someone please help me out as to how I can display the image (stacked image) without black background.


texts: set attributions to WMTS Layer from capabilities in Openlayers

texts:
I'm loading different layers from one WMTS source and want to display the correct attributions for the source of each layer.

In my code I define the attribution at the source-level for a ol.layer.Tile. The layer is later added to the map and display correct, but when I click the i-button on the map it only displays 'undefined'. In the console I don't get any errors.

var parser = new ol.format.WMTSCapabilities();

  fetch('http://wmts1.geoportail.lu/opendata/wmts/1.0.0/WMTSCapabilities.xml').then(function(response) {
    return response.text();
  }).then(function(text) {
    var result = parser.read(text);

    var hybrid_source = ol.source.WMTS.optionsFromCapabilities(result, {
          layer: 'hybrid',
          matrixSet: 'EPSG:3857',
          attributions: [new ol.Attribution({
                    html: 'data.public.lu Hybrid @ ACT Luxembourg'
                    })]
    });

    var hybrid_layer = new ol.layer.Tile({
          opacity: 1.0,
          title: 'Hybrid',
          baseLayer: true,
          visible: true,
          source: new ol.source.WMTS(/** @type {!olx.source.WMTSOptions} */ (hybrid_source))
        });
  });


texts:
I am using cartodb.js to import one of our maps via cartodb.createVis. I am working on adding custom functionality outside the scope of the web console's customizations (for example, custom search fields).

I want to write end-to end tests for cartodb's functionality. I want to select specific markers, trigger click events on them, and then expect DOM changes and map changes. I also want to be able to test for "markers" or map elements to appear or disappear

The problem is that I cannot use selectors to access individual markers or any of the functionality provided by the "visualization" outside of the search bar and layer selector. I have not been able to construct the appropriate CSS selectors, and selenium will not record click actions on map elements (except for double-clicks on map tiles, which cause a zoom, and aren't useful for me).

I can test for the existence of a tile by using a selector for its image. This doesn't work because:

a) any slight change to the visualization changes this link and breaks the test (I'm not sure the link is reliable even if we don't change the visualization)

b) I can't test for the existence of specific map elements (e.g., markers).

How can I write tests for our visualization?

Update: I tried using layer.trigger('featureClick' . . . as described here and here but that is failing silently. I'm not sure if it matters, but unlike that Gist I have multiple layers and my coordinates are not quite as precise.

I am also trying to trigger marker clicks using Selenium's clickAt method. For this I am using the example here I've tried:

td>clickAt
//div[@id='map']/div
601, 230

This fails silently.

clickAt
//img[contains(@src,'http://3.ashbu.cartocdn.com/documentation/api/v1/map/c60783e45b3de44d88d22b74986b35f6:0/6/31/24.png')]
601, 229

Gives:
TypeError: _d is undefined cartodb.uncompressed.js:14231


texts: C# - convert lat long coordinates to WKT format

texts:
How can I convert coordinates to WKT format like 39.8251136,-8.3824103 with 900913 system (spherical Mercator)?

My client is askying me something like this: SRID=900913;POINT(-939520.0088637
4744291.7338414) but I already search about this and I don't find anything.

I found this link: https://stackoverflow.com/questions/3561186/converting-lat-long-to-pointclass but I don't know where I can find the DLL with the SpatialReferenceEnvironmentClass class for I can test that code...


texts: Console application in .NET incorporating QGIS tools

texts: I am working on a project where I want to build a standalone application in .NET using QGIS tools. Is there any option like the ArcObjects SDK for .NET? Please help me on this.


texts: Grid size of a Digital Elevation Model developed from topographic surveyed elevation contour line?

texts:
I have generated a Digital Elevation Model from elevation contour line. I got these elevation contour line at 20 cm interval from DGPS survey.

What would be the grid size of generated DEM in "m".
For example in SRTM the spatial resolution or grid size is 90 m that means for 90m X 90m we have one elevation value.


texts: Draw polygon on load not working - OpenLayers

texts:
The following code should have 3 features:

Draw a polygon on load
On draw end send coordinates
On modify end send coordinates

2&3 work as expected. I added the drawn polygon last, there is no error on the console, but no polygon is drawn.

Also, it has no errors building in npm.

import Map from 'ol/Map.js';
      import View from 'ol/View.js';
      import {Draw, Modify, Snap} from 'ol/interaction.js';
      import {Tile as TileLayer, Vector as VectorLayer} from 'ol/layer.js';
      import {OSM,TileJSON, Vector as VectorSource} from 'ol/source.js';
      import {Circle as CircleStyle, Fill, Stroke, Style} from 'ol/style.js';

      import Feature from 'ol/Feature';
      import Polygon from 'ol/geom/Polygon';
    import Point from 'ol/geom/Point';
    import Transform from 'ol/Proj';

      var raster = new TileLayer({
        source: new OSM()
      });

      var source = new VectorSource();
      var vector = new VectorLayer({
        source: source,
        style: new Style({
          fill: new Fill({
            color: 'rgba(255, 255, 255, 0.2)'
          }),
          stroke: new Stroke({
            color: '#ffcc33',
            width: 2
          }),
          image: new CircleStyle({
            radius: 7,
            fill: new Fill({
              color: '#ffcc33'
            })
          })
        })
      });
      var feature = new Feature({
        geometry: new Polygon([
            [-13041991.514129914,8126540.3085058555],
            [-13766003.046047103,7167714.225696605],
            [-12562578.472725287,6668733.305050973],
            [-11486345.114470007,7099226.648353085],
            [-11730943.60498257,8439626.376361936],
            [-13041991.514129914,8126540.3085058555]
          ])
    });
    var vectorSource= new VectorSource({
        features: [feature ]
    });

    var vectorLayer = new VectorLayer({
        source: vectorSource
    });

      var map = new Map({
        layers: [raster, vector],
        target: 'map',
        view: new View({
          center: [-11000000, 4600000],
          zoom: 4
        })
      });

      var modify = new Modify({source: source});
      modify.on('modifyend',function(e){
        console.log("feature id is",e.features.getArray()[0].getGeometry().getCoordinates()[0]);
        });
      map.addInteraction(modify);
      console.log(map);
      var draw, snap; // global so we can remove them later
      var typeSelect = document.getElementById('type');

      function addInteractions() {
        draw = new Draw({
          source: source,
          type: typeSelect.value
        });
        draw.on('drawend',function(e){
            console.log(e.feature.getGeometry().getCoordinates()[0]);
        })
        map.addInteraction(draw);
        snap = new Snap({source: source});
        map.addInteraction(snap);
      }

      /**
       * Handle change event.
       */
      typeSelect.onchange = function() {
        map.removeInteraction(draw);
        map.removeInteraction(snap);
        addInteractions();
      };

      addInteractions();


texts: Object given is not a Python function - gdal.RasterizeLayer

texts:
I'm running into this error and I'm not quite sure what the error means. Novice here, operating out of jupyter environment, and trying to burn the vector layer of roads (road_lyr) into the raster road_ds.

gdal.RasterizeLayer( road_ds, [1], road_lyr, burn_values=[1], callback=gdal.TermProgress)

where my road_ds is

road_ds = gtiff_driver.Create('roads_prox.tif', in_band.XSize, in_band.YSize)
road_ds.SetProjection(in_ds.GetProjection())
road_ds.SetGeoTransform(in_ds.GetGeoTransform())

and road_lyr is

shp_ds = ogr.Open(folder)
road_lyr = shp_ds.GetLayerByName(roads_ln)
road_lyr.GetExtent()

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-16-b8a413aa8f73> in <module>()
----> 1 gdal.RasterizeLayer( road_ds, [1], road_lyr, burn_values=[1], callback=gdal.TermProgress)
C:\Anaconda3\lib\site-packages\osgeo\gdal.py in RasterizeLayer(*args, **kwargs)
   1622 def RasterizeLayer(*args, **kwargs):
   1623     """RasterizeLayer(Dataset dataset, int bands, Layer layer, void * pfnTransformer=None, void * pTransformArg=None, int burn_values=0, char ** options=None, GDALProgressFunc callback=0, void * callback_data=None) -> int"""
-> 1624     return _gdal.RasterizeLayer(*args, **kwargs)
   1625
   1626 def Polygonize(*args, **kwargs):
RuntimeError: Object given is not a Python function

How do I fix this?


texts: Generate GeoServer WMS layer from multiple shapefiles

texts:
I have multiple shapefile layers that I intend to use to create a WMS basemap with GeoServer.

What is the best approach for this?

I would first generate raster image from my shapefile layers, then I would generate raster basemap for multiple zoom levels: https://docs.qgis.org/2.8/en/docs/user_manual/processing_algs/gdalogr/gdal_miscellaneous/overviews.html

After that, I would load basemap to postgis database and serve it with geoserver.

Is this a good approach?


texts: Calculate the sum of a field from shapefiles and copy these values to a new table

texts:
I have a number of shapefiles, the attribute table of each  one has a FLOAT field called "Relative_W". I need to create a table that displays the summation of "Relative_W" for all the shapefiles. ANY HELP?

the shape files attribute as this

Unfortunately, the new table is always empty

here is my code.

# Create the Table to be filled with the summations
arcpy.env.workspace=r'C:\ModelBuilder\AAA\New Folder (17)\Flash Flood Risk Assessment Projcet\FFRA_WRRI.gdb'
fctble=arcpy.CreateTable_management(arcpy.env.workspace,"Scenarious_Comparison")
arcpy.AddField_management(fctble, "Scenario", "TEXT", field_length = 8)
arcpy.AddField_management(fctble, "Risk_Indicator", "FLOAT")

for i in range (0,2):
    sum = 0
    fc ='Risk_Scenario_%d' %i # the shape files that I need to read and sum its value
    with arcpy.da.SearchCursor(fc, "Relative_W") as cursor: # Relative_W the field of needed data
        for row in cursor:
            if row[0] is not None:
                sum = sum + row[0] # This is the summation value.
                                   # NOW, I need to add this value to the generated empty table "fctble"
    with arcpy.da.UpdateCursor(fctble,["Scenario","Risk_Indicator"]) as cursor:
        for row in cursor:
                cursor.updateRow([i,summed_total])


texts: Leaflet.MultiTileLayer : problem with my TMS : Y coordinates are bad

texts:
I have installed leaflet to create a web site for consulting data created by
firefighters in France.

I want to test the plugin leaflet.multitilelayer  for defining a TileLayer in terms of several TileServer urls. I have tested this plugin with data of others: it works (the example is here : http://makina-corpus.com/blog/metier/2013/infrastructure-webmapping-minimaliste-partie-3-leaflet). But when I have tested with our data (with server of spatial data GeoServer and MapServer to create TMS), it doesn't work. It seems that the coordinates are note the good, because I have compared with the simple L.tilayer.

Here is our code :

 var PLAN  =
  L.TileLayer.multi({
    17: { // this layer is shown at zoomlevels between maxZoom of
          // previous layer and 17 (ie 17)
        url: 'http://xxx.XXX.XXX.XXX/geowebcache/service/tms/1.0.0/[our_data]/{z}/{x}/{y}.png',
    },
     19: { // this layer is shown at zoomlevels between minZoom and 19 (ie 18,19)
        url: 'httpxxx.XXX.XXX.XXX/geowebcache/service/tms/1.0.0/[our_data]/{z}/{x}/{y}.png',
    }
}, {
minZoom: 17,
maxZoom: 19,
});

When I use the command L.tileLayer for these TMS, it works well.

Would anyone ever had this problem?


texts: How do you determine a shapefile's charset?

texts: Is there a way to determine the charset used for a given shapefile?


texts: Custom styling for cluster images in OpenLayers 3?

texts:
new with OpenLayers 3.

While it is possible to generate simple polygons and circles (like in the cluster example) I wish to generate a more complex Cluster style that varies depending on the features contained in the cluster.

For instance, I would like to generate a dynamic pie chart based on values. It is possible to render this client-side using HTML5 canvas function (arc()) and friends, but I cannot find a way to generate such a drawing using generic canvas primitives within OL3's styles.

Can anybody point me in the right direction? If impossible, I'll "enhance" OL3 with the missing functionality...


texts: Is it possible to add a weight factor to an OD-Matrix?

texts: Is it possible to add a weight factor to OD-Matrices? I am looking to compute a matrix which shows origin points going to the nearest destination, based on this destination's size (with preference given to bigger destinations), within a given timeframe. Thanks in advance.


texts: Is Tomcat or Jetty a better container for Geoserver on Windows

texts:
My Understanding is that Geoserver can be run in the default Jetty Environment or in Tomcat on Windows Server 2008 R2
GeoWebCache uses Tomcat as its container in its installation instructions.
Which is better for Geoserver, Tomcat or Jetty and why?
UPDATE (Thanks @BradHards)

Application is stand alone Geoserver talking via Private LAN to a Standalone Postgres/PostGIS Server delivering up to 4 WMS tiled layers to the internet
Up-time needs to be 99%
Load is undetermined at this time. Deployment of additional servers is expected when load reaches a predetermined limit, which we are trying to assess.
Environment: Windows Server 2008 (32 Bit), Public Facing. Currently running Geoserver on Jetty
Updates: only when they will provide performance benefits or additional features that will improve the WMS Layers delivered.


texts: Extract raster by mask output does not completely overlap mask?

texts:
I have a flow accumulation grid and a shapefile containing catchment areas. Now I want to do some calculations on them so I need to extract data from cells in the raster adjacent to certain catchment areas.

In that process I try to use Extract by Mask on my raster, with one of the catchment polygons as mask (in a loop later on, but now I just need it to work once).

The problem is that the extracted raster is shifted half a cell to the left, compared to both the original raster and the masking polygon. This is messing up the rest of my process, since I try to convert the raster to points and intersect with other polygons later on. And unless the extracted raster is overlapping the original one, I wont get the right polygon.

The original raster and the polygons are aligning perfectly, since the polygons are created from the raster. I have tried to check if there somehow have been some kind of change in the projections, but can't find any. I have also tried in new ArcMap documents, restarting everything several times (this solves the problem surprisingly often) and executing the tools from the python prompt instead of calling my script. I have even tried to execute the Extract by Mask tool from the Toolbox, with same result.

This is the (very simple) code that I use in the python prompt in ArcMap:

inpgs = r'C:\GIS\Catchm_ex.shp'
flowacc = r'C:\GIS\flowacc'
arcpy.MakeFeatureLayer_management(inpgs, 'currentMask', '"HydroID" = 269634')
outRaster = arcpy.sa.ExtractByMask(flowacc, 'currentMask')

I have some sample files here

I use Windows 7, ArcInfo 10.0 and Python 2.6.

I have tried Extract by Polygon now, and that seems to work fine. Is it just my machine malfunctioning on ExtractByMask?
Only drawback is that I'm having problems using that function in my script, since I can't work out how to get the coordinates of my masking polygon (only covering 9 cells, with a point in center if that helps).

I have now got it working with Extract by polygon by using the code below (modified from Esri's support pages). Unfortunately, this does not solve my problem. It appears as if the output (from this tool as well as Extract by mask) is sometimes shifted half a cell to the left, and sometimes not. Mostly it is, but occasionally it actually is located where it should be.

So now I'm really confused. Is this just happening to me, on my computer? If anyone else have experienced the same, please write a comment so I know that I'm not all alone :) To me this problem seems to be a too big of an issue to be a bug, which leaves me with this being a problem that I have created myself along the way.

Any ideas on how to avoid it?

def ReadPoints(infc):
    try:
        import arcpy, numpy
    except ImportError:
        sys.exit('Unable to import arcpy')

    # Identify the geometry field
    desc = arcpy.Describe(infc)
    shapefieldname = desc.ShapeFieldName
    points = []

    # Create search cursor
    rows = arcpy.SearchCursor(infc)

    # Enter for loop for each feature/row
    for row in rows:
        # Create the geometry object
        feat = row.getValue(shapefieldname)
        partnum = 0

        # Step through each part of the feature
        for part in feat:
            # Step through each vertex in the feature
            for pnt in feat.getPart(partnum):
                if pnt:
                    points.append([int(pnt.X), int(pnt.Y)])

    pt = numpy.array(points)
    pts = []
    for i in range(len(pt)):
        t = pt[i]
        pts.append(arcpy.Point(t[0],t[1]))
    return pts


texts: Is there a QGIS equivalent to ArcGIS's Make OD Cost Matrix Layer (Network Analyst)?

texts:
I'm looking to create a OD cost distance matrix for a network data set I have created. Is there a way of doing this in QGIS which is akin to Make OD Cost Matrix Layer within ArcGIS's Network Analyst.

My question is similar to this but there was no satisfactory answer given then. That question was asked nearly months ago and am wondering if there is a potential solution.


texts:
The one-to-many relations in QGIS is OK with forms, works great.
But how to labels or chart data with this one-to-many relations?


texts: Remove last 5 characters with Field Calculator in ArcGIS Desktop?

texts:
I wonder how to delete a word (or the last 5 characters) using the field calculator in ArcGIS 10.2.
I have a column such as:

Paris City;
London City;
Berlin City

I need to remove "City", to be left with Paris, London, Berlin.

Any tips on doing this automatically?