texts: How to solve UnicodeDecodeError in Python when executing Make Valid LWGEOM algorithm in QGIS?

texts:
I'm trying to fix some invalid geometries using the Make Valid LWGEOM algorithm in QGIS 2.8.3 processing toolbox but a Python error occurs:

Errore durante l'esecuzione di codice Python:

Traceback (most recent call last):
  File "C:/Users/Angelo/.qgis2/python/plugins\processing\gui\AlgorithmDialog.py", line 200, in accept
    if runalg(self.alg, self):
  File "C:/Users/Angelo/.qgis2/python/plugins\processing\gui\AlgorithmExecutor.py", line 51, in runalg
    alg.execute(progress)
  File "C:/Users/Angelo/.qgis2/python/plugins\processing\core\GeoAlgorithm.py", line 248, in execute
    str(e) + self.tr('\nSee log for more details'))
UnicodeDecodeError: 'ascii' codec can't decode byte 0xe8 in position 19: ordinal not in range(128)

Versione Python:
2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]

Versione di QGIS:
2.8.3-Wien Wien, 313a6bc

Percorso Python: ['C:/Users/Angelo/.qgis2/python/plugins\processinglwgeomprovider', 'C:/Users/Angelo/.qgis2/python/plugins\processing', 'C:/PROGRA~1/QGISWI~1/apps/qgis-ltr/./python', u'C:/Users/Angelo/.qgis2/python', u'C:/Users/Angelo/.qgis2/python/plugins', 'C:/PROGRA~1/QGISWI~1/apps/qgis-ltr/./python/plugins', 'C:\PROGRA~1\QGISWI~1\apps\Python27\lib\site-packages\matplotlib-1.3.1-py2.7-win-amd64.egg', 'C:\PROGRA~1\QGISWI~1\apps\Python27\lib\site-packages\nose-1.3.3-py2.7.egg', 'C:\PROGRA~1\QGISWI~1\apps\Python27\lib\site-packages\tornado-4.0.1-py2.7-win-amd64.egg', 'C:\PROGRA~1\QGISWI~1\apps\Python27\lib\site-packages\backports.ssl_match_hostname-3.4.0.2-py2.7.egg', 'C:\PROGRA~1\QGISWI~1\apps\Python27\lib\site-packages\certifi-14.05.14-py2.7.egg', 'C:\PROGRA~1\QGISWI~1\bin\python27.zip', 'C:\PROGRA~1\QGISWI~1\apps\Python27\DLLs', 'C:\PROGRA~1\QGISWI~1\apps\Python27\lib', 'C:\PROGRA~1\QGISWI~1\apps\Python27\lib\plat-win', 'C:\PROGRA~1\QGISWI~1\apps\Python27\lib\lib-tk', 'C:\PROGRA~1\QGISWI~1\bin', 'C:\PROGRA~1\QGISWI~1\apps\Python27', 'C:\PROGRA~1\QGISWI~1\apps\Python27\lib\site-packages', 'C:\PROGRA~1\QGISWI~1\apps\Python27\lib\site-packages\PIL', 'C:\PROGRA~1\QGISWI~1\apps\Python27\lib\site-packages\jinja2-2.7.2-py2.7.egg', 'C:\PROGRA~1\QGISWI~1\apps\Python27\lib\site-packages\markupsafe-0.23-py2.7-win-amd64.egg', 'C:\PROGRA~1\QGISWI~1\apps\Python27\lib\site-packages\pytz-2012j-py2.7.egg', 'C:\PROGRA~1\QGISWI~1\apps\Python27\lib\site-packages\win32', 'C:\PROGRA~1\QGISWI~1\apps\Python27\lib\site-packages\win32\lib', 'C:\PROGRA~1\QGISWI~1\apps\Python27\lib\site-packages\Pythonwin', 'C:\PROGRA~1\QGISWI~1\apps\Python27\lib\site-packages\Shapely-1.2.18-py2.7-win-amd64.egg', 'C:\PROGRA~1\QGISWI~1\apps\Python27\lib\site-packages\six-1.3.0-py2.7.egg', 'C:\PROGRA~1\QGISWI~1\apps\Python27\lib\site-packages\wx-2.8-msw-unicode', 'C:\PROGRA~1\QGISWI~1\apps\Python27\lib\site-packages\xlrd-0.9.2-py2.7.egg', 'C:\PROGRA~1\QGISWI~1\apps\Python27\lib\site-packages\xlwt-0.7.5-py2.7.egg', u'C:/Users/Angelo/.qgis2//python', 'C:\Users\Angelo\.qgis2\python\plugins\mmqgis/forms', 'C:\PROGRA~1\QGISWI~1\apps\qgis-ltr\python\plugins\fTools\tools']

How can I fix it? Or which other method could I use to fix the invalid geometries? I have to fix some shapefiles which may contain hundreds of invalid ones so I can't do it manually... I also tried the same procedure in QGIS 2.10.1 but the same error occurs...


texts:
I use QFIELD on Android 5.1.1 (Samsung Galaxy), I tried several versions of QFIELD without solving my problem:
when I click on a vector I just created, the attribute window opens on the right side of the screen and never slides back ... and hides the pencil button which allows to digitize!
Has someone had this problem before ?
I created my project on QGIS 2.12.3


texts: R: Return the UTM zone that a WGS84 point belongs to

texts:
I'm looking for an R function that might work similarly to this website that will tell you what UTM zone your WGS84 point belongs to.

http://www.rcn.montana.edu/resources/converter.aspx


texts: Is there a way in QGIS to dissolve a vector by two or more fields?


texts: CRS and QGIS georeferencing for Cassini projection

texts:
I have scans of OS 6" series maps which I understand are Cassini projection.  According to an article my Mugnier, Photographic & remote sensing, Oct 2003 they had  the typical scale factor equal to unity, no false origin, and a single triangulation station as the projection origin.  From this I take it that the proj4 string should be something along the lines of:

 +proj=cass +lat_0=53.9619966389 +lon_0=-1.0804806111 +x_0=0 +y_0=0 +units=m +no_defs

where lat_0 and lon_0 values are the coordinates of the origin, in this case York Minster.
If I then use the QGIS georeferencer with a linear transform with control points taken from the printed coordinates of the corners the input and output set to this CRS I expect it to generate a world file for the scan with the coordinates are the top left of the sheet, x & y rotations of zero and x and y pixel values which are equal except that y is negative.  In fact the pixel values are quite different so that when QGIS opens the file the proportions are grossly distorted.
Clearly something in my understanding of this is wrong, but what?

These are points very near the corners of the printed map.  The actual scan file is somewhat larger as it includes the entire sheet.  Projections are not printed.  The 6" county series had different projections for each county or, I think, in some cases for small groups of counties.  For Yorkshire this was, as above, centred on York Minster.  This particular sheet is 272.

-1.72916666666666652,53.72888888888888914
-1.87361111111111112,53.52027777777777828
-1.72833333333333328,53.52111111111111086


texts: T-SQL Invalid Geometry Created by Intersection

texts:
I have two sets of geometries (A and B) in SQL Server 2016. I want to determine the overlap areas of those geometries from set A that intersect set B.

Most of the time this code works:

  SELECT a.geom.STIntersection(b.geom)
  FROM A a
      JOIN B b ON a.geom.STIntersects(b.geom) = 1

But in a few cases, I get the following error:

  Msg 6522, Level 16, State 1, Line 1
  A .NET Framework error occurred during execution of user-defined routine or aggregate "geometry":
  System.Runtime.InteropServices.COMException: Error HRESULT E_FAIL has been returned from a call to a COM component.
  System.Runtime.InteropServices.COMException:
     at System.Runtime.InteropServices.Marshal.ThrowExceptionForHRInternal(Int32 errorCode, IntPtr errorInfo)
     at Microsoft.SqlServer.Types.GLNativeMethods.ThrowExceptionForHr(GL_HResult errorCode)
     at Microsoft.SqlServer.Types.GLNativeMethods.Combine(CombineMode combineMode, GeoData g1, GeoData g2)
     at Microsoft.SqlServer.Types.SqlGeometry.STIntersection(SqlGeometry other)

I have isolated the offending geometries by iterating through one at a time. It looks like the problem arises when two geometries touch but do not overlap. However, SQL shows that STTouches = 0. I suspect what is happening is that some kind of defective geometry is being created, like maybe a geometry collection?

Both geometries are valid. Adding .MakeValid() does not affect the results. Both geometries are polygons or multipolygons.

My question is, how can I prevent my query from erroring out? I am happy to just skip these cases if I can identify them somehow. I'd like to avoid iterating through all of them one by one, as this is very inefficient.


texts: Do rasters have to be of same extent for the QGIS Raster Calculator?

texts: I am trying to use the Raster Calculator but, I get a raster with Nan values. My rasters are in the same projection system but, one of them, it is a classified landsat map which covers only some part of the other map that I am trying to use. Can I actually do that?  They should be of the same extent?


texts: How to work with GeoServer and OpenLayers?

texts:
I am new to GeoServer and am a bit lost at the moment.

I have added shapefiles to GeoServer and created a layer containing all the shapefiles that I would like to display on my externally hosted web site.  When I go to layer preview in GeoServer it allows you to display an OpenLayers preview in a pop up. I have copied the page source from the preview into the HTML of my hosted web site. This seems to work fine when I open my web page when GeoServer is running.  As soon as I turn off GeoServer I can no longer view these data on my web page.

Can anyone help me?


texts: Converting shp file to wkt polygon string using ArcPy?

texts:
I used the following code to convert .shp file to .wkt.

The code was successful but where can I find the .wkt file (the print row ouput)?

Is it possible to add a path for the wkt files to be stored in?

import arcpy
InFile = r"C:\TEST.shp"
for row in arcpy.da.SearchCursor(InFile,["SHAPE@WKT"]):
    print row[0]


texts: How I can use feature layer DefinitionExpression dynamically by using the javascript API 3.17?

texts:
My Flow is like below:

When area name are selected from the dropdown, the map layer is filtered by using DefinitionExpression of feature layer.

layer.setDefinitionExpression('area_name = xyz');
Till that this is working perfectly but there are no any layer refreshment by applying setDefinitionExpression second time. It mean that there are no any action happened after selecting area from dropdown second time.

Do you have any idea ?


texts: How to create dynamic point_coord variables for Neighborhood Selection?

texts:
I'd like to us the arcpy.GANeighborhoodSelection_ga in a loop where the point coords would come from a different source. What I tried already:
If I define the coords, it works:

pointcoord = "689383.6885 3973775.2178"
arcpy.GANeighborhoodSelection_ga(inFile2, "neighbor", pointcoord, "1", "1","20", "20", "0", "One sector")
Result 'neighbor'

When I try to use variable for the coords I got error:

pointcoord = x, y
print pointcoord
(689383.6885, 3973775.2178)
arcpy.GANeighborhoodSelection_ga(inFile2, "neighbor",  pointcoord , "1", "1","20", "20", "0", "One sector")
Runtime error
Traceback (most recent call last):
      File "<string>", line 1, in <module>
      File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\ga.py", line 1075, in GANeighborhoodSelection
        raise e
    RuntimeError: Object: Error in executing tool

pointcoord = "\"" + x + ' ' + y + "\""
Runtime error
Traceback (most recent call last):
  File "<string>", line 1, in <module>
TypeError: cannot concatenate 'str' and 'float' objects

pointcoord = "\"" + str(x) + ' ' + str(y) + "\""
print pointcoord
"689383.6885 3973775.2178"
arcpy.GANeighborhoodSelection_ga(inFile2, "neighbor", pointcoord, "1", "1","20", "20", "0", "One sector")
Runtime error  Traceback (most recent call last):   File "<string>", line 1, in <module>   File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\ga.py", line 1075, in GANeighborhoodSelection     raise e ExecuteError: ERROR 000622: Failed to execute (Neighborhood Selection). Parameters are not valid. ERROR 000628: Cannot set input into parameter point_coord.


texts: How to create a plot area polygon geometry into multiple plot area polygon geometry using OpenLayers

texts: I am using Mapnik, Mapproxy and OpenLayers 6 for the development of GIS based web-application. Now I need to divide the selected polygon area(land plot)into multiple as per user requirement using the GIS web application. The user want to divide the polygon plot area  into multiple along with textual information  by selecting in the GIS map How to develop the polygon splitting features as shaded in the attached picture?


texts: GeoServer OS MasterMap WMS topoarea error SLD

texts:
I'm getting an error when displaying Ordnance Survey MasterMap data through GeoServer as a WMS in QGIS.

The error is:

<?xml version="1.0" encoding="UTF-8"?>
<ServiceExceptionReport version="1.3.0" xmlns="http://www.opengis.net/ogc" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.opengis.net/ogc http://localhost:8080/geoserver/schemas/wms/1.3.0/exceptions_1_3_0.xsd">
    <ServiceException>
  java.lang.RuntimeException: org.xml.sax.SAXException: Entity resolution disallowed for SYSTEM org.xml.sax.SAXException: Entity resolution disallowed for SYSTEM Entity resolution disallowed for SYSTEM
    </ServiceException>
</ServiceExceptionReport>

Exception report:

<xsd:schema targetNamespace="http://www.opengis.net/ogc" elementFormDefault="qualified" version="1.3.0 2010-02-03">
<!--

      WMS is an OGC Standard.
      Copyright (c) 2004,2010 Open Geospatial Consortium, Inc. All Rights Reserved.
      To obtain additional rights of use, visit http://www.opengeospatial.org/legal/ -->

    <xsd:element name="ServiceExceptionReport">
        <xsd:complexType>
            <xsd:sequence>
                <xsd:element name="ServiceException" type="ogc:ServiceExceptionType" minOccurs="0" maxOccurs="unbounded"/>
            </xsd:sequence>
            <xsd:attribute name="version" type="xsd:string" fixed="1.3.0"/>
        </xsd:complexType>
    </xsd:element>
    <xsd:complexType name="ServiceExceptionType">
        <xsd:simpleContent>
            <xsd:extension base="xsd:string">
                <xsd:attribute name="code" type="xsd:string"/>
                <xsd:attribute name="locator" type="xsd:string"/>
            </xsd:extension>
        </xsd:simpleContent>
    </xsd:complexType>
</xsd:schema>

I'm very new at Geoserver and not sure about the error at all. When researching the error a lot of people get SAXParseException but can't find something similar to what I get.

I think it's something to do with the colour hex in the SLD file I'm using, sample below:

 <FeatureTypeStyle>
    <Rule>
      <Name>Unclassified - 1:0 to 1:4,000</Name>
      <ogc:Filter>
        <ogc:PropertyIsEqualTo>
          <ogc:PropertyName>style_code</ogc:PropertyName>
          <ogc:Literal>99</ogc:Literal>
        </ogc:PropertyIsEqualTo>
      </ogc:Filter>
      <MinScaleDenominator>0</MinScaleDenominator>
      <MaxScaleDenominator>4000</MaxScaleDenominator>
      <PolygonSymbolizer>
        <Fill>
          <CssParameter name="fill">#F7F5F0</CssParameter>
        </Fill>
      </PolygonSymbolizer>
    </Rule>
  </FeatureTypeStyle>


texts: Generating an irregular polygon buffer around a line in QGIS

texts:
I want to create a cartographic effect to represent a hiking trail. I thought that it may look cool to generate an irregular buffer around the line representing the trails so it looks more "natural".

I think it can be done with the Geometry Generator, as it can generate a simple buffer around the line, but I don't know if is possible to generate an irregular shape.

Does anyone have a clue? I am also open to other ideas to map hiking trails.

EDIT: Here is a (terrible) sketch of what I would like to do.


texts: Placemark name not shown after converting from .shp to kml using OGR

texts:
I've converted a bunch of .shp files to .kml using ogr2ogr and I've noticed that the placemark names in Google Earth are not reflected in the newly made kml files. I have also used ogr2gui and found that the placemark names do appear when made with them. With that said, ogr2ogr seems to be the better tool for automatic batch conversion of the files, and hence my desire to stick with it.

Edit: I have been told by Inventis (developers of Ogr2Gui) that by default the shp's column headers are replace by NameField and Description. Then, the driver fills those two tags with the attributes of the DBF.

In this case, could this be a driver issue? should I be building LIBKML and using that instead of the default KML?

I cannot use the 'NameField=' and "Description=" since the header names vary with each shp file.

Can someone help?


texts: QGIS Python Standalone Script using qgis.rastercalculator

texts:
I have been attempting to use the raster calculator algorithm in a standalone python script and have been trouble getting it to work. My issue is with the expression.

    precip = QgsRasterLayer(
    "./Rasters/Georef_Precip.tif", "Georef_Precip"
    )

    reclass_tmp = QgsRasterLayer(
    "./Rasters/Reclass_Temp.tif", "Reclass_Temp"
    )

    georef_rain = processing.run(
    "qgis:rastercalculator",
    {
        'EXPRESSION':f'{reclass_tmp}@1 * {precip}@1',
        'LAYERS':[precip],
        'CELLSIZE':0,
        'EXTENT':'-15052048.99614927,-12097374.584049538,5572589.970826002,7733078.012284913 [EPSG:3857]',
        'CRS':QgsCoordinateReferenceSystem('EPSG:3857'),
        'OUTPUT':'./Rasters/New_precip.tif'
    }
    )

Using the above expression gives an error of "raise QgsProcessingException(msg)
_core.QgsProcessingException: There were errors executing the algorithm."

    test = "./Rasters/Reclass_Temp.tif"
    precip_test = "./Rasters/Georef_Precip.tif"

    georef_rain = processing.run(
    "qgis:rastercalculator",
    {
        'EXPRESSION':f'{test}@1 * {precip_test}@1',
        'LAYERS':[precip],
        'CELLSIZE':0,
        'EXTENT':'-15052048.99614927,-12097374.584049538,5572589.970826002,7733078.012284913 [EPSG:3857]',
        'CRS':QgsCoordinateReferenceSystem('EPSG:3857'),
        'OUTPUT':'./Rasters/New_precip.tif'
    }
    )

And, if I use this second expression it runs successfully but the outputted raster is empty.

Can you help fix my error?


texts: How to get JSON representation from JTS Geometry?

texts:
I'm using JTS (http://www.vividsolutions.com/jts/JTSHome.htm) and need to get a JSON representation of a Point :

Hibernate beam:

@Type(type = "org.hibernate.spatial.GeometryType")
@Column(name="the_geom")
private Point theGeom;

Table:

CREATE TABLE cities
(
  id_city serial NOT NULL,
  name text,
  the_geom geometry,
  id_user integer,
  CONSTRAINT cities_pkey PRIMARY KEY (id_city),
  CONSTRAINT fk2de48dcf3c0b0090 FOREIGN KEY (id_user)
      REFERENCES users (id_user) MATCH SIMPLE
      ON UPDATE NO ACTION ON DELETE NO ACTION
)

I need something like:

jsonString = someClass.convertToJson( myBeam.getTheGeom() )


texts: Saving output from vector-tools as temporary file or straight into input layer?

texts: In mapinfo it is possible to save an output without creating new shp. Either by saving into the input- layer or to a temporary file. For example if you want to ”union” two layers you may want to get the output into one of the input layers. Anyone know of smiliar possibilities in qgis?


texts: Group and union polygons that share a border in R

texts:
I am attempting to union all polygons that: 1) share any portion of a border, or 2) have within a certain threshold distance (e.g., 50m) of each other, within a SpatialPolygons object.

for_cov_sp_0
class       : SpatialPolygons
features    : 1815
extent      : 313966, 751096, 462788, 781598  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=50 +datum=WGS84 +units=m +no_defs +ellps=WGS84
+towgs84=0,0,0

For example polygon 1 is within 50m of polygon 2 and polygon 2 is within 50m of polygon 3, so I would want to union all 3 together.

I have attempted using the poly2nb() function from the 'spdep' package with a snap distance to account for the threshold distance described above. I have been able to obtain a data frame of the following structure: each column and column is the position (#) of a polygon from my original SpatialPolygons object. The value at each cell [polygon #, polygon #] is a 1 if these two polygons are neighbors (share a border or borders are within 50m) and 0 if these two polygons do not.

nb <- poly2nb(for_cov_sp_0, snap = 50)
nb_df <- as.data.frame(nb2mat(nb, style="B", zero.policy=TRUE))
colnames(nb_df) <- as.character(seq(1, length(for_cov_sp_0)))
rownames(nb_df) <- as.character(seq(1, length(for_cov_sp_0)))

nb_df[1:10, 1:10]
   1 2 3 4 5 6 7 8 9 10
1  0 1 0 0 0 0 0 0 0  0
2  1 0 0 0 0 0 0 0 0  0
3  0 0 0 0 0 0 0 0 0  0
4  0 0 0 0 0 0 0 0 0  0
5  0 0 0 0 0 0 0 0 0  0
6  0 0 0 0 0 0 0 0 0  0
7  0 0 0 0 0 0 0 1 0  0
8  0 0 0 0 0 0 1 0 0  0
9  0 0 0 0 0 0 0 0 0  1
10 0 0 0 0 0 0 0 0 1  0

The final steps that I cannot figure out are: 1) some way to determine ALL polygons that are neighbors (e.g., polygon 6 and 7 are neighbors, and polygon 7 and 8 are neighbors, so polygons 6, 7, and 8 should all belong to the same group, and 2) create an ID for each group of polygons that are neighbors on which to apply a spatial union.

UPDATE: Roger Bivand suggested that I explore the 'sf' package, and an issue post on this package's GitHub repository explored a similar problem and answered the first portion of my question: https://github.com/r-spatial/sf/issues/234

However, I still cannot determine an automated way to union the polygons that are listed as neighbors.

tst <- for_cov_sf_3 %>% mutate(nbs = st_intersects(.))
tst
Simple feature collection with 3922 features and 1 field
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 319066 ymin: 460238 xmax: 737926 ymax: 775118
epsg (SRID):    32650
proj4string:    +proj=utm +zone=50 +datum=WGS84 +units=m +no_defs
First 20 features:
      nbs                       geometry
1       1 POLYGON((516766 775118, 516...
2       2 POLYGON((515686 775058, 515...
3       3 POLYGON((514516 774938, 514...
4       4 POLYGON((514336 774878, 514...
5       5 POLYGON((510916 773318, 510...
6       6 POLYGON((510166 772928, 510...
7       7 POLYGON((512806 772688, 512...
8       8 POLYGON((514426 772268, 514...
9       9 POLYGON((513346 772208, 513...
10 10, 11 POLYGON((514486 772028, 514...
11 10, 11 POLYGON((514426 771968, 514...
12     12 POLYGON((520966 769688, 521...
13     13 POLYGON((507166 768428, 507...
14     14 POLYGON((507346 768368, 507...
15     15 POLYGON((507436 768308, 507...
16     16 POLYGON((508006 768248, 508...
17     17 POLYGON((507886 768188, 507...
18     18 POLYGON((509296 768128, 509...
19     19 POLYGON((523306 768128, 523...
20     20 POLYGON((504826 768008, 504...


texts:
I'm diving into using the HTML popup in ArcGIS 10 and I'm trying to figure out how I can get an image to appear. I have a point file that contains the field "Coll" that has a image name with no extension ("100015"). All the images are stored in a directory ("T:\SEFIS2010_GVdata\SEFIS Video"). I can make an image appear in the popup using this XSL code

<xsl:when test="FieldName[starts-with(., 'Coll')]">
    <img src="T:\SEFIS2010_GVdata\SEFIS Video\100015.jpg"/>
</xsl:when>

but what is the syntax to append the directory, the FieldName, and the ".jpg" extension?


texts: Free map service for private web applications

texts:
Background

Looking to display a map (for internal use only) on a web page. The application must be able to add custom layers (without sending the data to the map service host):

Municipal zones (defined by shapefiles)
Custom zones (defined by polygons, possibly in terms of Municipal zones)
Site locations (defined by lat/long)
Additional city data (defined by shape files)

Because it is a private venture, it eliminates many map services, including: Google Maps, Bing, Yahoo! Maps, MapQuest, and other pay-for-use services.

Update #1

Perhaps MapQuest is a possibility.

Update #2

Requirements that I forgot to mention:

There are additional data (cities, towns, hamlets) not part of OpenStreetMap. These cities must be included without giving the data to the map provider.
There are various layers (such as employee locations) that also must be displayed, without telling the map provider where to find the data (i.e., the data will be behind a firewall).

Current Development

Currently, I have not been able to find a server that meets my requirements, so I have been following the instructions across various web sites to install a Mapnik web server:

http://dbsgeo.com/foss4g2010/html/_static/pdf/osm-rendering-foss4g2010.pdf
http://osgeo.org.cn/l18n/geodjango/install.html
http://wiki.openstreetmap.org/wiki/Mapnik/Installation
http://www.postgresonline.com/journal/archives/180-pgAdmin113plugins_postgis.html
http://blog.geoserver.org/2009/01/30/geoserver-and-openstreetmap/
http://workshops.opengeo.org/stack-intro/postgis.html
Bulk load multiple shapefiles into PostGIS

At this point, I have various pieces of the data:

City and municipality data (shape files loaded into PostGIS)
Natural Earth data (shape files; not loaded)
Mapnik stylesheets (various files; not loaded)

GeoServer

I tried installing GeoServer (prior to Mapnik) but found the documentation to be lacking, the web-based administration tool to be cumbersome to configure, and the responsiveness on a 2 GHz machine to be painfully slow (10 seconds to render a tile with a single layer: roads), and the default aesthetic bland. (I know I can use stylesheets to liven up the presentation, but that means having to find [or create!] suitable files, install them, test them, and so forth.)

Requirements

The technical requirements include:

Cities, roadways, bodies of water (lakes, rivers, oceans)
Political boundaries (provinces, countries)
High quality (visually similar to Google Maps, MapQuest, etc.)
High uptime
900913 projection (similar to Google Maps)
Responsive tiles (under 0.1 seconds per tile)
Excellent coverage of street data for Canada
Useful for slippy map (compatible with OpenLayers)
10,000 tile requests per day
Free for commercial (or private) use, or less than $2000 per year?

Good Example

Here is an excellent example (OpenLayers) of what I'm looking to achieve:

Bad Examples

An example of what I'm not looking to achieve:

An example of incomplete data:

Related Links

http://atlas.nrcan.gc.ca/site/english/dataservices/web_map_service.html
http://columbo.nrlssc.navy.mil/ogcwms/servlet/WMSServlet
http://wiki.openstreetmap.org/wiki/WMS
http://geoserver.org/display/GEOS/Available+WMS+and+WFS+servers
http://www.skylab-mobilesystems.com/en/wms_serverlist.html
http://www.ogcnetwork.net/servicelist
http://ows-search-engine.appspot.com/ - [broken]

Question

What WMS, WMTS, or WFS servers offer the required functionality?


texts: How can I dissolve adjacent polygons based on minimum attribute value?

texts:
I'm trying to merge the polygons (administrative boundaries) in a dataset untill they reach a minimum value. Lets say that the polygon dataset consists of the number of people per polygon and I need each polygon to contain at least 100 people. In that case I want the polygons with a value less than 100 to merge with an adjecent polygon, preferably in this order:

merge polygon with one adjecent polygon that contains less than 100 people
if adjecent polygons of less than 100 people are not available, merge with a polygon that contains 100 or more people.
if after the merge the minimum number is not reached yet, repeat the process untill no more adjecent polygons are available.
The end result should be a polygon dataset in which as many polygons as possible have at least 100 people per polygon.

I hope this discription is clear enough. Let me know if you need more info.

I'm using ArcMap 10.2.1. I already tried the dissolve tool, but this doesn't have the option to dissolve based on a minimum value. I saw that during an edit session in ArcMap you can merge polygons manually, but I found no way to do this in a model or a python script. If necessary I'm open to use other tools such as QGIS or SpatiaLite.


texts:
I am trying to calculate the cost of network to homes. In the image i have cable (brown line) between utility boxes(pink/blue points) and home feed (black dash) to the home.
I know the cost for each element (line, point) but i want to attach this cost to the homes it goes on to feed (tree branch)or attached the home ID to the features Is what im describing doable?


texts: Generating multiple visibility bearings in ArcGIS as automated procedure?

texts:
As a relative newbie to GIS, I'm quite sure there will be a way to do this in GIS (probably most efficiently by writing a python script!) but (to me anyway) this is a fairly complex thing to do and I have no idea where to start.  I know I am going to have to learn python and I'm happy to do this if it's going to allow me to achieve what I need. Here's what I'm trying to achieve (and this is not a one off, it will be a regular process, so if I can spend a bit more time at the beginning to automate as much of this as possible so much the better):

I have a shapefile of points, which are viewpoints in a particular area.

I then have multiple shapefiles (could be up to 50+) with other sets of points in each file (proposed development points,) each of which I want to buffer at (potentially different) distances (stored in the attribute table for each points file) to show different potential development extent for each.

What I'd then like to get out of GIS (I think using terrain / viewsheds info,) are the bearings from each point in the viewpoints shapefile, to the visible extents (left to right, i.e. 270-276 degrees) of EACH of my 50+ shapefile buffer files (note, not each point within each shapefile, just furthest left and furthest right from each viewpoint if possible?) The bearings would need to tangent the outer edges of the outermost buffers to show full potential extents, rather than just to the points themselves.

Also, to make matters even more complicated, I only want to give bearings of what can actually be seen? So if the furthest left point in one particular file is outwith the viewshed for that shapefile (I would need to run viewsheds for each of my 50+ files,) then the visibility bearings wouldn't include this, they would just go from the points that can be seen.

Finally, I'd need GIS to give me distances from each viewpoint to the nearest point of each development, so that would be the nearest point in each of my 50+ development shapefiles.

I have ArcGIS 10 with 3D Analyst extension.

Can anyone help? Does anyone have ANY idea if this is something that is realistic to achieve in GIS, or any tips in terms of how I can do this?

[EDIT]: Ok, so I realise that this is a complicated procedure, (though I'm pleased it seems to be doable!) so I'm splitting it into indivudual tasks, some I know how to do and some I'm sticking, but I think some of the sticking points are individual queries, so any additional questions while I work this out, I will link below...

Obtaining bearings to polygon extents in ArcGIS?

I'll also add some diagrams and further breakdowns as requested.


texts: MODIS DAAC2Disk fails with "Unable to obtain ECHO token"

texts:
I have been using the tool DAAC2Disk to easily download NDVI using an extent. However recently the tool fails with the error: "Unable to obtain ECHO token"

Has anyone encountered such error? Or know where to ask about it?


texts:
I'm new to QGIS, and have never used the Python console.
I have created a complex layer in QGis using an experimental, plugin which I had to manually install. It contains 91506 features; lines going around the country. I'm trying to run Q-chainage plugin to put points along these lines at regular intervals. I get the following error:

An error has occurred while executing Python code:
AttributeError: 'NoneType' object has no attribute 'asPoint'
Traceback (most recent call last):
  File "C:/Users/j6525/.qgis2/python/plugins\qchainage\qchainagedialog.py", line 124, in accept
    decimal)
  File "C:/Users/j6525/.qgis2/python/plugins\qchainage\chainagetool.py", line 199, in points_along_line
    divide)
  File "C:/Users/j6525/.qgis2/python/plugins\qchainage\chainagetool.py", line 69, in create_points_at
    point = QgsGeometry.fromPoint(point.asPoint())
AttributeError: 'NoneType' object has no attribute 'asPoint'

Python version: 2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]
QGIS version: 2.18.14 Las Palmas, ff83b9a479
Python Path:
•   C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins\processing
•   C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python
•   C:/Users/j6525/.qgis2/python
•   C:/Users/j6525/.qgis2/python/plugins
•   C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins
•   C:\PROGRA~1\QGIS2~1.18\bin\python27.zip
•   C:\PROGRA~1\QGIS2~1.18\apps\Python27\DLLs
•   C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib
•   C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\plat-win
•   C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\lib-tk
•   C:\PROGRA~1\QGIS2~1.18\bin
•   C:\PROGRA~1\QGIS2~1.18\apps\Python27
•   C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages
•   C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\jinja2-2.7.2-py2.7.egg
•   C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\markupsafe-0.23-py2.7-win-amd64.egg
•   C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\win32
•   C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\win32\lib
•   C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\Pythonwin
•   C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\Shapely-1.2.18-py2.7-win-amd64.egg
•   C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\wx-2.8-msw-unicode
•   C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\xlrd-0.9.2-py2.7.egg
•   C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\xlwt-0.7.5-py2.7.egg
•   C:/Users/j6525/.qgis2//python
C:/Users/j6525/Desktop/GIS assessment/02_fields_of_conflict/02_fields_of_conflict


texts: How to prevent Qgis clipper to clamp nodata values?

texts:
I'm trying to crop some old georeferenced cartography in order to have a continuous map (every map has its blank borders that overlaps).
Usign raster->extract->clipper and setting no value to -32768 as the original maps, I got the same message:

for band 1, destination nodata value has been clamped to 0, the original value being out of range.

the problem is all the blacks in the map are translated to nodata. If I set no data value to something greater than 255, it clamps to 255 and all the whites are hidden. Is there a way to fix that?
the commandline is like that:

gdalwarp -dstnodata -32768 -q -cutline "borders.shp" -crop_to_cutline -of GTiff "source.tif" "destination.tif"


texts: vb.net arcmap add-in log4net log still follow the configuration before update

texts:
I met a strange thing in my log4net for vb.net arcmap add-in.

At first, I follow the link http://www.codeproject.com/Articles/42276/How-to-Use-Log-Net-with-VB-NET-A-Simple-Step-By-St to write logs for my arcmap add-in project. There are 2 files related to the log file:

(1) Log4NetAssembly1.exe.log4net, the configuration for log4net is:

<log4net>
<appender name="ExampleAppender" type="log4net.Appender.RollingFileAppender">
  <file value="C:\Log4NetExamples\Example.log" />
  <appendToFile value="false" />
  <rollingStyle value="Size" />
  <maxSizeRollBackups value="3" />
  <maximumFileSize value="100KB" />
  <staticLogFileName value="true" />
  <layout type="log4net.Layout.PatternLayout">
    <conversionPattern value="%date [%thread] %-5level -
            [%logger] %message%newline" />
  </layout>
</appender>
<logger name="dev">
  <level value="ALL" />
  <appender-ref ref="ExampleAppender" />
</logger>

(2) app.config

<listeners>
    <add
        name="textWriterTraceListener"
        type="System.Diagnostics.TextWriterTraceListener"
        initializeData="C:\Log4NetExamples\Example.log" />
  </listeners>

It works properly with log file C:\Log4NetExamples\Example.log.

Now I want to create log file daily and want each log file with the date of creation, so I updated my code as follows:

(1) Log4NetAssembly1.exe.log4net

<log4net>
<appender name="RollingLogFileAppenderOutput" type="log4net.Appender.RollingFileAppender">
  <file value="C:\Log4NetExamples\" />
  <appendToFile value="true" />
  <rollingStyle value="Date" />
  <datePattern value="yyyyMMdd'.log'" />
  <staticLogFileName value="true" />
  <layout type="log4net.Layout.PatternLayout">
    <conversionPattern value="%date [%thread] %-5level -
            [%logger] %message%newline" />
  </layout>
</appender>
<logger name="dev">
  <level value="ALL" />
  <appender-ref ref="RollingLogFileAppenderOutput" />
</logger>

(2) app.config

<listeners>
    <add
        name="textWriterTraceListener"
        type="System.Diagnostics.TextWriterTraceListener"
        initializeData="C:\Log4NetExamples" />
  </listeners>

The strange thing is after I changed to the daily logging design, my log file still written to C:\Log4NetExamples\Example.log, which means it still follow the previous configuration of log4net. I really cannot figure out why, is there anywhere remember the old configuration? I search the whole vb.net project and cannot find "Example.log" at all.

Anybody knows what happen? Please help me, thanks.


texts: Creating Color Ramp programmatically using ArcPy in ArcGIS Pro?

texts:
I am aware of accessing color ramp using the ColorRamp method of the arcpy.mp module. However, this only gets existing color ramps.

Is there any way of creating new color ramp programmatically, for example, Yellow to Black?


texts: How to count occurrences of one field grouped by values of another?

texts:
Is there a way of calculating the number of unique text strings in a field in ArcGIS 10.1? I have a number of repeated strings in this field, but would like to calculate the number of unique strings in a field based on groupings from another field.

E.g  - Current table:

Field 1    Field 2
Apple      Green
Apple      Red
Apple      Red
Citrus     Yellow
Citrus     Green
Citrus     Orange
Citrus     Orange

Result:

Field 1    Field 2
Apple      2
Citrus     3


texts: Size of the output files cannot match after converted from polygon to raster

texts:
Basically, the emodisfile is the raw data file, and shapefile is one of the CLU shapefiles. I need to convert the shapefile from polygon to raster and then to ASCII file based on the emodisfile.
So I set the extent, cell size and coordinate system to be exactly the same as the emodisfile. After that, I convert the shapefile and the output ASCII file is shpfipstest.TXT.

And I also convert the emodisfile to ASCII file, which is called emodisCoordtest.txt.
Since I convert the shapefile all based on the settings of emodisfile, then they should have the same extent and Coordinate System.

However, according to the results, they cannot match together and have different number of columns and rows.

While I also tested it with ArcMap, it had the expected result.
But when it comes to Python and ArcPy, it seems something wrong.

Any ideas? Below is a part of the code:

emodisfile = workingDir+"\\US_eMTH_NDVI.2014.365-006.QKM.VI_NDVI.005.2014010002626.tif"
shapefile = workingDir+"\\clu_public_a_ar039\\clu_public_a_ar039.shp"
env.workspace = workingDir
env.overwriteOutput = True
arcpy.env.coincidentPoints="MAX"
arcpy.env.extent = emodisfile
arcpy.env.cellSize = emodisfile
arcpy.env.cartographicCoordinateSystem = emodisfile
outraster = workingDir+"\shapeFIPS"
outascii = workingDir+"\shpfipstest.TXT"
arcpy.PolygonToRaster_conversion(in_features=shapefile, value_field="FIPS",
   out_rasterdataset=outraster, cell_assignment="CELL_CENTER")
arcpy.RasterToASCII_conversion(in_raster=outraster, out_ascii_file=outascii)
emodisAscii = workingdir+"emodisCoordtest.txt"
arcpy.RasterToASCII_conversion(
    in_raster=emodisfile, out_ascii_file=emodisAscii)

I have changed some settings into:
des = arcpy.Describe(emodisfile)
arcpy.env.extent = des.extent
arcpy.env.cellSize = des.children[0].meanCellHeight

And right now the files have the same numbers of rows and columns, and origin coordinates, but a problem came up that the outraster file did not have an attribute table, while it should contain "FIPS" value in it.
Therefore, in the outascii file, all values are -9999.


texts: How to stop QGIS from looking for the wrong PostgreSQL host address in an existing workproject?

texts:
It looks like a saved QGIS project workfile may be storing the old PostgreSQL host IP access details. Updates made to the PostgreSQL host IP in the main QGIS window are not being used in the project workfile.  The result is that when opening my project PostgreSQL looks for the wrong host and then comes up with a "time-out" and is stuck in a loop...

I attempted to set up a new PostgreSQL file with a revised host address etc and tested this successfully before accessing my project workfile.  However, I still get a time out error with PostgreSQL, as within the project workfile the pointer is still set at the old host.


texts: OGR UnionCascaded gives strange lines (flat polygons)

texts:
I'm using OGR to perform a union using UnionCascaded. All my polygons have been created by GDAL polygonize which ensures that all polygons follow pixel edges and should therefore align nicely. But I get these strange flat polygons looking like lines:

They are positioned at locations where two polygons met before they where unioned.

Can you recognise this pattern and what lead to it?


texts: Python script not delivering expected output for Vegetation indices (NDVI, CLgreen) using ArcPy and ArcGIS

texts:
I have a script that was written in Python in PyCharm.

The script is something meant to run once then all the indices I list in the python can run from PyCharm rather than opening the ArcMap or ArcGIS Pro software everytime.

When running the script I get the NDVI that is expected, however the CLG index does not provide a reading as expected.

The Clgreen using the Arcpy provides an output, while the Clgreen2 using the ArcGIS fails in giving an output.

import arcpy
import arcgis
from datetime import date

output = r'C:\Users\john\Desktop\testing\RedEdge\Exports'
raster = r'C:\Users\john\Desktop\testing\RedEdge\Exports\RedEdge.tif'

today = date.today()
date = today.strftime("%d-%m-%Y")

#formulas - indices
NDVI = arcpy.ia.NDVI(raster, nir_band_id=5, red_band_id=3)
Clgreen = arcpy.ia.CIg(raster,nir_band_id=5, green_band_id=2)
Clgreen2 = arcgis.raster.cig(raster,band_indexes="5 2", astype=None)

#save
Clgreen.save(output+'\\'+date+'_' +'clgreen.tif')
Clgreen2.save(output+'\\'+date+'_' +'clgreen2.tif')
NDVI.save(output+'\\'+date+'_'+'ndvi2.tif')

print('scripting was successfull ! ')

Image 1 is Clgreen in imageviewer and image 2 is in ArcGIS PRO.

So the questions I have is why is the output of the Clgreen not shown
in ArcMap or ArcGIS pro, but shows fine in imageviewer
What would the reason be that the Clgreen2 is failing?


texts: How to protect my application ESRI FGDB data from authorized access using other applications?

texts:
I have ArcObjects app which uses file geodatabase for storing data and performing CRUD operations. This file geodatabase  needs to be protected from reading data in other applications such as ArcMap, Civil 3d and other 3rd party's.
Is there any way to do this?

Thanks


texts: Joining table with multiple-values-field in QGIS

texts:
I have a vegetation vector layer with vegetation type mapping. I want to join this vector layer with a CSV-file (no geometry) which includes threatened species occurrence.

The vector layer has a "vegetation class" field which provides a single vegetation classification for each polygon. The CSV-file has a species list with a "Class" field which indicates all the vegetation classes (comma separated) the species is known to occur.

When I create a join between the vector layer and the CSV-table data, there doesn't appear to be anywhere for me to indicate the join field contains multiple values and I, therefore, receive NULL results when inspecting the join.

How can I indicate that the value from the target field ("vegetation class"), which is a single value, needs to search each row in the join field ("Class") which contains multiple comma separated values?


texts:
i want to display raster images using a color palette in a leaflet map. The images are tiled (gdal2tiles) from a single .tif to multiple .png files with values from 0 to 100.

I add the image using the tileLayer function.

Is there a way to use a color palette in leaflet to display the values (e.g: 0=blue, 100=green)?

if not

is there a way to transform grayscale .png (or the input .tif) to a RGBA png using a predefined color palette?


texts: Geoserver ECQL Intersects returns empty result

texts:
I am trying to understand ECQL and the intersects filter.
This filter entered on the geoserver layer preview returns my WMS as ecpected

Intersects(the_geom,POINT(1916881.37823 5762168.62373))

Now I try to use this filter on a URL to retrieve the WFS for the feature:

http://myserver/fwsys/wms?service=WFS&version=1.0&request=GetFeature&typeName=fwsys:fwsys_region&CQL_FILTER=Intersects(the_geom,POINT(1916881.37823%205762168.62373))

The browser returns this result

<wfs:FeatureCollection xsi:schemaLocation="http://fire.niwa.co.nz/fwsys_region http://gs-test.niwa.co.nz/fwsys/wfs?service=WFS&version=1.0.0&request=DescribeFeatureType&typeName=fwsys%3Afwsys_region http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.0.0/WFS-basic.xsd">
    <gml:boundedBy>
        <gml:null>unknown</gml:null>
    </gml:boundedBy>
</wfs:FeatureCollection>

Instead of my expected WFS...

I have tried the same URL but with a polygon filter and it worked.

What am I doing wrong?

Updated this questions wuth details about the layer.

The layer I am trying to use is in EPSG:2193

PROJCS["NZGD2000 / New Zealand Transverse Mercator 2000",
    GEOGCS["NZGD2000",
        DATUM["New Zealand Geodetic Datum 2000",
            SPHEROID["GRS 1980", 6378137.0, 298.257222101, AUTHORITY["EPSG","7019"]],
            TOWGS84[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
            AUTHORITY["EPSG","6167"]],
        PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],
        UNIT["degree", 0.017453292519943295],
        AXIS["Geodetic longitude", EAST],
        AXIS["Geodetic latitude", NORTH],
        AUTHORITY["EPSG","4167"]],
    PROJECTION["Transverse_Mercator", AUTHORITY["EPSG","9807"]],
    PARAMETER["central_meridian", 173.0],
    PARAMETER["latitude_of_origin", 0.0],
    PARAMETER["scale_factor", 0.9996],
    PARAMETER["false_easting", 1600000.0],
    PARAMETER["false_northing", 10000000.0],
    UNIT["m", 1.0],
    AXIS["Easting", EAST],
    AXIS["Northing", NORTH],
    AUTHORITY["EPSG","2193"]]


texts:
I'm working with Geoserver 2.12.1 (Windows Installer Version) and I had to enable CORS before for getting access when I attempted to do a GetFeatureInfo.

I did the steps to enable CORS:

Edit webapps/geoserver/WEB-INF/web.xml:

<filter>
    <filter-name>cross-origin</filter-name>
    <filter-class>org.eclipse.jetty.servlets.CrossOriginFilter</filter-class>
  <init-param>
      <param-name>allowedOrigins</param-name>
      <param-value>*</param-value>
  </init-param>
  <init-param>
      <param-name>allowedMethods</param-name>
      <param-value>GET,POST,DELETE,PUT,HEAD,OPTIONS</param-value>
  </init-param>
  <init-param>
      <param-name>allowedHeaders</param-name>
      <param-value>origin, content-type, accept</param-value>
  </init-param>
</filter>

<!--
  THIS FILTER MUST BE THE FIRST ONE, otherwise we end up with ruined chars in the input from the GUI
  See the "Note" in the Tomcat character encoding guide:
  http://wiki.apache.org/tomcat/FAQ/CharacterEncoding
-->
<filter-mapping>
  <filter-name>Set Character Encoding</filter-name>
  <url-pattern>/*</url-pattern>
</filter-mapping>

Now, the GetFeatureInfo for WMS Service is working, but, I recently installed WPS Plugin and when I do some request to:

http://localhost:8080/geoserver/wps

It's not working. It's showing the following error:

Failed to load resource: the server responded with a status of 403 (Forbidden)
(index):1 Failed to load http://localhost:8080/geoserver/wps: Response to preflight request doesn't pass access control check: No 'Access-Control-Allow-Origin' header is present on the requested resource. Origin 'http://localhost:3000' is therefore not allowed access. The response had HTTP status code 403.

But I have it configured to accept all methods.

What could be the reason?


texts: QGIS Data defined labels within a single layer - prioritize display

texts:
I am creating a world map where I want to make sure that certain countries with an incident are ALWAYS labelled automatically but the other countries are labelled ONLY when there is space to fit the label.

I have created a display field in the table and populated it with

10 - countries where the name MUST be displayed
1 - other countries where the label should show if it doesn't overlap.

Based on QGIS label priority settings I have created a statement

I understand that this won't work but I need to be able to say

Case when
 "Display"=10 then "CNTRY_NAME" - must display
 "Display"=1 then "CNTRYY_NAME" - should display when not overlapping

How do I achieve this in v3.12 without having to duplicate the layer and creating a filter where the labels are shown even when overlapping (filter display=10) and do not overlap (filter display=1)


texts:
I have created many rasters in Google Earth Engine that represent the NDVI. I want to apply statistics on them. I want to count how many pixels I have in every raster, and then I want to count how many of the pixels have values higher than 0.3 and 0.5.

This is the code I wrote, I have tried to use the ee.Reducer.count(), but I get error messege- ''Invalid argument specified for ee.List(): EPSG:4326'' and also I don't know how to make it count values that are bigger than 0.5.

var agri_7=agri_1;
//2000-2002
var dataset = ee.ImageCollection('MODIS/006/MOD13A1')
.filter(ee.Filter.date('2000-03-01', '2002-03-31'))
.select('NDVI');
// filter the image collection to contain only images from March
var onlyMarch = dataset.filter(ee.Filter.calendarRange(3, 3,'month'));
print(onlyMarch);
// make a composite (mean image) of the images of March and clip to the geometry
var meanImage = onlyMarch.mean().clip(agri_7).divide(10000);
print(meanImage);

var mean20002002 = meanImage.reduceRegion({
  reducer: ee.Reducer.mean(),
  geometry: agri_7,
  scale:500
});

//var ndvi = dataset.select('NDVI')
var ndviParams = {min: 0, max: 1, palette: [
    'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
    '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
    '012E01', '011D01', '011301'
  ]};

Map.centerObject(agri_7, 9);
Map.addLayer(meanImage,ndviParams,'NDVI20002002');
print(mean20002002,'mean NDVI 2000-2002');

//count the number of total pixels
var c = meanImage.reduceRegion({
  reducer: ee.Reducer.count(),
  geometry: agri_7,
  crs: 'EPSG:4326',
  crsTransform: 'EPSG:4326',
  });

print(c,'Total Pixels 2000 2002');


texts: Is any spatial utility or methods are available for creating contour or isopleths on Oracle spatial? I want to create contours from point geometry and show it on a map (using mapserver) of the site. Please share me any idea/information related to this.


texts: OpenDataCube OWS issue exposing Sentinel 5 tile

texts:
I am trying to index and expose a Sentinel 5 file (netcdf) with OpendataCube OWS (gdalinfo and ncdump available at the end of this post). The Sentinel-5 file is in the NetCDF4 format and uses hierarchy groups to store variables.
Adding a product and a dataset seems to work fine (no exception appears in the console).

I also modified the ows_cfg.py to expose the dataset (attached to this post). Note: custom layer folder at line 1040.
Unfortunately when I execute 'datacube-ows-update --views' I encounter the following error:
psycopg2.errors.InternalError_: parse error - invalid geometry                                                                                HINT:  "POLYGON((  , " <-- parse error at position 13 within geometry

Attachments available at: https://drive.google.com/drive/folders/1yiRQfQw_Yzu17Km6B6El08sxEVfBZmVY?usp=sharing


texts: Setting GeoServer's WMS services SRS limitation via REST Api

texts:
On the GeoServer admin page, it is possible to restrict the Spatial Reference Systems (SRS) that can be requested in WMS GetMap requests by putting a comma separated list of SRS numbers in the input box at
http://localhost:8080/geoserver/web/?wicket:bookmarkablePage=:org.geoserver.wms.web.WMSAdminPage
(set server and port according to your server).

In consequence, only the specified SRS numbers can be used for the SRS parameter in GetMap requests (all others will not get a valid response). Besides, limiting GeoServer to only a few SRS that are really needed by your application will shorten the GetCapabilities document significantly because a very long list of all available SRS is included by default when there is no limit set.

Is there any possibility to set this limit via GeoServer's REST Api? Unfortunately, I could not find any example on this in the documentation.


texts: How gdal_rasterize works - how it interpolates / resample values?

texts:
I have huge number of points saved as shapefile.

I created an empty raster file with low spatial resolution. Next, I applied gdal_rasterize on mentioned .shp and used low resolution raster as destination.

It works, but since I didn't provided an resampling method, I wonder which algorithm has been applied to prepare raster for me.
Unfortunately, I can't find this information within the GDAL documentation: http://www.gdal.org/gdal_rasterize.html


texts: I have a pretty robust dataset of .gpx files i'm animating with the torque layer style. I'd like to be able to color the map by type using a column...for me its 'car' vs 'flight' vs 'train. This works with simple colors, but havent had luck with the torque and css.


texts: How to match air photo to existing vector layers on Mac OS?

texts:
I'm trying to add an air photo to a project and match it to existing data. CanVec data is mostly based on the old CIA dataset of the 60's and is hopelessly out of date, but the roads are up to date - the roads can easily be seen on the air photo. So, I need to use the roads to line the photo up. This will be an iterative process, of course.

I read of a Georeferencer plugin for QGIS, but can't persuade QGIS to find it (it doesn't seem to have the "3rd party" option shown by online tutorials). Will it do what I need it to? If so, how can I install it?


texts:
This is about the lower part of TesterNoType

In table 1 I have ex

Filename, Type, Date
Parks,    1,    2018-01-22
Parks,    4,    2018-01-22

Table 2 aka new information table I have information about 1000 parks
In table 3 aka old information table I have 988 parks, some have stayed the same and should not be changed, some have old info about the type of the park

I have another workspace where I send the filename as Feature to Read and reads the table and writes to table 2.

In the old one with no types I only wrote the reader from the new and truncated the old and voila, but here I can't truncate, I want to both add all new rows and for every row that has type 1 for example (same obj id can happen in old and new or not).

I only want to replace every single row in table 3 that has type 1 and replace it with the information from table 2 instead. All that not has type 1 should still be there.

I'm doing a very simplified version of this now where I only read Parks, 1
and should read parks table from the old GDB and write to another GDB where I have the old version of Parks, but the writing part where I only want to replace where attribute type=1 doesn't seem to be working.

Should I do a SQL script somewhere or only read rows where type=1 or find another way to do this?


texts: Code-completing/auto-completing ArcPy scripts using Notepad++ (external editor)?

texts:
I have been scripting in Python for ArcGIS 10, and so far I have just been using Notepad++.

I like how the Python command line interface code-completes arcpy functionality for you, showing you your options, and I was wondering if there is a way to add this feature to an external editor.


texts:
I have 2 JP2 aerial imagery files that are incorrectly sized. They are 100x larger than the other aerial images but correctly geo-referenced (starting x, y co-ordinates).

I believe the issue is their defined pixel size ( 12.5, -12.5 ) as the other aerial imagery files have pixel sizes of: ( 0.125, -0.125 ).

I have opened the JP2 files up in notepad and changed these parameters but it corrupts the file or saves it as text. Dos anyone know what Raster tools in QGIS I can use to remedy this or any other techniques?

Thanks for your time.


texts: Changing projection of geotiff file (raster data) in QGIS?

texts:
After downloading an SRTM geotiff (DEM) from USGS's website, I would like to convert it from GCS to a projected coordinate system. After opening it in QGIS:

Is that enough if I simply just click on 'save as...' and set CRS from WGS84 to the preferred PCS?
Or should it be done like 'Raster > Projections > Warp'?

(I feel like it is not that easy and there is a trick somewhere...)


texts: How to create Google map fixed marker with postal address

texts:
I use this script to show a Google fixed map It's marker in the center always

and it works just fine. The issue is how to add postal address that shows in the marker content?


texts: Statistics on classified raster overlaid with a polygon vector

texts: I have a classified raster file and a polygon vector file. I need to overlay polygon and find out area occupied under each classification in each polygon. I am not familiar with R. Can you please suggest a solution in QGIS?


texts: Cannot Install QGIS Mapserver On Ubuntu 12.04

texts:
I'm trying to install qgis-mapserver on my ubuntu 12.04. I just installed qgis desktop 2.6.0-brighton successfully. Here is my steps when installed the desktop :

sudo apt-get install python-software-properties
sudo add-apt-repository ppa:ubuntugis/ubuntugis-unstable
sudo apt-get update
sudo apt-get install qgis python-qgis qgis-plugin-grass

But when I tried to install the mapserver :

sudo apt-get install qgis-mapserver libapache2-mod-fcgid

I got this error :

The following packages have unmet dependencies:
qgis-mapserver : Depends: qgis-providers (= 2.4.0-0precise2) but 2.6.0-0precise1 is to be installed


texts:
In ArcMap, I can change the selection outline from the default RGB (0, 255, 255) to whatever color I want by going to Selection > Selection Options > Selection Tools Settings.  In the example below, I have set my selection color to a red outline.

I'm trying to do the same in QGIS, but the entire polygon is shaded red, instead of with just a red outline like I can do in ArcMap.
Is it possible to achieve similar selection display in QGIS as I can in ArcMap?


texts:
I have a TIFF file (Landsat8 tile) and a polygon of a burned area, I would like to plot only the polygon appropriately within the bounds of the TIFF and save that as an image however I am struggling.
The polygon has been uploaded as a GitHub gist here.
The image TIFF can be found here.
Plotting the TIFF with rasterio and matplotlib:
with rasterio.open(f'{input_dir}/{band2}', 'r', height=1024, width=1024) as t:
    f, ax = plt.subplots(figsize=(10,10))
    show(t.read(1), ax=ax)
    print(t.crs)

The CRS of this image is: EPSG:32610
Plotting the polygon produces the following:
f, ax = plt.subplots()
gpd.GeoSeries(boundary[0], crs='EPSG:32610').plot(ax=ax, color='red')

The CRS doesn't seem to be correct because this polygon is now way larger than the image it should be in the bounds of. When plotting the polygon directly onto the tile (just to see what that would look like) produces the following:
with rasterio.open(f'{input_dir}/{band2}', 'r', height=1024, width=1024) as t:
    f, ax = plt.subplots(figsize=(10,10))
    show(t.read(1), ax=ax)
    gpd.GeoSeries(boundary[0], crs='EPSG:32610').plot(ax=ax, color='red')

The end result I am hoping for is just the polygon where it should exist in the bounds of the tile on a blank background.
Edit
The original CRS of the polygon was EPSG:4326. However if I then try the following:
gpd.GeoSeries(boundary[0], crs='EPSG:4326').to_crs(epsg=32610)

I get the polygon as seen above.
0    POLYGON ((515747.044 4334707.851, 515863.285 4...
dtype: geometry

I've uploaded the original polygon in the EPSG:4326 here.


texts: Setting colour on point cloud data from raster

texts:
I was reading this article: 14 ways to take charge of lidar data.

Number 12 states that you can change the colour of your point-cloud data by using a raster image. Does anyone know the process for doing this within ArcGIS?

I am looking for a way to make the visualization of lidar data really pop and I have high resolution aerial photography that I could use.


texts:
I'm trying to create a globe view in QGIS 3.16.1 with an Azimuthal Orthographic projection. This is my projection codes:
+proj=ortho +lat_0=32.4279 +lon_0=53.688 +x_0=0 +y_0=0 +R=6371000 +units=m +no_defs +type=crs

After the new projection are defined and there was no error, after applying this new projection, my map does not appear spherical (globe veiw). What is the problem?


texts: How to use NASA Digital Elevation model data?

texts:
I've heard it is possible to use NASA elevation model to determine approximate elevation at (lat;lng) point.
This might be really useful for my upcoming project since Google Elevation API has weird limitations...and I don't really like USGS Elevation Query Web Service because of its speed. I'd want to have elevation query tool/framework on local computer, which works with local file.

I did my homework, and found GeoTiff files might be downloaded with no problem. But I've failed to find any useful documentation or framework.

Can somebody share a link to the specs, or example or tutorial on how to work with GeoTiff files?
Is there any frameworks to get this done (.NET,Java,C++,C - doesn't really matter)?

According to what I see in the GeoTiff files and its header files, I assume each tiff represents small area and pixel color (16 bit value) indicates elevation there in meters. Is that correct?


texts: Are there any ground truth data available for San Francisco?

texts:
I'm working on the algorithm in this paper and I've modified it a bit. Now in order to quantify the accuracy of my algorithm with this algorithm, I need ground truth data of SanFrancisco.

Which has at least five classes:

buildings (if detailed, it'll be better)
sea (if detailed better)
fir trees of the vegetated area
spruce trees of the vegetated area
pine trees of the vegetated areas
(Generally, it would be better to have a ground truth as detailed as possible).
Do you know such ground truth or any useful rasterized map to be able to be able to make a ground truth with optical imagery?

Maybe this question is a better fit for opendata.


texts: Opening template project with PyQGIS in QGIS 3?

texts:
What do I need to do to open a project-file as a template with PyQGIS in QGIS 3?

I am able to open a project like this:

project = QgsProject.instance()
project.read('C:/path/to/my/project.qgs')

But how do I open this "project.qgs" as a template-project?


texts: Run code on pressing a key in Pythonaddins

texts:
I tried making an addin having below code where I want the user to press a predefined key and the "required_code" should run on the selected features.

import arcpy
import pythonaddins

class Key1(object):
    """Implementation for Key Listener_addin.Key1 (Button)"""
    def __init__(self):
        self.enabled = True
        self.checked = False
    def onKeyUp(self, keycode, shift):
        print("Process Started")
        print("Key pressed: " + keycode)
        self.checked = True
        #required_code goes here with messagebox
        print("Process Completed")
        self.checked = False

However The above code doesn't do anything, not even printing first line.

How could I enable this addin to run the code on pressing a key?
Also how can I print which key is pressed by user?

Could you correct me on what I am doing wrong?

UPDATE:

I have tried to define keycode using msvcrt but it is still not working. My updated code is:

import arcpy
import pythonaddins
from msvcrt import getch

class Key1(object):
    """Implementation for Key_Listener_addin.Key1 (Button)"""
    def __init__(self):
        Key1.enabled = True
        Key1.checked = False
        keycode = msvcrt.getch()
    def onClick(self):
        if Key1.checked == True:
            Key1.checked = False
        elif Key1.checked == False:
            Key1.checked = True
    def onKeyUp(self, keycode, shift):
        if Key1.checked == True:
            if keycode == 85:
                print(keycode)

85 is the ASCII value for 'U'.


texts: Converting QGIS Algorithm from history to a plugin / tool

texts:
Rather than going through the process of creating a plugin is there a way to take algorithms from the history and put them together to create a quick tool like in ArcMap?

For example I have

processing.run("native:buffer", {'INPUT':QgsProcessingFeatureSourceDefinition('_DCDB_July2017_a5475545_fb58_4c0a_b257_588c1ea8c0d4', True),'DISTANCE':1000,'SEGMENTS':5,'END_CAP_STYLE':0,'JOIN_STYLE':0,'MITER_LIMIT':2,'DISSOLVE':False,'OUTPUT':'TEMPORARY_OUTPUT'})

and

processing.run("native:selectbylocation", {'INPUT':'D:/Cadastre/_DCDB_July2017.tab','PREDICATE':[0],'INTERSECT':'Polygon?crs=EPSG:28355&field=OBJECTID:integer&field=LOT:string(5)&field=PLAN:string(10)&field=LOTPLAN:string(15)&field=SEG_NUM:integer&field=PAR_NUM:integer&field=SEGPAR:integer&field=PAR_IND:integer&field=LOT_AREA:double&field=EXCL_AREA:double&field=LOT_VOLUME:double&field=SURV_IND:string(1)&field=TENURE:string(40)&field=PRC:integer&field=PARISH:string(20)&field=COUNTY:string(16)&field=LAC:integer&field=SHIRE_NAME:string(40)&field=FEAT_NAME:string(60)&field=ALIAS_NAME:string(254)&field=LOC:integer&field=LOCALITY:string(30)&field=PARCEL_TYP:string(24)&field=COVER_TYP:string(10)&field=ACC_CODE:string(40)&field=CA_AREA_SQM:double&field=O_SHAPE_AREA:double&field=O_SHAPE_LEN:double&field=Shape_Length:double&field=Shape_Area:double&uid={61b8a064-9170-42ea-b1bd-117e1d914dfb}','METHOD':0})

I want to be able to select a block, click on the button and enter in the buffer distance and it should add a layer with this.

In graphical modeller I can't get it to only buffer on a selection -otherwise this would have been a possible solution.

--- Update based on How to run QGIS algorithms with selected features? I tried the following

import processing
input_vlayer = 'D:\Cadastre\_DCDB_July2017.tab'
intersect_vlayer = 'Some path'

input_vlayer = iface.addVectorLayer(input_vlayer,  'input_vlayer', 'ogr')
intersect_vlayer = iface.addVectorLayer(intersect_vlayer,
                                       'intersect_vlayer', 'ogr')

for i in range(5):

    # Select the features
    input_vlayer.select(i)
    # Set input params
    params = {'INPUT':input_vlayer,'DISTANCE':1000,'SEGMENTS':5,'END_CAP_STYLE':0,'JOIN_STYLE':0,'MITER_LIMIT':2,'DISSOLVE':False,'OUTPUT':'TEMPORARY_OUTPUT'}
    result = processing.run("native:buffer", params)

But it runs on the whole dataset.


texts: What is a good way to retrieve geographic data along with images and mp3 files?

texts:
I'm working in the Django backend of an audio guide for tourism with mobile devices. Basically, the client (Phonegap app for Android and iOS) can follow a route over a map (Leaflet client) and see images or listen audio files in certain points of the route. Now, I want to provide offline service. So, the client can choose a route and download it, meaning download all the POIs, including images an audio files.

Of course, I can send the coordinates of the POIs in a json, or a gpx file. But what concerns me is the audio files and the images. Some random thoughts:

The devices can store the coordinates in SQLite, and the audio files and images in the SD card
Store the audio files and images as BLOBs in SQLite is not an option, because some files can be really big (like 100MB)
My first approach is to create a JSON with the coordinate points and the name of the files related with them, and compress it along with the mp3 and images. Then, send that zipped file to the clients. Sounds... gross.

Is there any recommended way to do this? Send geolocated audio and images, basically.


texts: Transfer data from access mdb to personal geodatabase (OLE, Python scripting)

texts:
I am very new in Python. And I guess I need to use lots of Python here.

I have MS Access 2003-2007 database (mdb format) and personal geodatabase (mdb). I created a tool in Access that imports data from personal geodatabase to access mdb, checks it and updates changes.

Now I am working on implementing backwards functionality (transfer data from access mdb to personal geodatabase). For now, I can implement that manually using OLE connection.

How I do it:

1.in ArcCatalog create OLE connection,

2.Connect to my access database backend,

3.Choose a Table I want to use from this access database backend,

4.Open ArcMap, add Feature class from personal geodatabase to the map,

5.Perform join of Feature class and Table,

6.Then using Field calculator update all necessary fields.

For now in my vision I can create custom tool for join Feature class and Table.

Do you have any suggestions how I can automate this process, at least half of this process.

Here is a code to join Table to Feature Class

import arcpy
from arcpy import env
env.overwriteOutput = 1
env.workspace = r"D:\Working\Tools\ToolShare\Scripts"

arcpy.MakeFeatureLayer_management("Crossings", "tempLayer")
arcpy.AddJoin_management("tempLayer", "UNIQUE_NAME", "tblCrossings", "UNIQUE_NAM")

When I open the feature class, there are no changes...
When I tried to use arcpy.JoinField_management("tempLayer", "UNIQUE_NAME", "tblCrossings", "UNIQUE_NAM")instead of arcpy.AddJoin_managementfields were added to the feature class, but permanently (I don't need that).

After this part I guess I need to use UpdateCursor (Field Calculator point). But I don't know how correctly it can be used here.

Here are my thoughts, please help me to organize them.

Join and calculate field script

import arcpy
import sys
import os
from arcpy import env

env.overwriteOutput = 1
scriptPath = sys.path[0]
arcpy.AddMessage("Script folder: " + scriptPath)
toolSharePath = os.path.dirname(scriptPath)
arcpy.AddMessage("ToolShare folder: " + toolSharePath)

inFeatureClass = arcpy.GetParameterAsText(0)
inTable = arcpy.GetParameterAsText(1)

arcpy.MakeFeatureLayer_management("Crossings", "tempLayer")

arcpy.AddJoin_management("tempLayer", "UNIQUE_NAME", "tblCrossings", "UNIQUE_NAM")

arcpy.CalculateField_management("tempLayer", "Status_", "!Status!", "PYTHON_9.3", "")

arcpy.RemoveJoin_management("tempLayer", "")

Add feature class script

import arcpy
from arcpy import env

inFeatureClass = arcpy.GetParameterAsText(0)
mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd,"*")[0]
newlayer = arcpy.mapping.Layer(inFeatureClass)
arcpy.mapping.AddLayer(df, newlayer,"BOTTOM")

Add feature class, join and calculate together

import arcpy
import sys
import os
from arcpy import env

env.overwriteOutput = 1

scriptPath = sys.path[0]
arcpy.AddMessage("Script folder: " + scriptPath)
toolSharePath = os.path.dirname(scriptPath)
arcpy.AddMessage("ToolShare folder: " + toolSharePath)

inFeatureClass = arcpy.GetParameterAsText(0)  # Feature class from personal geodatabase
inTable = arcpy.GetParameterAsText(1)  # Table from access database (mdb)

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd,"*")[0]
newlayer = arcpy.mapping.Layer(inFeatureClass)
arcpy.mapping.AddLayer(df, newlayer,"BOTTOM")

arcpy.MakeFeatureLayer_management("Crossings", "tempLayer")
arcpy.AddJoin_management("tempLayer", "UNIQUE_NAME", "tblCrossings", "UNIQUE_NAM")
arcpy.CalculateField_management("tempLayer", "Status_", "!Status!", "PYTHON_9.3", "")
arcpy.RemoveJoin_management("tempLayer", "")


texts: Using latitude,longitude of point and distance in degree to get another point?

texts: Suppose I have latitude and longitude of a point and a distance in degree, how can I get to the other point? (Suppose the point to find is aligned horizontally or vertically)


texts: Export the legend into its 100% size using Print Composer

texts:
I am trying to export the legend into an image using Print Composer into its 100% size as it is in the composer window. To be sure it is the intended size I want, I clicked Zoom to 100%:

But upon clicking Export as Image, I am forced to specify the image width and resolution, and there is no option for exporting into the 100% legend size as viewed in the Composer window:

Clicking Crop to Content doesn't help as I need to specify the resolution in dpi, which I need to painstakingly readjust just to get the 100% size. I also tried using 100 dpi as the export resolution but the resulting image is not its full size.

I had thought there should be a simpler solution than manually readjusting the dpi, but I don't know how. How can I export the legend into its 100% size?


texts: Is MATSim useful tool for Taxi fleet simulation?

texts:
I am working on a university project and I seek some advice:

My project consists in modeling the taxi fleet of a major city and investigating if replacing a percentage of the fleet with electric vehicles would make the taxi driver adapt their driving habits. Secondly, I would like to get valuable data of where to set up charging devices for taxis exclusively throughout the city.

My plan was to model it with MATSim, a multi agent traffic simulation software.

The research is based on a big amount of GPS data from 100 taxis over a time period of almost a year. So far, I have analyzed the data with Matlab and I will be able to export agent plans soon. In this case, I am modeling every taxi driver as an agent. I am not introducing a new mode „taxi“ because I can’t model the demand. But as I was working, I realized that maybe MATSim is not even suitable and now I am questioning my plan.

Do you think MATSim will be suitable?


texts:
I have a NetworkX graph corresponding to a mix of road and telecommunication network of a town, and different sets of nodes (of variable size) representing the location of network devices. I have to find an optimal path connecting all devices along the network and I tought the steiner_tree algorithmn available in NetworkX should do the trick.

Most of the time I get results within a reasonable time and with a reasonable use of resources (i.e. RAM). Sometimes however, when working with graph of few thousands nodes (around 5000) the process takes longer and eats a great amount of RAM.

I am searching for ways to reduce the complexity of the graph and I found the contracted_nodes function in networkx.algorithms.minors module. I use it to "contract" any node with a degree of 2, obtaining a single edge out of a sequence of consecutive edges.

The contracted_nodes has an optional self_loop parameter with the effect of preserving contracted edges as self loops.

Which can be the effects, if any, of these self loops on the steiner_tree algorithm?


texts: What are and how do I delete "spurious" artifacts?

texts: I dissolved the inner polygons of a shp file and was left with some lines inside. I'm not sure how they got there, since the original file does not seem to have them. I've also tried deleting them, but have not been able to do so.


texts: Shape of polygons does not adjust to new position after translating in QGIS

texts:
I want to translate the polygon in this shapefile (PMNM_py.shp, http://sanctuaries.noaa.gov/library/imast/pmnm_py.zip) from its initial position (NW of Hawaii) to central Europa to give somebody living in Europa an idea of the massive scale of this marine sanctuary using QGIS.

I can move it via editing mode/view->select features/edit->move feature(s) to its new position but I don't see the shape getting distorted (This should happen because I move the poly to another location in the same projection).

What am i doing wrong?


texts:
I am installing Postgis on my Debian server, following these instructions:

http://postgis.net/install/ (Ubuntu / Debian)

https://wiki.postgresql.org/wiki/Apt

However, after following the instructions it looks  like I have installed PostGres (9.4), but no PostGIS.

Do I have to install something more that was omitted in the instructions?

Can I somehow add PostGIS as a PostGres Plugin / Extension?


texts: Get size of a QgsLayoutItemPage in PyQGIS

texts:
I work on QGIS plugin in python, in my plugin I load several print layout in a project print layout manager. These print  layouts are either A4 or A3. I would like to get the size value (A4 or A3) of each print layout of my project print layout manager. I tried this code but it only returns the size in millimeters of the page.
project = QgsProject.instance()
composition = QgsPrintLayout(project)
document = QDomDocument()
template_file = open('C:/Users/Thysma/Documents/Stage/Composeurs/A4 Portrait.qpt')
template_content = template_file.read()
template_file.close()
document.setContent(template_content)
composition.loadFromTemplate(document, QgsReadWriteContext())
project.layoutManager().addLayout(composition)

layout = project.layoutManager().layoutByName('A4 Portrait')
page_size = layout.pageSize()
print (page_size)


texts:
I'm looking to spatial join points to their k nearest neighbors, with attributes summarized (mean, sum, sd, etc.).

I have looked at Join by Location (Summary) and Nearest Neighbor Analysis but the former doesn't allow nearness and the latter doesn't include a feature's attributes.

I imagine this is a duplicate but I've been unable to find anything.


texts: I have a polygon feature class in which the polygons frequently overlap each other. I want to find all of the unique combinations of overlapping polygons within that feature class. I believe I want to perform a Union on the feature class, but I also want to be able to concatenate the names so that each polygon has a unique name made up of the concatenated name of all the polygons that make it up. Any ideas on how to do this?


texts: Styling circle with radius depending on GeoJSON property in OpenLayers 3?

texts:
I want to represent the weight of each point of my geojson collection by a circle with a scalable radius.

Here an sample of my GeoJSON file:

"features": [
{ "type": "Feature", "properties": { "id": 20401, "radius": 0.032878 }, "geometry": { "type": "Point", "coordinates": [ 6.569957243078316, 46.744609055776721 ] } },
{ "type": "Feature", "properties": { "id": 20402, "radius": 0.024677 }, "geometry": { "type": "Point", "coordinates": [ 6.56996417375066, 46.744609055776721 ] } },

Is there a way to extract the radius property?

actually I represent my points like this

var points = new ol.layer.Vector({
         source: new ol.source.Vector({
         projection : 'EPSG:4326',
         url: 'weeds.geojson',
         format: new ol.format.GeoJSON(),
      })
  });

It's the default representation who appear on the map.


texts: Strange behavior for QgsProject snapping options - Python QGIS 2.8

texts:
I have some strange behavior regarding the snapping options function    setSnapSettingsForLayer in the QgsProject class.

I am setting the snapping parameters in the following way:

QgsProject.instance().setSnapSettingsForLayer(myLayer.id(), True, QgsSnapper.SnapToVertexAndSegment, QgsTolerance.Pixels, 10, False)

The snapping does not seem activated (there is no pink cross on my line layer that I called). When I open the snapping options via the GUI, what I first see is the snapping mode combo box on "Current Layer" with the snapping turned off. When switching to Advanced snapping, the layer I called in the method does indeed have snapping turned on with the correct parameters. Now, when clicking okay, the snapping is activated.

Any input on how I can get the snapping options turned on without having to navigate through the GUI?

Some more information, in case it has to do with the way the layer has been referenced--the layer is defined in this way:

layers = iface.mapCanvas().layers()
for layer in layers:
    if layer.name() == 'MyLayer':
        myLayer = layer

I am now testing it this way...but it doesn't seem to make a difference.

layers = QgsMapLayerRegistry.instance().mapLayers().values()
for layer in layers:
    if layer.name() == 'MyLayer':
        myLayer = layer

What I think is the problem:

I believe it has to do with the way the Current Layer snapping mode is set up. My code works okay when the project is saved with this mode on "sur un sommet et un segment" (endpoints and segments, in english, I believe). When the project has been saved with this mode on "Off" however, the advanced settings (set programmatically) do not override these settings. Is this correct? If yes, is there any way around it? Such as programmatically setting the "Current Layer" mode as well to ensure that it will be on the correct setting.

I found this QgsSnappingUtils::SnapToMapMode in the API, and am trying it this way to see if it gets me where I need (forcing the project to use snapping options based on the advanced mode, and not on current layer):

QgsSnappingUtils.setSnapToMapMode(QgsSnappingUtils.SnapAdvanced)

I am getting a TypeError though:

TypeError:QgsSnappingUtils.setSnapToMapMode(QgsSnappingUtils.SnapToMapMode): first argument of unbound method must have type 'QgsSnappingUtils'


texts: Does ArcGIS.com support dynamic map services (or only feature services and tiled services)?

texts:
Using ArcGIS Server it is possible to publish a map service, which has a URL like http://services.arcgisonline.com/ArcGIS/rest/services/Demographics/USA_1990-2000_Population_Change/MapServer/1

In this case, the "1" refers to the particular polygon layer. It's possible to create a featureLayer from this REST endpoint.

When I publish a map to ArcGIS.com the options are "Feature Access" or "Tiled Mapping". This means that the REST endpoint has /FeatureServer/ rather than /MapServer/ in the URL.

Am I correct in the assumption that ArcGIS.com doesn't support dynamic map services, and only supports feature services and tiled services?


texts: Color polygons using SQL Database based Rules in Openlayers

texts:
I want the rules to be applied according to values in a database. Is it possible to include external database information that is not in the shapefile?

In spoken words the original rule is "Color all vectors where the property:order07 is less than value:04000 with the color #ffeecc"

and I want to change it into something like "Color the vectors where the property:fromDatabase xy: is value:xxyy with the color #xxxxxx"


texts: Python Script to iterate through folder and copy all feature classes in MXDs to new file geodatabases?

texts:
I have a folder and would like to iterate through each MXD and take each SDE feature class that is in it and place in a new folder and file geodatabase, with a new MXD using the file geodatabase paths.

folder with original MXD name
new MXD in new folder with original MXD name
new file geodatabase with all feature classes from MXD

Each new MXD should have the above associated with it.

I have a solid understanding of each function that I would like to write aside from iterating and copying MXDs and their features, how should this be done with arcpy?

This script iterates through a folder and all MXDs and prints the MXD and feature classes used in the MXD.  I would like to further this by accessing each feature class from arcpy.mapping and sending those to a file geodatabase. And making a new MXD with the new feature classes.

import arcpy, os, datetime

folderPath = 'C:\MXD_test'

#Loop through each MXD file
for root, dirs, files in os.walk(folderPath):
    for file in files: # files is a list of files in the current directory
        if file.lower().endswith(".mxd"):
            fullpath = os.path.join(root, file) # root is the current directory
            #Reference MXD
            mxd = arcpy.mapping.MapDocument(fullpath)

            DFList = arcpy.mapping.ListDataFrames(mxd)
            for df in DFList:
                # Format output values
                if df.description == "":
                    descValue = "None"
                else:
                    descValue = df.description
                # ==== Note the new descName variable
                descName = df.name

                lyrList = arcpy.mapping.ListLayers(mxd, "", df)
                for lyr in lyrList:
                    lyrName = lyr.name
                    if lyr.supports("dataSource"):
                        lyrDatasource = lyr.dataSource
                    else:
                        lyrDatasource = "N/A"

                    print lyrDatasource
                    print fullpath


texts:
When I'm dissolving a shapefile:

ogr2ogr "b:\out.shp" "b:\Tele.shp" -dialect sqlite -sql "SELECT ST_Union(geometry), nev, jar FROM Tele GROUP BY jar"

I got error message like:

  Warning 1: One or several characters couldn't be converted correctly from UTF-8 to ISO-8859-1.

I don't understand this because Tele.cpg contains "UTF-8"

When trying this:

ogr2ogr "b:\out.shp" "b:\Tele.shp" -dialect sqlite -sql "SELECT ST_Union(geometry), nev, jar FROM Tele GROUP BY jar" ENCODING=UTF-8

the error message is:

  Warning 1: layer names ignored in combination with -sql.
  ERROR 1: Failed to create field name 'nev' : cannot convert to "UTF-8"

So what should I do to keep my strange Hungarian characters?


texts: Calculating Koppen climate classification in Google Earth Engine seems wrong

texts:
For a paper I am writing I need to know the extend of all drylands on earth. Instead of getting third party data I decided to calculate this myself so I keep in control of all steps of the process.

Online the definition of drylands according to Koppen can be easily found:

  Multiply the average annual temperature in Celsius by 20, then add

  (a) 280 if 70% or more of the total precipitation is in the spring and
  summer months (April–September in the Northern Hemisphere, or
  October–March in the Southern), or

  (b) 140 if 30%–70% of the total
  precipitation is received during the spring and summer, or

  (c) 0 if
  less than 30% of the total precipitation is received during the spring
  and summer.

  If the annual precipitation is less than 50% of this
  threshold, the classification is BW (arid: desert climate); if it is
  in the range of 50%–100% of the threshold, the classification is BS
  (semi-arid: steppe climate)

I tried implementing these factors and calculating the arid and semi-arid areas but I am pretty sure that France, the Amazon rain forest, Hungary and even Thailand are no drylands

I have also taken the median of the extend of the TRMM dataset (2000 - 2017) to eliminate inter annual variability but this did not seem to make any difference.

The code I am using is:

function ClassifyYear(precipitationdata, temperaturedata, year){

  var data = precipitationdata.filterDate(year + "-01-01", (year+1)+"-01-01");
  var Temp = temperaturedata.filterDate(year + "-01-01", (year+1)+"-01-01");
  var AvTemp = Temp.mean().subtract(273.15); //Average Annual temperature

  var PPTns = data.filterDate(year+"-03-20", year+"-09-22").sum(); //Precipitation in the summer of the Northern Hemisphere
  var TotPPT = data.sum();
  var PPTss = TotPPT.subtract(PPTns); //Precipitation in the summer of the Southern Hemisphere

  var mask30 = PPTns.lte(TotPPT.multiply(0.3)); //Masking pixels where summer PPT is less then 30% of the total PPT
  var mask70 = PPTns.lte(TotPPT.multiply(0.7)).subtract(mask30); //Masking pixels where summer PPT is less then 70% of the total PPT
  var mask100 = mask70.remap([0,1], [1,0]); //Masking pixels where summer PPT is more then 70% of total PPT

  var Threshold = AvTemp.multiply(20).add(mask70.multiply(140)).add(mask100.multiply(280)); //Calculating Threshold PPT
  var KoppenDry = TotPPT.lte(Threshold); //If PPT is less then Threshold -> Dryland
  var KoppenArid = TotPPT.lte(Threshold.multiply(0.5)).rename("arid"); //if PPT less then 0.5*Threshold -> Arid
  var KoppenSemi = KoppenDry.subtract(KoppenArid).rename("semi-arid"); //Dryland - Arid = Semi-Arid
  var Koppen = KoppenArid.addBands(KoppenSemi);

  return Koppen;
}

var list = [];
for (var year = 2000; year <= 2017;year++){ //Taking the median climate classification

  var result = ClassifyYear(TRMM.select("precipitation"),GLDAS.select("Tair_f_inst"), year);
  list.push(result);
}
print(list);

Map.addLayer(ee.ImageCollection(list).median());

Does anyone see a mistake that I am missing or knows more about the Koppen classification and finds what wrong with the calculation? The results I am getting:


texts:
I'm trying some code to get the building footprints of Barcelona:

import osmnx as ox
import matplotlib.pyplot as plt
%matplotlib inline

place_name = "Barcelona, Barcelona, Catalonia, Spain"
graph = ox.graph_from_place(place_name)
buildings = ox.footprints_from_place(place_name)

When I run the code it yields the following error message:

TopologicalError: The operation 'GEOSDifference_r' could not be performed. Likely cause is invalidity of the geometry <shapely.geometry.polygon.Polygon object at 0x00000262079B6550>

I've tried updating the OSMnx version but it seems not to be the problem.
How can I solve this issue?


texts: Calculate percentage of each biome on each national state in QGIS

texts:
I have one shapefile of Brazilian states (black lines) and another shapefile of Brazilian biomes (colored areas).

I want to add, on the states' attribute table, the percentage of each state's area occupied by each biome. The table would become something like this:
  State.Name    Amazon  Cerrado  Caatinga  Atlantic.Forest ...
Amazonas           100        0         0                0
Tocantins           10       90         0                0
Minas Gerais         0       50        10               40
Rio de Janeiro       0        0         0              100
...

The names of the biomes would come from a field in the biomes shapefile.
I have found methods for single polygons (not for a table with dozens, or thousands, if done for municipalities -- which is my actual goal). Also found a method to calculate the area if the biomes were a raster, but here they are a vector, and I'd rather not rasterize it. Found a Python script too, but I don't dominate that language, and I think this should be simpler than that.
How can I do that in QGIS?


texts: Different results from cs2cs and gdaltransform

texts:
I expected to get the same result when performing a transformation using cs2cs and gdaltransform. However, the results differ:

echo "53.70015N 9.12875E" | cs2cs EPSG:4326 EPSG:3857
1016207.80      7113571.58 0.00

echo "53.70015N 9.12875E" | gdaltransform -rpc -s_srs EPSG:4326 -t_srs EPSG:3857
5977873.35352241 1020534.697294 0

cs2cs -v
Rel. 6.2.1, November 1st, 2019

Gdal version is 3.0.2 using the same Proj4 lib.

Maybe I have completely misunderstood the two programs. When I use the results from cs2cs as corners for gdal_translate, I get the expected results, however:
wget http://b.tile.openstreetmap.org/8/134/82.png
gdal_translate -of PNG -outsize 256 0 -projwin 939825 7200381 1095810 7045677 -projwin_srs EPSG:3857 /usr/local/src/gdal-3.0.2/frmts/wms/frmt_wms_openstreetmap_tms.xml hemmoor_82.png

Will produce the same output.


texts: One way replication - Parent Database read only access

texts:
I am planning to setup one way replication from parent to child . My team only has read access to the parent database (SQL server) and child database is our own managed sql server.

Is it possible to setup one way replica with read only access to parent? If not, what are my other options?


texts: Matching Landsat 8 with previous Landsat images

texts:
I am working with Landsat 8 in ENVI. My image's dimensions are 8921 x 13131. I need to apply a forest mask to this image with dimensions of 8805 x 12321. I try to resize an image, but it doesn't work. (After that I need to calculate forest disturbance index for 2013 and subtract it from DI for 2000, so I really need  an 2013 image (Landsat 8) to be resized).

I resize it using basic tool - spatial subset - subset using - file - and then pick my forest mask. However, as a result I get a file of 8806 x 12321 instead of 8805. If I adjust the number of samples manually, it gives me 8805. But then I try to apply a mask and that's what I get (see a screenshot).  They just don't match. And when I try to overlay my mask (which I made based on Landsat 7 image) with my Landsat 8 image, those two images have different frames.

Does anybody have any good advice?


texts: Drop Marker with 3857 projection in leafletjs

texts:
Is it possible in leaflet to drop marker with 3857 projection on map?

I found a solution where i first translate a point to 4326 using proj4js and than drop it on a map.

var source = new Proj4js.Proj('EPSG:3857');
    var dest = new Proj4js.Proj('EPSG:4326');
    var p = new Proj4js.Point(-12307.601047734, 6711475.1283642);
    Proj4js.transform(source, dest, p);

new L.Marker([p.y, p.x], {bounceOnAdd: true}).addTo(map);

Is there any other way in core leaflet or leaflet plugin to drop marker without using proj4js library?


texts:
I am designing a QGIS project for collecting data in the field using Qfield or INPUT as mobile app. So the question is how to know if the surveyor is standing within the parcel when collecting data.
I thought that the solution is to calculate the distance between the surveyor's position and the centroid of the parcel, but I couldn't find the formulas to capture the surveyor's position


texts: PyQGIS point buffer with OpenLayers base map

texts:
Using PyQGIS, I desperately need help on creating a buffer around a point in a node layer and then selecting features from a trails layer that fall within the buffer. When I load the shape files without OpenLayers, it works. If I load a base map from OpenLayers like Bing Aerial or Google Satellite, the buffer no longer finds any features from the trails layer.

Also, zoom to location works great without OpenLayers, but when I bring in the map, it zooms in so far that all I see is a blank screen, and I have to manually zoom out quite a bit to see the shape files.

The base map is critical for the project or I would just do without.

Here is the code for selecting the trail features within the node buffer. If the feature we are looking for is not part of the selection, the code then zooms to the buffer bounding box around the node so the user can manually select a new trail. The base map is necessary in order to manually make a proper selection.

def check_trail(self, node_id):

    # Check if trail record has changed.

    found = 0
    canvas = self.iface.mapCanvas()

    request = QgsFeatureRequest().setFilterFid(node_id)
    feature = self.nodes_layer.getFeatures(request).next()

    pntGeom = feature.geometry()
    pntBuf = pntGeom.buffer(canvas.mapUnitsPerPixel()*30, 0)
    rect = pntBuf.boundingBox()

    # Convert the rectangle coordinates to the layer CRS.
    layerRect = canvas.mapRenderer().mapToLayerCoordinates(self.trailfile, rect)

    # save currently select trailfile
    trail_current = self.trailfile.selectedFeaturesIds()

    # remove current trailfile selection
    self.trailfile.removeSelection()

    # Select the features that are in the rectangle and return their fids.
    self.trailfile.select(layerRect, False)

    # get selected trails
    trails_found = self.trailfile.selectedFeaturesIds()

    for trail in trails_found:
        self.trailfile.removeSelection()
        self.trailfile.select(trail)
        temp_trail_name = self.trailfile.selectedFeatures()[0].attribute('name')

        # compare buffered trail name with the most recent trail name obtained
        if temp_trail_name == self.trails[-1]["name"]:
            # add trail
            found = -1
            break

    if not found:
        # Zoom to location where there is a new trail
        self.trailfile.removeSelection()
        selectList=[node_id]
        self.nodes_layer.setSelectedFeatures(selectList)
        box = self.nodes_layer.boundingBoxOfSelected()
        self.iface.mapCanvas().setExtent(box)
        self.iface.mapCanvas().refresh()
        self.trailfile.select(trail_current)

        # select trail to attach to segment
        while not self.selectFeatureResult and not self.cancelProdUnits:
            # Set the canvas click response to select a trail.
            self.clickResponse = select_trail
            self.selectTrail.show()
            self.selectTrail.exec_()
        self.selectFeatureResult = 0


texts: Geospatially match cells of different rasters - possible?

texts:
I have 2 rasters that were originally vector polygons. Raster2 was created from a selection of raster1. After converting both to raster, their cells don't match up (i.e. attribute info for each cell of raster1 doesn't match up with the coinciding cell for raster2). I want to be able to select a cell in raster1 with the value of 3 and know raster2's cell with the value of 3 is in the same location. Is there a way to match them?

Note: both have the same cell size.


texts: Open Layers 5 - Dynamically toggle zoom buttons

texts:
I would like to disable the zoom in button if a user tries to zoom in past a certain level. I would like to disable the zoom out button if a user tries to zoom out past a certain level.

I see that there are minZoom, zoom and maxZoom properties in a View object. Those do correctly prevent a user from zooming past where they should go. However the zoom in and zoom out buttons still look like they work because they are not disabled.

Does anyone have a code snippet, or could point to an example of how buttons on the ol map can be programmatically disabled?


texts: texture fill doesn't display in ArcScene

texts:
My texture fill doesn't display properly in ArcScene but it does in ArcMap. What am I doing wrong? Help!

ArcMap:

ArcScene:

The same texture fill has been selected in both ArcMap and ArcScene. Why does one show up and the other get blotted out?


texts:
I convert Bavarien NAS-Data (DHDN 3 Degree Gauszone 4) via ogr2ogr into shapefiles:

"C:\Program Files (x86)\QGIS Wien\bin\ogr2ogr.exe" "ESRI Shapefile" -skipfailures Entenhausen.xml

This seems to work perfectly, but on loading Data from other sources I found out that the X-Coordinate value in the converted files is 6 digits long (548521,69) whereas the Y-Value is 7 digits long (5376490,59) - the leading digit of the x value (wich specifies the Gauss-Krüger zone 4) present in data from other sources is missing.

Is there any (easy to follow) possibility either in ogr2ogr or qgis to reproject the x-values? I think ogr2org would be better as there are quite a lot of shapefiles.


texts:
I have created a java function to plot the points on the map, and I need to pass the x-y coordinates and the taskevent as the parameters. I would like to ask what should be the input parameter for taskevent??

public void plotCase(TaskEvent event){

  event.getWebContext().getWebMap();

}

Thanks a lot!!