texts:
I am creating a Plugin, where on addFeature button, first user will get a dialog box to select the type of feature user wants to draw. On the basis of users choice, other editing constraints will get activated as defined for a layer. Currently I am able to create a dialog box with radio buttons for users selection and things are working correctly for first run, if user is following below steps in order.

Start editing (Toggle Editing)
Add Feature button
Dialog box for users choice
Editing by the user as per applied constraints
Stop Editing (Toggle Editing)

The above steps are working fine for first run, but if user starts another editing session, then the dialog box is getting for 2 times and so on for further sessions.
While working on this, I found layer.editingStarted.connect(function), this signal is getting continued even after complete execution of function. I tried to disconnect this but not able to do it. Kindly suggest how I can disconnect the same or any other workaround for this problem.
My Code and Dialog box are attached for reference. Please check and suggest.
lyr=self.iface.activeLayer()
lyr.editingStarted.connect(lambda:fiber_add_feature(db,lyr))
def fiber_add_feature(self,db,vlyr):
    if vlyr.name()=="Fiber Cable":
        iface.actionAddFeature().triggered.connect(lambda: self.fiber_route_type(db,vlyr))
def fiber_route_type(self,db,lyr):
    if lyr.name()=="Fiber Cable":
        self.route_dlg=FiberRouteTypeDialog()
        self.duct_require_dlg=DuctRequirementDialog()
        if lyr.isEditable():
            ar="Aerial"
            ugr="Under Ground Route"
            Duct_presnt=True
            self.route_dlg.show()
            route_result = self.route_dlg.exec_()
            if route_result:
                if self.route_dlg.ar_rb.isChecked():
                    fiber = fiberEditTool(self.canvas, self.iface,db,ar, not(Duct_presnt))
                    self.canvas.setMapTool(fiber)
                elif self.route_dlg.ugr_rb.isChecked():
                    self.duct_require_dlg.show()
                    duct_result = self.duct_require_dlg.exec_()
                    if duct_result:
                        if self.duct_require_dlg.cb_duct.isChecked():
                            fiber = fiberEditTool(self.canvas, self.iface,db,ugr,Duct_presnt)
                            self.canvas.setMapTool(fiber)
                        elif self.duct_require_dlg.cb_wo_duct.isChecked():
                            fiber = fiberEditTool(self.canvas, self.iface,db,ugr,not(Duct_presnt))
                            self.canvas.setMapTool(fiber)


texts: How to add a raster to ArcMap if file name is unknown?

texts:
I have a raster file for which I know the the workspace, but not the filename (it is a temporary file generated by operations from my program that randomizes its name on creation, i.e. "ras457"). I cannot store the raster at time of generation for various reasons, so is there a way to programatically grab the raster from the folder without knowing its exact name and store it in an IGeoDataset?

The raster is stored in a temporary folder created when I run the program, so it will be the only item in the folder.

One route I have considered exploring is using .NET's directory info class to get the names of the contents of the folder and extracting it from there, but there may be a simpler way to accomplish this that I have not considered.


texts: Showing multiple attributes in legend of ArcMap?

texts:
I have two questions regarding creating the legend in ArcMap 10.2.2:

In the attribute table of my shapefile I change the number of decimals from 6 to 2; I can see the changes in the attribute table but in the legend and symbology still the full values are shown.
So the question is how can I have the number of decimals that I need in symbology and  legend.
I created a legend to show two fields by using " categories" option of symbology.
As you see the values are not in a column (Image).
I want to know How can I have two  nice individual columns for my text and values.


texts: Clip and Zip Coordinate system options Geocortex

texts:
I have the standard geoprocessing service, Clip & Zip, which references a spatial reference folder containing various coordinate systems (in prj).
The service works in geocortex via the standard 'Extract.xaml' file. However I need the ability to Change coordinate systems on export. The reason being, ESRI Base maps are on WGS so the geocortex instances / map services are WGS.

Anyone had any luck getting clip n zip working with the ability to select a coordinate system via geocortex? We are happy to have the data exporting in its proprietary coordinate system but would prefer the ability to choose a coordinate system, just not the web viewer coordinate system.


texts: Why does PostGIS ST_Distance return 0.0?


texts: Using Himawari 8 in Leaflet.js from NASA GIBS

texts:
I'm currently developing a Leaflet based map webpage and want to display satellite imagery, specifically from the Himawari 8 satellite. NASA have the GIBS api which appears perfect for my purposes. Using this Leaflet plugin I can get it to work for the MODIS (and all other polar orbitals) fine as here:
https://jsfiddle.net/80c71t2d/1/
However, this plugin predates the geostationary satellite feeds. If I try to use the geostationary satellites (i.e. Himawari 8) it seems to almost work at times but either no images are displayed or I get weird fuzzy/blocky images:
https://jsfiddle.net/2fxdyo1j/1/
Note the only difference between these two is in the bottom of the JavaScript (sorry I'm new to making a jsfiddle so I wasn't sure how to keep it tidier; this is a stripped down version of my page).
MODIS
    var GIBS_Polar_Modis = new L.GIBSLayer('MODIS_Terra_CorrectedReflectance_TrueColor', {
    date: new Date(),
    transparent: true
    }).addTo(map);

Himawari
    var H8_B13 = new L.GIBSLayer('Himawari_AHI_Band13_Clean_Infrared', {
    date: new Date(),
    transparent: true
    }).addTo(map);

This is the specification of the layers:
    L.GIBS_LAYERS = {
      "Himawari_AHI_Band13_Clean_Infrared": {
        "title": "Himawari_AHI_Band13_Clean_Infrared",
        "template": "http://gibs.earthdata.nasa.gov/wmts/epsg3857/best/Himawari_AHI_Band13_Clean_Infrared/default/{Time}/{TileMatrixSet}/{TileMatrix}/{TileRow}/{TileCol}.png",
        "zoom": 6,
        "date": true},
      "MODIS_Terra_CorrectedReflectance_TrueColor": {
        "title": "MODIS_Terra_CorrectedReflectance_TrueColor",
        "template": "http://gibs.earthdata.nasa.gov/wmts/epsg3857/best/MODIS_Terra_CorrectedReflectance_TrueColor/default/{Time}/{TileMatrixSet}/{TileMatrix}/{TileRow}/{TileCol}.jpeg",
        "zoom": 9,
        "date": true},
    }

Given the plugin predates the geostationary satellite feed, I had to work out what to put in. I can see in the browser console that files are getting obtained from the NASA API as I expect but certainly nothing that looks correct to me.


texts: Using a feature iterator Value variable in a SQL expression of select by attribute

texts: I am using Iterate feature selection to select dealers one by one as input features to the Near tool. I'm using Select by Attributes to select sales points from that same dealer. How do I sync the Select by Attribute (dlr_no=%value%) with the current feature selection that is being selected by the iterator?


texts: How to configure modal widgets in ArcGIS Server JavaScript API version 4?

texts:
The ArcGIS Server JavaScript API version 4 has the concept of widgets, such as the LayerList or Legend widgets.

Is it possible to set these widgets to be modal, whereby opening any widget closes all other widgets in the same mapView?

I can't see a way to obtain a list of all widgets on a mapView, or even to close a widget once you have a handle on it.


texts:
I am pretty new in working with WFS services and need some help to understand the basics.
This is the URL of the WFS I want to work with:
https://haleconnect.com/ows/services/org.732.341f2791-919e-49de-8d86-3b18e040c430_wfs?SERVICE=WFS&REQUEST=GetCapabilities&VERSION=2.0.0
When I connect it there is no problem, but when I want to add any of the listed layers I get an error message saying: Layer is not valid: The layer pagingEnabled='true' preferCoordinatesForWfsT11='false' srsname='EPSG:4326' typename='tn-a:AerodromeNode' url='https://haleconnect.com/ows/services/org.732.341f2791-919e-49de-8d86-3b18e040c430_wfs' version='auto' is not a valid layer and can not be added to the map. Reason:
What could be the reason for that?


texts: How to get a list of supported gdal formats within python

texts:
Is there a way to get a list of installed gdal formats on a system with python? On the command line this can be done easily with:

gdalinfo --formats

...but can it be done in python?


texts: Why is the "Contour" tool returning contours with an odd rectangular pattern?

texts:
I am generating 2' contours made from a smoothed DEM that is 3.2808 foot (1m) horizontal resolution. The result shows these odd squares that are about 80' on a side. The smoothing was done using 8 neighbors (EDIT: I thought I used focal majority but I checked my notes and I actually did use a Focal Mean. I was also wrong about the 8 neighbors. It was a 9x9 neighborhood.). The DEM doesn't show these squares (see second screenshot below).

I drew one of the squares with the drawing tool to demonstrate, but I think they are pretty obvious.
Q) Can anyone provide any explanation for this, or a potential solution?

This result is from the "Contour" tool in "Spatial Analyst Tools -> Surface -> Contour" and I'm using 10.3.1

Second screenshot showing the smoothed DEM (same general area but zoomed in a little)


texts:
Finally have gotten profecient at maplex labelingg and tweaking each option however, I am running into the issue of some labels appearing too close to each other and they are always duplicates. However, I want to keep duplicating street names every 3 inches. Any ideas?

Referring to the Fort Meade Road in yellow 2 times.


texts: Having CARTO legend items act as Radio Button?

texts:
Currently, I have a web map with multiple layers.

I would like the final user to be able to select one layer can be turned on at a time.

Currently, multiple layers can be turned on, but only the layer on top is visible.


texts: Using clip tool (ArcGIS) for iteration/looping 2 GDBs in : one for clipped data, second - for clipping features

texts:
I have 2 GDBs. The first one (let's name it X.gdb) contains feature classes (A,B,C...n) that are to be clipped. The second gdb (Y.gdb) contains fc (let's name it border1, border2, border3...) with only one polygonal feature (municipal district) in each. These are clipping features. What I need to do is to create model/python script, that creates separate gdb for every clipping feature (border1.gdb, border2.gdb, border3.gdb...) with clipped features from X.gdb

What kind of mistake did I make?


texts: QGIS saving analysis output as Postgis table

texts:
i am using QGIS and Postgis together. I can load data in Postgis and than i import it to QGIS. But i am trying to save an analysis output (example: buffer analysis, voronoi polygon or SEXTANTE analysis tools) directly in Postgis. Is there anyway to make this?

For instance in ArcGIS you can save an analysis output in a .GDB .MDB geodatabase but in QGIS  firstly i save an output as .SHP than i load it to Postgis or in OpenJump after analysis you must save layer in Postgis table. So there is one more step !


texts:
geolocation: function () {
    var self = this;
    self.geo = new ol.Geolocation({
        trackingOptions: {
            enableHighAccuracy: true
        },
        projection: self.Map.olMap.getView().getProjection()
    });
},

I have this code but if i want to get current position, i get "undefined" message.Why?

geo.getPosition() //return undefined


texts: Creating one DEM from multiple ASC files

texts:
I am struggling with this several days and now its time to ask how to do it in QGIS.

I have very detailed 27 ASC grid files with overall size of 9 GB created from 1m contours of large area (600 m2).
The goal is to merge all this data into one DEM for easier handling (output will be smaller i think ?)
and later to create a slope, hillshade color shaded reliefs ,etc for whole area easier.

Available Data:

I have simplified contours for whole area in SQL Sever which are much less in size than original ones.
ASC files are created from original contour data and I would like to create one DEM from that.

I have tried next QGIS/GRASS solutions with simplified contours of a whole area:

How to generate a DEM from a contour Shapefile?
Creating DEM from contours in QGIS?
http://grasswiki.osgeo.org/wiki/Contour_lines_to_DEM

but I had only success in rasterizing contours (v.to.rast.attr). That creates a 218 MB TIF with dimensions of 37945 x 30143.
But, in the next step when I start 'r.surf.contour' it just hangs forever (last try was 24 hour long and then I quit, although process was very active and consumed almost all memory (6GB RAM). Is this too large dataset and this is normal execution time on modern PC configuration ?

Other way would be to use ASC files and to merge it to one DEM. If someone could give me information how to do it efficently or some detailed workflow would be great.

I would say that the purpose of DEM will be to create a shaded relief as a background for printing 1:5000 maps.

UPDATE:

Well I have converted all 27 ASC to GeoTiff and also realized there are 3 more grids (.grd) which has been made by GS Surfer 12. I converted them also in GeoTiff.
Now when I try to build VRT next thing happens:

The first of 27 GeoTifs are joined into VRT but I've lost all elevations. VRT have only one value (167, which is the lowest value from all dem's included).
Those 3 files built from Surfer GRD have somehow different characteristics and it wont join with others.

**

> Warning 6: gdalbuildvrt does not support heterogenous band
> characteristics. Skipping D:/DEM (GeoTIFF)/OF_IZO_28_dem.tif

**

gdalinfo for one tif in original group:

Driver: GTiff/GeoTIFF
Files: D:/DEM (GeoTIFF)/OF_IZO_27_dem.tif
       D:/DEM (GeoTIFF)/OF_IZO_27_dem.tif.aux.xml
Size is 5001, 5001
Coordinate System is `'
Origin = (7584999.500000000000000,4795000.500000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Image Structure Metadata:
  COMPRESSION=PACKBITS
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 7584999.500, 4795000.500)
Lower Left  ( 7584999.500, 4789999.500)
Upper Right ( 7590000.500, 4795000.500)
Lower Right ( 7590000.500, 4789999.500)
Center      ( 7587500.000, 4792500.000)
Band 1 Block=5001x1 Type=Float32, ColorInterp=Gray
  Min=243.000 Max=834.072
  Minimum=243.000, Maximum=834.072, Mean=503.092, StdDev=130.621
  NoData Value=-1
  Metadata:
    STATISTICS_MAXIMUM=834.07208251953
    STATISTICS_MEAN=503.0916431759
    STATISTICS_MINIMUM=243
    STATISTICS_STDDEV=130.62100147903

gdalinfo for one tif made from Surfer .grd>

Driver: GTiff/GeoTIFF
Files: D:\DEM (GeoTIFF)\OF_IZO_31_dem.tif
       D:\DEM (GeoTIFF)\OF_IZO_31_dem.tif.aux.xml
Size is 2125, 3512
Coordinate System is `'
Origin = (7587875.569015772100000,4790000.499983053700000)
Pixel Size = (0.999968455743971,-0.999966106522373)
Image Structure Metadata:
  COMPRESSION=PACKBITS
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 7587875.569, 4790000.500)
Lower Left  ( 7587875.569, 4786488.619)
Upper Right ( 7590000.502, 4790000.500)
Lower Right ( 7590000.502, 4786488.619)
Center      ( 7588938.036, 4788244.560)
Band 1 Block=2125x1 Type=Float64, ColorInterp=Gray
  Min=409.000 Max=1239.595
  Minimum=409.000, Maximum=1239.595, Mean=679.431, StdDev=154.121
  NoData Value=1.7014100000000001e+038
  Metadata:
    STATISTICS_MAXIMUM=1239.5948358501
    STATISTICS_MEAN=679.43099985765
    STATISTICS_MINIMUM=408.9999998738
    STATISTICS_STDDEV=154.12086745843

How I would correct this and get all grids in VRT with preserved elevations ?


texts:
Any idea how to add a tile layer (XYZ-Layer) using PyQGIS 101 to a map?

Something like osmLayer = QgsXYZLayer(uri.uri(), "OpenStreetMap", "XYZ")

"XYZ" is a raw guess.

Running QGIS 3.6.0


texts: Using Single Polygon for entire dataframe in PySAL?

texts:
I'm trying to calculate the queen contiguity weights using pysal's queen_from_shapefile() helper method based on this page

Not being a GIS specialist I wanted to check my intuition on this Q&A site first before proceeding further in my project. I'm using Shapely, GeoPandas and PySAL. Code follows, then question:

import pandas as pd
import geopandas as gpd
from shapely.geometry import Point, Polygon
import requests
import json
import time
import sys
import pysal as ps

testdf = (pd.DataFrame({'postCode': ["NE18 0AD", "TF4 3AL", "PL3 4PJ", "DD8 4BZ", "ML8 4FJ",
                                   "TW20 2DQ", "CA26 3AS", "M3 4DT", "PA72 6JH", "EC3N 1LA"],
                     'hasAttribute': [0, 1, 0, 0, 1, 0, 1, 0, 0, 0]}))
//random postcodes from here: https://www.doogal.co.uk/PostcodeGenerator.php
postcodelist = list(testdf.postCode)
url = 'http://api.postcodes.io/postcodes'
headers = {'content-type': 'application/json'}
lats = []
longs = []
payload = {'postcodes' : postcodelist }
r = requests.post(url, data=json.dumps(payload),  headers=headers)

time.sleep(0.1)
if r.status_code is 200:
    returnData = r.json()
    sys.stdout.write('.')
else:
    print("Query failed, status: %d".format(r.status_code))

result = returnData['result']
# Lets comprehend the list to get the lat and log data
# We need an if statement to catch cases where the postcode cant be found
lats += [ result[idx]['result']['latitude'] if result[idx]['result'] is not None else np.nan for idx in range(len(result)) ]
longs += [ result[idx]['result']['longitude'] if result[idx]['result'] is not None else np.nan for idx in range(len(result)) ]
testdf['lat'], testdf['long'] = lats, longs

polygon_geom = Polygon(zip(list(testdf.long), list(testdf.lat)))
polygon_list = [polygon_geom for i in range(testdf.shape[0])]
crs = {'init': 'epsg:4326'}
polygondf = gpd.GeoDataFrame(testdf, crs=crs, geometry=polygon_list)
polygondf.to_file("testpoly.shp")
qW = ps.queen_from_shapefile("testpoly.shp")

Is it correct to assume the same Polygon for all rows of the dataframe for calculating queen weights and, if not, what is the correct approach?


texts:
I am using ArcMap 10.4.

I’ve been trying to implement a python add-in that runs when Arcmap is opened.

When the script does run it joins a layer with a feature service layer. I have everything working if you run it in the Arcmap python window. When I take the next step and make it into an add-in it will not work. I added a sleep time to it to delay running until arc is fully open in hopes that it would fix the problem, but no such luck.

I performed a test with a text file and it will update the text file when arc is opened so I know that it is indeed running, it’s just not doing what I want it to. Any ideas?

Here is what I have:

import arcpy
import pythonaddins

class RoadJoin(object):
    """Implementation for RoadJoin_addin.RoadJoin (Extension)"""
    def __init__(self):
        # For performance considerations, please remove all unused methods in this class.
        self.enabled = True
    def startup(self):
        pass
        # Local variables:
Roads_2015 = "Roads 2015"
Roads_2015__2_ = Roads_2015
All_roads = "private\\transportation_FA_wa\\All roads"

time.sleep(30)
        # Process: Add Join
arcpy.AddJoin_management(Roads_2015, "BRUCE_ID", All_roads, "BRUCE_ID", "KEEP_ALL")


texts: Creating polygon from visible map extent using ModelBuilder?

texts:
Is there a way to create a polygon using the visible extent of a data frame in ModelBuilder?

I'm using ArcGIS Desktop 10.0.


texts: Which ArcGIS product gives full 3D terrain navigation and (flight) simulation?

texts:
I looked at the ArcGIS ESRI website and find that 100% of the screenshots depicts a 2D map with at most relief shading.

Is there any product ArcGIS or not that gives real 3D terrain visualisation capable of things like flight simulation?

Thanks.


texts: Measuring distance traveled, in 3D using ArcGIS Pro?

texts:
I'm working in ArcGIS pro, I have tracking data from a bird's flight path and want to know how far the bird flew. I started with Google Earth KML Layer and imported to ArcGIS Pro. This layer has data for each point along the bird's flight path at regular intervals. In the attribute table, each point has an X, Y location and altitude height.

We want to measure the total distance the bird flew. The bird flew up and down, in circles, loops, etc. So we are not just measuring in 2 dimensions, but the altitude is a factor as well.

We looked at the Generate Near Table tool but this does not account for altitude.
We can measure distances in 3D using the measure tool in ArcGlobe. However, we have over 6000 points.
Lastly, we looked at the near 3D tool, but this seems to measure between feature classes, instead of items within the same feature class.

Can anyone recommend another tool or workflow for us to use?


texts: qgis2web not taking my rule-based styles

texts:
It is the first time I post here. Is qgis2web working with rule-based styling from QGIS? I've created a few rules to have a dynamic style: polygons when the scale is under 12,500 and points beyond that scale. However, when exporting it to webmap it is not taking the styles. Am I missing something out?

Please see below how I have styled it. Unfortunately, I don't think creating a point layer is going to work for this project. Just to give you some context... my polygon layer is originally stored in a Oracle database, published in Geoserver and added to QGIS as WFS. We will be updating the polygons so I don't want to be updating the points every time something changes.


texts:
I just installed QGIS server using the instructions provided here. I am pretty sure that I have installed QGIS properly, but I am not sure. I am using Windows 7, and XAMPP. I typed the following into my browser:

http://localhost:8081/qgis/qgis_mapserv.fcgi.exe?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities&

And I get a page that begins like this:

<WMS_Capabilities xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns="http://www.opengis.net/wms" version="1.3.0" xsi:schemaLocation="http://www.opengis.net/wms http://qgis.org/wms_1_3_0.xsd">
<Service>
<Name>WMS</Name>
<!--  Human-readable title for pick lists  -->
<Title>QGIS mapserver</Title>
<!--
 Narrative description providing additional information
-->
<Abstract>A WMS service with QGIS mapserver</Abstract>
<KeywordList>
<Keyword>QGIS mapserver</Keyword>
</KeywordList>
<!--
 Top-level web address of service or service provider.
 See also OnlineResource elements under <DCPType>.
-->
<OnlineResource xmlns:xlink="http://www.w3.org/1999/xlink" xlink:type="simple" xlink:href="http://www.sourcepole.ch/"/>
<!--  Contact information  -->

I am pretty sure this is what I am supposed to get, so I assumed that what I was doing was correct. After that, I placed a QGIS file in C:\OSGeo4W\apps\qgis\bin and tried to view this by putting:

http://www.localhost:8801/qgis/qgis_mapserv.fcgi.exe?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities&map=C:/OSGeo4W/apps/qgis/bin/LakeCountyBlight.qgs&

in my browser.

I wasn't sure if that was what I was supposed to do, so I went into QGIS > vector > add WMS/WMTS > new > URL > localhost:8801/qgis/qgis_mapserv.fcgi.exe?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities&map=C:/OSGeo4W/apps/qgis/bin/LakeCountyBlight.qgs

When I did that, I got the error message:

Failed to download capabilities:
Download of capabilities failed: Host www.localhost not found

I'm not sure why nothing is happening. One thing that might be a problem is that my computer has both IIS, and Apache on it. I set IIS and Apache to different ports by going into the httpd.conf files. However, I also had to change the httpd.conf file for QGIS server to something different for Xampp to get it to work. I was wondering if that was the problem, and if so, how would I fix it?


texts: Which one is best open source GIS software where we can use  Google/bing/yahoo maps as a basemap? i want to update my shp files using any of these basemap.


texts: Alternatives to ArcGIS' "RegionGroup" a raster

texts:
I'm using the ArcGIS Python 2.7 interpreter in my PyCharm as I need several algorithms of the ArcPy package. However, I am having trouble with the RegionGroup of ArcGIS 10.7, as it sometimes crashes my loop with the exit code -1073741819 (0xC0000005).
I have already asked a question on how to fix that (Process finished with exit code -1073741819 (0xC0000005) from RegionGroup (ArcGIS 10.7)), but I am assuming that it is rather difficult to solve (at least Google seems to think so).
I need to stay with the 2.7 Python in my PyCharm, but maybe there are some alternatives to region group a raster?


texts: How to create a two sided ticks labels grid in layout (ArcMap)?

texts:
I'm making a planing map. I want to add to the final layout a specific  kind of grid:
The grid needs to be with a 6 character label that is separated in to two labels. Each one of the labels needs to be from both side of the grid's outer ticks. Like this example:

I've tried to change almost all of grid's properties but still could not able to figure out how to do it.

My result so far isn't satisfying:

It's important to say that I do not have the ArcGIS Production Mapping extension.
Is there any way to create this kind of grid without this extension? If so, how?


texts: Trying to install QGIS - missing 1.11 GDAL

texts:
I just spent the last n hours working on this, and still haven't gotten anywhere, so I'll throw it to the group. I see that other people have asked about this here and here, and I'd be happy to take down this Q if that gets answered.

My specs:

Mac 10.9.4
Mavericks OS

I'm trying to install QGIS. I'm working with the KyngChaos package. It gives me an error on installation, saying, "Failure: The GDAL 1.11 framework is required." I've since tried to brew install gdal, which installed the 1.10 framework. QGIS installation still fails. I then did brew update, and tried again (which seemed promising). Still fails. I then tried to follow this guy's approach, and brew install qgis. This kept my Terminal busy for about 3-4 hours... and still failed. Though I swear I saw a green-fonted gdal library get installed in all that Terminal chatter.

FWIW, the brew install qgis gave me a slightly different error message:

[ 50%] Built target qgis_gui
make: *** [all] Error 2

No idea what this means. I'm trying to find gdal in my computer via the Terminal: mostly doing which gdal -a, though this doesn't return anything. So I'm not sure how to even check that GDAL 1.11 got installed.


texts:
In my GeoServer I want to cache tiles for a layer. But when I am trying to seed tiles for a layer I can see the status as stated the image below.

  Tiles completed is -1.

I am refreshing but it still seems -1.

I cannot understand what the problem is or it is caching tiles or not. Also I cannot find any log in Tomcat log file (catalina.out).


texts: Setting raster resolution in meters when working in a Geographic Coordinate System

texts:
I work in WGS 84 (Geographic Coordinate System) in ArcGIS 10.1. I need to create a raster specifying its cell size in meters, and output it in WGS84.

Is it a good solution (towards a less shape/area distortion of the final result) to proceed in this way?: Create the raster from the shape which is in a Geographic CRS (WGS84), having set before the CRS of the Data Frame to a projected one (WGS84 UTM) in order to be able to give cell size in meters, which be a small value in relation with the actual required resolution in order to get accurate result when the expected resampling happens; then, once the raster is in a projected CRS, reproject it to WGS84 (the resampling is expected that involves a low shape/área distortion)?

IMPORTANT for this implementation to be practical: In the last step (reproject the raster to WGS84), when specifying the Cell size of the resulting projected raster, the units of the default value that appears are decimal degrees. Are this value coincident with the Cell size in meters of the raster from which I am projecting to WGS 84?

Another related doubt to solve: Which would be the proper projected coordinate system through the raster would pass, better than WGS84 UTM? A equal area one? an equal distance one? Or would UTM (which is a equal shape one) be the best?


texts: How to output points using TerraScan based on selected shape and text of block boundary

texts:
For those of you who are familiar with TerraScan, does anyone know how to save/output LAS files based on selected shape and text of block boundaries (tile index)?

The TerraScan User Guide is pretty useless.  I've tried selecting shapes and text, and running a macro with FnScanOutput(), using variable #block.las, Inside fence only.  I've tried creating a project with only the desired files, and running the macro without "inside fence only" selected on the whole project.  I can't seem to get anything to work.


texts: Grass : how to create attribute row for a new feature?

texts:
In a python script, I've created a vector and associated attribute table :

grass.run_command('v.edit', map='graph', type='line', tool='create')
grass.run_command('v.db.addtable', map='graph', columns='cat iNTEGER, adj INTEGER')
grass.run_command('v.db.connect', flags='p', map='graph')

then I create features in this map with :

grass.read_command('v.edit', flags='n', tool='add', map='graph',type='line', input=fichtemp)

But I can't find a way to create a row in the attribute table for a feature added.
How can I get the cat|ID of the feature just created?


texts: Train Images Classifier - OTB Spatial reference of input and samples don't match

texts: I am running the Train Images Classifier on OTB 7.1.0. I used QGIS to upload a raster images and then created a vector polygon file for my classes. I ensured that both the raster and vector files are in the same projected coordinate system. I checked both in R and QGIS. When I run the Train Images Classifier tool I receive an error (image 1). I look at both files coordinate system below the error, but they are an exact match (image 2). Image 3 contains the beginning of the analysis. Could anyone shed some light on this? Let me know if you need any additional information.


texts:
I can see that it's possible to load in csv data to R (as per here: Load a CSV as Simple Features R), however, it appears to me that there is a limitation on the range of delimiters allowed given the direct use of GDAL (COMMA/SEMICOLON/TAB/SPACE, see "Layer Creation options" here).

So with sp() I could use sep as follows:

pow_pointX <- read.csv("./data/poi_2015_12_scot06340459.csv", sep="|")

Given that the pipe character "|" is default for Ordnance Survey datasets distributed in the UK, this seems like a serious limitation. Has anyone found a workaround that doesn't require the use of sp()?


texts: Searching planet_osm_point by longitude and latitude

texts:
I used osm2pgsql to import data into my PostgreSQL database and it put gave me the following 4 tables, planet_osm_point, planet_osm_line, planet_osm_polygon and planet_osm_roads.

I have looked at the schema for osm2pgsql and inside the planet_osm_point table there is a column 'way' of type point. I tried to query it as an array like it says in the documentation here (http://www.postgresql.org/docs/current/static/functions-geometry.html):

  It is possible to access the two component numbers of a point as
  though the point were an array with indexes 0 and 1. For example, if
  t.p is a point column then SELECT p[0] FROM t retrieves the X
  coordinate and UPDATE t SET p[1] = ... changes the Y coordinate. In
  the same way, a value of type box or lseg can be treated as an array
  of two point values.

Anyone know the correct query the planet_osm tables from osm2pgsql based on long/lat?
For example if I wanted to retrieve the information for Big Ben (London) I would query along the lines of:

    SELECT * FROM planet_osm_point
    WHERE ST_Y(ST_Transform(way, 4326)) = '-0.12'
      AND ST_X(ST_Transform(way, 4326)) = '51.5';

There is no error message, it comes out with the table names, but no data is output either. (Its an empty table). What am I doing wrong here?

Perhaps it is too precise, i.e. there are no coordinates that match it exactly? Perhaps I need to try being more general, say 51 < long < 52 and -0.5 < lat < 0. Although this would throw back multiple results so would not be quite as accurate.


texts:
The overriden function is
public override void DrawDynamicLayer(esriDynamicDrawPhase DynamicDrawPhase, IDisplay Display, IDynamicDisplay dynamicDisplay)

What's the difference between dynamicDisplay and Display?  They both have similar functions like "DrawPolygon" and "DrawText"

and then how about if I cast my map to an ActiveView and access ScreenDisplay, how are all of these display types different?  How should I know which one to use?


texts:
I have thousand of points generated along river channels. I am trying to connect these points based on distance. I tried "Join attributes by nearest Points" and "geometry by expression" in QGIS. The method was mentioned here by Babel. It works somehow. However, there are many gaps between the line segments and some lines are repeated drew.
How can I use QGIS to fill the gap and remove the repeated line segments?


texts:
I've definitely got at least a couple of grey hairs while trying to convert some dwg. files to  shapefiles. I've seen that this subject has been disscused a lot, but hope someone will have the patience to help me.

The procedure I'm following is this:

I add the dwg files to ArcMap.
TOC- right clik-Convert CAD Feature Layer
Data Management Tools-Projections- Define Projection.

The conversion happened but then I realized things are not that simple, looked over the files again and apparently I have some "lost points"....

My data is acctually in that middle point and it looks like this:

I definitely have to georeference the data before converting it to shapefile. But how do I do that? I tried to create some control points, but I'm probably not doing it the right way cause nothing happens. Plus I have no idea where those lost points should actually be.

Forgot to mention that the dwg files have no coordinate system.


texts: Writing custom GUI map application without installing QGIS

texts:
I would like to know if it is possible to develop standalone map applications with QGIS without installing QGIS software or embedding QGIS. Most instructions ask to either install QGIS or embed in to the final executable.
Is it possible to call QGIS functionalities by using shared libs? According to this (Standalone python app that uses qgis.core without installing qgis) answer there's one example application, https://github.com/DMS-Aus/Roam which requires py2exe to be installed.

Is there a promising way to achieve the objective without installing QGIS?


texts: Efficient way to resample, crop, and mask a raster image using R

texts:
Does anyone know an efficient way to resample an image? I am currently trying to crop a large image and clip two regions from a shapefile. From this clipped image, I would like to increase the pixel density by a factor of 10. I would like to then mask this clipped and resampled image. Is there an efficient way to do this?

I am using an image with a resolution of 250m, and when I plot the .kml the pixels are nowhere close to the edges of my shapefile, so I found resampling and increasing the pixel density, then masking this image is the best way so far. The problem is the polygons I am clipping from are extremely far away from each other, so this means that right now the area in between has to all be resampled, and it is taking an extremely long time.

img <- crop(Canada,spdf)
proj4string(img) <- CRS("+proj=lcc +ellps=GRS80 +lat_1=49 +lat_2=77 +lon_0=-95")

factor <- .1
incol <- ncol(img)
inrow <- nrow(img)
resample_img <- raster(ncol=(incol/factor), nrow=(inrow/factor))
extent(resample_img) <- extent(img)
final_img <- resample(img,resample_img,method="ngb")
proj4string(maybe)<- CRS("+proj=lcc +ellps=GRS80 +lat_1=49 +lat_2=77 +lon_0=-95")

final_img <- mask(final_img,spdf)


texts: Is there a way to get at C's OGRLayer.Intersection() or OGRLayer.Clip() in Python?

texts: I am working with OGR in Python and noticed that the C libraries have useful Intersection() and Clip() functions as part of the Layer class. Is there any way to get at these functions in Python? I know that these functions exist on the Geometry level but I'm looking specifically at the Layer. Thanks!


texts: Rendering a map with Python works inside QGIS, but not in standalone script

texts:
I have an existing QGIS project and I want to write a standalone Python script that outputs a png file of the map. I work on Windows 10. This is the code that I use (basically taken from the PyQGIS cookbook):

import os
from qgis.core import *
from PyQt5 import *
from PyQt5.QtSvg import *
from PyQt5.Qt import *

app = QgsApplication([], True)
app.setPrefixPath(r"C:/Program Files/QGIS 3.10/apps/qgis", True)
app.initQgis()

project = QgsProject.instance()
project.read("filename.qgz")
layer = project.mapLayersByName("myLayer")[0]

options = QgsMapSettings()
options.setLayers([layer])
options.setBackgroundColor(QColor(255, 255, 255))
options.setOutputSize(QSize(800, 600))
options.setExtent(layer.extent())
render = QgsMapRendererParallelJob(options)
image_location = os.path.join(os.getcwd(), "render.png")

def finished():
    img = render.renderedImage()
    img.save(image_location, "png")
    print("saved")

render.finished.connect(finished)
render.start()

Executing these commands line by line in the Python console inside QGIS works perfectly fine and it outputs a png file with the map in a fraction of a second. However running this from the Python interpreter, it does not output anything. When calling render.isActive(), it returns True.

Is this a bug? Or do I have to do things differently when I want to render a map in a standalone script?


texts: Using Google Maps, the mouse cursor and the dot on the map not inline

texts:
I'm trying to implement the example http://openlayers.org/dev/examples/drag-feature.html .

When I'm implementing this with a Google map, I have a strange bug: the mouse pointer and the dot on the map are not in the same place.

Now if I change the line that in bold: map.addLayers([gmap,vectors]); to map.addLayers([wms,vectors]); everything works, (in other words I don't use Google Maps).

Can someone tell me what is happening?

Here is my page and code:

<title>OpenLayers: Google Layer Example</title>

OpenLayers: Google Layer Example


texts: Counting member features in ArcGIS table?

texts: I'm working on Arcgis 10.3 and have a table as shown below.  I need to count each member of the Name field according (and not according) to Class Field.


texts: Effect on using different spheroids system

texts:
Recently I am investigating case of grounding vessels, one of the findings is the paper chart used for passage planning was using Bessel 1841 spheroids, whilst the GPS system was using WGS 1984 spheroids. The waypoints' coordinates planned in the paper chart then inserted into the GPS (using exactly same coordinate values).

The navigating officer continuously maintains the ship track following the waypoint setup in the GPS not the paper chart. However the ship track went close to a shallow reef and then grounded.

Why could two spheroid systems result in shifting of ship position?


texts: Error handling for comboboxes in a Python Addin using ArcGIS for Desktop 10.5

texts:
I'm creating an addin that contains two comboboxes and a button, the latter of which adds data to the current map based on user selections in both comboboxes. Both comboboxes have pre-defined lists extracted from a folder of lyr files. If the user hasn't selected from either comboxes and clicks the button, I would like it to recognise this and show a dialogue box. Currently, if the button is clicked on startup, the global parameters litho and gde can't be defined. If items are selected from the comboboxes but then the text is deleted, then the button pressed, the last selection is remembered and previous data added to the map. I'd like the selection to reset to "" if the text is deleted in the box.

class fmb_toolAddLayers(object):
def __init__(self):
    self.enabled = True
    self.checked = False

def onClick(self):
    try:
        mxd = arcpy.mapping.MapDocument("CURRENT")
        df = arcpy.mapping.ListDataFrames(mxd)[0]
        gde_lyr_name = gde + " - FMB Paleo Map.lyr"
        litho_lyr_name = litho + " - FMB Lithopies.lyr"
        gde_layer = arcpy.mapping.Layer(os.path.join(gde_folder, gde_lyr_name))
        litho_layer = arcpy.mapping.Layer(os.path.join(litho_folder, litho_lyr_name))
        try:
            if not gde == "" and litho == "":
                arcpy.mapping.AddLayer(df, gde_layer, "TOP")
            elif gde == "" and not litho == "":
                arcpy.mapping.AddLayer(df, litho_layer, "TOP")
            elif not gde == "" and not litho == "":
                arcpy.mapping.AddLayer(df, gde_layer, "TOP")
                arcpy.mapping.AddLayer(df, litho_layer, "TOP")
        except Exception, ex:
            pythonaddins.MessageBox("Please select a layer(s) \n\n" + str(ex), "Error", 0)
        pass
        del mxd, gde_layer, litho_layer
    except Exception, ex:
        pythonaddins.MessageBox("Error adding layers to map \n\n" + str(ex), "Error")
        pass

# ------------------------------------------------------------------------------
#  TOOLBAR - COMBO BOXES
# ------------------------------------------------------------------------------

class fmb_comboGDE(object):
    def __init__(self):
        for root, dirs, files in os.walk(gde_folder):
            gde_files = [f.replace(" - FMB Paleo Map.lyr", "") for f in files]
        self.items = gde_files
        self.editable = False
        self.enabled = True
        self.dropdownWidth = 'WWWWWWWWWWWW'
        self.width = 'WWWWWWWWWWWW'
    def onSelChange(self, gde_selection):
        global gde
        gde = gde_selection
    def onEditChange(self, text):
        pass
    def onFocus(self, focused):
        pass
    def onEnter(self):
        pass
    def refresh(self):
        pass

class fmb_comboLithopies(object):
    def __init__(self):
        for root, dirs, files in os.walk(litho_folder):
            litho_files = [f.replace(" - FMB Lithopies.lyr", "") for f in files]
        self.items = litho_files
        self.editable = False
        self.enabled = True
        self.dropdownWidth = 'WWWWWW'
        self.width = 'WWWWWW'
    def onSelChange(self, litho_selection):
        global litho
        litho = litho_selection
    def onEditChange(self, text):
        pass
    def onFocus(self, focused):
        pass
    def onEnter(self):
        pass
    def refresh(self):
        pass


texts: I just want somebody to help me with a link so that I can download global climatic data variables like temperature, precipitation, cloud cover etc


texts: Regular expression delimiter for text file in QGIS 3

texts:
I have a text file containing coordinates but also comment lines such as :

#file description
#description point 1
id1 X1 Y1 Z1
#description point 2
id2 X2 Y2 Z2

and I want to import this file in QGIS 3.4. I need to build a regular expression for the delimiter that will allow QGIS to read lines with coordinates and pass comments beginning with #.


texts: Synchronizing Changes From File Geodatabase to SDE Geodatabase

texts: I have an SDE geodatabase that has a one way replication to a file geodatabase, where I will make edits to the feature classes offline somewhere. I know that this scenario doesn't support one-way or two-way replication, so is there anyway I can push the updates I made in the file geodatabase back to the SDE geodatabase?


texts: Use Python to Find duplicate values in a feature Class and populate a field

texts:
So I found this python script on another post here and was trying to adjust it for my needs.  I'm a very novice python user so I'm struggling with how to modify the script.  I have a Feature Class stored in a feature Dataset that I want to search a field for duplicate values and populate a new field with Y for duplicate or N for none.  The below script is what I found that looks like it will work once I find a way to drill down into my file geodatabase.

        from arcpy import *

        inShapefile = pointsShapefile

        checkField = "xyCombine"
        updateField = "dplicate"

       #List of values found once
        occursOnce = []
       #list of values found twice
        occursTwice = []

     cursor = da.SearchCursor (inShapefile, [checkField])
          for row in cursor:
          #Check value is not null
       if row[0]:
    #If not already found to occur twice, proceed
    if not row[0] in occursTwice:
        #If hasn't occured once yet
        if not row[0] in occursOnce:
            #Add to occurs once list
            occursOnce.append (row[0])
        #If value has already been found once
        else:
            #Add to occurs twice list (duplicates)
            occursTwice.append (row[0])
    del cursor

     cursor = da.UpdateCursor (inShapefile, [checkField, updateField])
     for row in cursor:
      #Check value is not null
        if row[0]:
           #check if value in occursTwice list (i.e. is duplicate)
    if row[0] in occursTwice:
        row[1] = "Y"
    else:
        row[1] = "N"
    cursor.updateRow(row)
   del cursor


texts: OpenLayers 3 doesn't want to draw programmatically a polygon

texts:
I am using OpenLayers 3 and I am trying to draw a polygon using given coordinates, but it doesn't want to draw it.

var source = new ol.source.Vector();

            var ring = [
             [3139880.24789847, 5961935.332187176], [3179627.5026067616, 5972025.01992082],
             [3146606.706387566, 5927997.291628557], [3186353.9610958574, 5939615.719927904]];

            draw = new ol.interaction.Draw({
                source: source,
                type: 'Polygon',
                geometryFunction: ring,
            });

            draw.on('drawend', function (e) {
                var id = guid();
                e.feature.featureID = id;
                e.feature.setProperties({
                    'id': id,
                    'name': 'Polygon',
                    'description': 'Some values'
                })
                map.removeInteraction(draw);
            });
            map.addInteraction(draw);


texts: How to build PostGIS database schema for flexibility and normalization?

texts: I'm creating a database that will be used by researchers in population genetics with different goals and sampling strategies. The database needs to accommodate locations as points, lines, or polygons and I'm puzzling over whether those should be kept in a single "locations" table, or in separate tables for "point_locations", "line_locations", and "polygon_locations". The single-table design is appealing because it would make it simple to enforce the required one-to-one relationship between samples and their locations with a foreign key constraint. I can imagine that there may be reasons to split locations out into separate tables for points, lines, and polygons, however. If I can keep all locations' spatial data in a single column, then which geography datatype should I use (I know geometry would be more flexible, but that's a topic for another post)? I'm thinking geography(polygon,4326).


texts: Joining table to layer ArcMap without matching fields?

texts:
I know how to join things, but only when they both have matching fields.  I have created my own fields before as well, but only when I have had small datasets (less than 100), so it was practical for me to fill in each field by hand.

Now I have the parcelbase data for Saint Louis City (over 128,000 parcels), and the March 2017 crime .csv data from SLMPD (approx 4,000 crimes in March 2017).  Both the parcelbase and SLMPD data have addresses, but they have them in different formats (one has street and the other has street direction and street - ex: Field St vs N Field St), so I can't join them since they won't match up in Validate.  The second problem is the addresses are two fields each - the address number (ex, 5400) and the street name (ex, Field St), and I don't know how to join based on matching two fields instead of just one.

My boss thinks this isn't something that can be done with Join at all (and I agree), but neither of us know how to match them up otherwise.  I do have X Y coordinates for the SLMPD data (not for the parcelbase, but the parcelbase is projected to _1927 and the SLMPD is just an attribute table from a .csv file).

Is this a Georeferencing thing?

How do you do one of those?

If it is not a Georeferencing issue, does anybody know what it is that I need to do to combine these two things?

We want to visually map crime in Saint Louis City for March 2017.


texts: Sum values within one raster (25m) using another raster (1km) ArcGIS Desktop 10.2

texts:
I am summarizing some population data using a 1km raster overlay. I converted population data, polygon layer, into a 25-meter raster that contains the value for population/m^2. That value is then multiplied by 625 m^2 (the area of each cell), then I figured that I could do a zonal summary using my 1km raster layer (which is snapped nicely to the 25-meter raster layer), but now I am stuck. If I could just add up the values of each 25-meter cell that occurs in a 1-km cell then I would only need to multiply by 625 to get the approximate total population in that cell... I could use a majority resample but this should give a better estimate.

I got some great advice on coding a solution, but I'd like to figure it out along this line of thinking if possible.


texts: In ArcGIS 10.2 I created multi-segment lines using the Editor. Once a line was drawn, I selected "Finish Sketch" and then "Saved Edits". My new line appears in the attribute table but not on the map. However, if I double click said line in the attribute table, the line appears as selected in the map but disappears upon removing my selection. Why don't my edits appear permenantly on the map?


texts: Splitting string into "alphabetic" or "numeric" in QGIS 3.10.1

texts:
What are the right functions and formulas with the field calculator to delete "alphabet" or "number" inside a field?
Let me be clear. I want to be able to delete "alphabetic letters" or "numbers" inside my records (see the image). I've used replace, left, substr functions, but they don't make any difference between "alphabet" and "number".


texts: Location/format of inserted pictures used for marker symbols in ArcGIS Desktop?

texts:
I have an old mxd that uses several hundred custom made picture marker symbols.  We don't know where the image files that were inserted as pictures marker symbols are, or what format they were in.  We're trying to search our files to find said images but Arcmap (10.3) doesn't provide that information anywhere we can find.  Inserted pictures give you the path, but not so with picture marker symbols for some reason.

Is there any way to tell the file name and/or format of these inserted picture marker symbols?


texts: How to correct the elevation of the raster?

texts:
While comparing the DEMs of relatively small area raster (1) obtained of UAV data and (2) EU-DEM elevation model, I spotted around 30 meter difference in elevation between them.

Since no GCP was used while harvesting data and I know that commercial grade UAV GPS chips are not able to deliver accurate Z values, I checked my findings via elevationfinder.com which gave same results as EU-DEM.

To correct the data, I've established Z values difference for points which are easiest to identify on both rasters and recalculated it using Raster Calculator. I am not sure however whether this is the proper way to do it. Can somebody confirm or suggest other, more scientific and reliable approach?


texts: Shapefiles with same coordinate system do not match

texts:
I have a problem with ArcGIS. I'd like to project two layers in the same coordinate system (WGS84/Mercator; 3395) but these don't match.

I have one map with sovereign states borders, drawn from Nüssli EurAtlas data. The spatial reference system for all layers in this dataset is the world Mercator projection of World Geodetic System (WGS84).

Additionally, I have data on pogroms with latitude and longitude, drawn from Koyama (https://onlinelibrary.wiley.com/doi/abs/10.1111/ecoj.12331).

Ideally, I'd like to be able to project the pogroms point layer and project these points on the Nüssli map. I have "implemented" two strategies to do so:

I started by projecting the map shapefile in WGS84/Mercator, and then projected the points (pogroms) layers, by adding data XY, specifying the latitude and longitude of the excel file, and finally defining 3395 as the coordinate system.
Given that this turned out not working, I tried something different. I projected both of these shapefiles in Europe Lambert Conformal Conic, before adding a new dataframe and then inserting the two shapefiles in Lambert coordinates. And the shapefiles keep not matching.

Here is a snapshot of the pogroms data.


texts: How to find Gaps between road segments in ArcGIS

texts: I have huge data and want to find out gaps between segments. Although I have used Near tool and found some results but not a clear picture.


texts:
I have created a mask with the NIR band to separate earth from the sea.

How can I delete little spot area (errors)  in the mask or change the  pixel value of these spot areas ?


texts: Inaccurate Result from using QGIS contour function on pgr_driving distance

texts:
I was wondering whether anybody else has experienced a similar problem when using the contour function in QGIS.

I've run PGRouting's pgr_drivingdistance function on some osm data, and applied the contour function to this, although this has produced an inaccurate output

I've used the clip function to clip the sections of highway that fall within this contour as shown below, however this is consequently inaccurate because of the the contouring.

Does anybody know of any solutions to this problem?

I'm using the contour plugin in QGIS, and these are my settings. I've played around with these settings but it doesn't seem to make any real difference to the outputs which I'm getting.


texts:
I would like to generate a .shp with five points using shapefile available through pyshp. However, when I loop through my list I get only the last value. This makes sense as I did not append any of the points to an empty list (e.g. new_shp = [])

import shapefile as sf
import os

filename2 = 'test/point10'
lis = [(33.21, -122.15, 'france'), (35.31, -122.15, 'germany'), (35.41, -123.15, 'Hawaii'), (30.51, -122.15, 'Philippines'),(32.30, -122.15, 'Texas')]
for l in lis:
    w = sf.Writer(sf.POINT)
    w.point(l[0], l[1])
    w.field('location')
    w.record(l[2], 'Point')
    w.save(filename2)

# create the PRJ file
prj = open("%s.prj" % filename2, "w")
epsg = 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]'
prj.write(epsg)
prj.close()

I am assumming that I have to append all the points somewhere in the for loop, but I don't understand how to actually append it. How do I effectively loop through lis to get a .shp file with all five points?


texts: Splitting up selected islands using ArcGIS for Desktop?

texts:
What
I'm trying to do is to split up that selected islands as single-part features/ polygons. But not all of them, you can see some little dots around the islands are like a part of the island. So I need to separate them by selecting. How can I do that?


texts: Simplifying road network for cartography?

texts:
How do I simplify and merge segments of a complex network for cartographic purposes?

As you can see, i have a lots of segments and i want to merge and simplify to looks neat on a scale of 1: 100 000. I use many arcgis tools to try to make it look neat but its never working reallly good....

I would like to get neat road segment like this picture. How can i do that?


texts: How to update position of markers with “real time gps coordinants”?

texts:
I've got real time gps positions for a few cars and I want to create a map with updating markers. My code works but it doesn't "update" the markers, instead it adds new objects with new coordinates to the leaflet map. After few minutes my map is full of markers. What I'm doing wrong? Here is my basic concept.
var intervalV = document.getElementById("intervalValue").value * 1000;
document.getElementById("setIntervalButton").onclick = startData;
function startData() {
DataInterval = window.setInterval(getNewData, intervalV);
};

function getNewData() {
    $.getJSON(server, {
        fun : "GetGpsData",
        userId : "user",
        sessionId : $("#sessionId").val()
    }, fillMap);
}
function fillMap(json) {
for (var i = 0; i < json.devicesData.length; i++) {
    var positions = json.devicesData[i].positions.length;
        var devicepostiion;
        if (json.devicesData[i].connected == false
        ) {
            var devicepostion = L.marker([json.devicesData[i].positions[positions - 1].lat, json.devicesData[i].positions[positions - 1].lon], {
                    icon : offlineCarIcon
                }, {
                    draggable : false
                }).addTo(map);
        } else {
            devicepostion = new L.marker(, {
                    icon : onlineCarIcon
                });            devicepostion.addTo(map).setLatLng([json.devicesData[i].positions[positions - 1].lat, json.devicesData[i].positions[positions - 1].lon]).update();

    }
}

}


texts: QGIS determine if line (layer 1) has a + or - offset to a line (layer 2)

texts:
I have a network of lines in layer 1 and a network of lines in layer 2, that always stay one side of the line in layer 1.
I'm wanting to essentially determine whether a line is left or right of the other based on the direction of digitization.


texts:
I've got a Shapefile containing contourlines which I would like to style differently according to their value. For that reason I created a field named "100s" where I would like to give values from 0=25m-lines, 1=100m-lines and 2=500m-lines.
For doing that I tried the "replace" function with this syntax:

NAME="4600" OR NAME="4700" OR NAME="4800"  replace('NULL','NULL','1')

but QGIS gives the following error: "syntax error, unexpected FUNCTION, expecting $end". When I then set "$end, the error message is the same. Any suggestions how to resolve that?
Thanks in advance,
Ralf


texts: Extracting polygons from closed lines in QGIS

texts: I have a line layer having lakes and rivers. I would like to extract only the lakes (which are closed curves) from the layer.


texts: Get xy coordinates of QgsMapCanvas in PyQGIS 3?

texts:
I am trying to get the QgsMapCanvas coordinates when they change using PyQt5 and PyQGIS 3. What I have used for that:

self.canvas.scaleChanged.connect(lambda: self.showScale(self.canvas.scale()))
self.canvas.xyCoordinates.connect(lambda: self.showXY(self.canvas.xyCoordinates()))

and functions definitions are

def showXY(self, p):
    # SLOT. Show coordinates
    self.lblXY.setText(str(p.x()) + " | " + str(p.y()))

def showScale(self, scale):
    # SLOT. Show scale
    self.lblScale.setText("Scale 1:" + str(scale))

I am able to get the scale but when the xyCoordiantes is connected to the slot, window goes unresponsive and force closes.
What am I doing wrong?


texts: spplot() does not work properly when called via a function

texts:
I have spatial data in R, let's say an object of formal class SpatialPointsDataFrame but the below problem is encountered with Raster too.

I can plot easily the data with spplot() and save a .png file. If I want to create and store a function to recall this operation (since I have to repeat this many times in a simulation), spplot() save an empty .pgn file UNLESS I use dev.off() outside the function.

# example data
trees <- data.frame(x=runif(100, 1,100), y=runif(100, 1,100), value=runif(100, 20,60))
trees <- SpatialPointsDataFrame(trees[,c("x","y")], trees)
proj4string(trees) <- CRS("+proj=utm +zone=30 ellps=WGS84")

# simple plotting, it works
png(filename="spplot0.png", width = 480, height = 480, units = "px")
spplot(trees, "value")
dev.off()

# create functions
map1 <- function(){
  png(filename="spplot1.png", width = 480, height = 480, units = "px")
  spplot(trees, "value")
  dev.off()
}

map2 <- function(){
  png(filename="spplot2.png", width = 480, height = 480, units = "px")
  spplot(trees, "value")
}

# carry out functions. map1() save an empty image.
map1()
map2();dev.off()

Why map1() does not work and map2()+dev.off() outside does?

I know there are alternative mapping options (plot, spplot, ggplot) but I wish to stick with spplot() for this job. I have already found a work-around solution (see function map2) but it baffles me, my code gets more cluttered and I fear to encounter more problems in the future.


texts: Change QGIS language instalated in Ubuntu 16.04

texts: I'm begginer Ubuntu 16.04. I installed QGIS following steps I found in youtube (https://www.youtube.com/watch?v=pEAdKmToULo). QGIS was installed in Portuguese. I want change QGIS language to English and if possible, I want be able to change QGIS language by clicks. It is possible?


texts: What does QGIS do with NaN values in a GeoTIFF?

texts:
I have a number of 32-bit-float single-band GeoTIFFs where about 30% of the values are NaN. I have confirmed the presence of these NaN values by reading the GeoTIFFs using GDAL and libTIFF. However, when I open these GeoTIFFs in QGIS as raster height layers, all the NaN values are somehow replaced with what appears to be valid ground heights.

Why is QGIS doing this and is there a way to force QGIS to load the raw values from a GeoTIFF?


texts: Arcpy: allow user input for predefined filter list of strings

texts:
I've created a small toolbox where I use a GPString input with a filter list:

param1 = arcpy.Parameter(
    displayName="OSM tag value",
    name="in_key",
    datatype="GPString",
    parameterType="Required",
    direction="Input",
    multiValue=True
)
param1.filter.list = self.getConfig("tags")

The function getConfig reads a JSON and fills the values of the String input. Yet I want to enable the user to use an own string for the input at param1.

Is there any chance to do this with the given approach or do I need to create a second input for a user defined GPString?

Here is the issue on GitHub.


texts: Reprojecting raster between laea and lon/lat alignment issues

texts:
This is a follow-up post to a recent question on StackOverflow. I'm now facing a GIS problem rather than a coding one.

I'm working with a DEM from this source, which has the following projection information:
+proj=laea +lat_0=-100 +lon_0=6370997 +x_0=45 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0.

I need to projectRaster(mydem, worldclim) to fit worldclim data from this source, which has the following projection information:
+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0.

However, when I reproject mydem using worldclim as a template, I run into major alignment issues. mydem is a map of North America in Lambert azumithal equal-area and worldclim is the full globe in lat/lon, so I'm wondering if there is something going on extent-wise that might be causing the problem.

Here is an image of the original mydem:

Here is an image of the original worldclim layer:

Here is an image of the reprojected DEM, calculated using

mydem.lonlat <- projectRaster(mydem, worldclim)

And here are the mydem.lonlat and worldclim layers overlain with one another:

I followed Dr. Hijman's approach in this post but I don't think rotation is the issue here. There are a few other related questions here and here but they are unanswered and still don't address my exact issue.

extent(mydem.lonlat) and extent(worldclim) are identical, as are res() and proj4string() for both layers... so they definitely align in plotting space; there's just clearly an issue with geographic space.


texts: How add two coverages with different bounding-box and different resolution using geotools

texts:
I'm developing a Java application with some spatial analysis operations. One of the operations is to add two coverages. But how can I add two coverages with different bounding-box and different resolution, using GeoTools?

Any idea?
Is there any sample where I can start?


texts: Incomplete intersection area using gIntersection in R

texts:
I have two SpatialPolygonsDataFrame, one of them is the deforested area in the Amazon and the other is the city of Vila Rica, which is in the same region.
I would like to calculate the total deforested area in Vila Rica, so I'm trying to delimit the deforested area inside the city.

I used the gIntersection function in rgeos package to compare these two SpatialPolygonsDataFrame, generating the following result.

Green is the deforested area in the Amazon and blue is the area defined by gIntersection function.

Apparently the function worked well with smaller polygons, identifying only areas inside the city, but for some reason didn't work for the largest polygon, which crosses the city.
I don't know why the function didn't identify the deforested area inside the city for the largest polygon, perhaps a problem in the function?

I've tried the Over method in SP package but it doesn't seem to work like gIntersection.

The shapefiles (Rdata extension) are available in https://www.dropbox.com/sh/1jgn3b3k8hviw1i/AABGBS7zAzoI2YiYPsMisbApa?dl=0.

The "cities" data refers to all the cities in the Amazon, not only for Vila Rica.
In reality my main goal is calculate the total deforested area in each city of the region. I used only Vila Rica in the example to simplify.

Here is the R code for the map:

library(rgdal);library(sp)
library(rgeos)

### Download Rdatas and save in workdirectory ########
setwd(<select wd>)

load("deforested.Rdata")
load("cities.Rdata")

proj4string(deforested) <- proj4string(cities)

### Subset Vila Rica from cities shape ########
VilaRica <- subset (cities, nome=='Vila Rica')

### Delimit the deforested areas which are included in the city ###
deforestedVR  <- over(deforested,VilaRica)
deforestedVR  <- deforested[!is.na(deforestedVR),]

### Selecting the intersections ###
intersection <- gIntersection(deforestedVR, VilaRica, byid=TRUE,drop_lower_td=TRUE)

### Plotting the map ###
plot(deforestedVR,col='green')
plot(VilaRica,lwd=2,add=TRUE)
plot(intersection,col='blue',add=TRUE)

I was wondering if the polygon class may be affecting the function, since for the simplest examples it works perfectly.

I separated the largest polygon and changed its class from SpatialPolygonsDataFrame to SpatialPolygons and did the same thing to the Vila Rica shape.
Using again the gIntersection function, now works perfectly as you can see in the map:

Here is the R code for the second map:

### Identify and subset the biggest polygon ###
nPolys <- sapply(deforestedVR@polygons, function(x)length(x@Polygons))
bigPolygon <- deforestedVR[which(nPolys==max(nPolys)),]

### Changing class ###
VilaRica <- SpatialPolygons(VilaRica@polygons[1])
bigPolygon <- SpatialPolygons(bigPolygon@polygons[1])

### Selecting the intersections ###
intersection2 <- gIntersection(bigPolygon,VilaRica)

### Plotting the map ###
plot(bigPolygon,col="green")
plot(VilaRica,lwd=2,add=TRUE)
plot(intersection2,col='blue',add=TRUE)

Why for the same polygon the function didn't work for SpatialPolygonsDataFrame class but worked for SpatialPolygons class?
Hoping you all can help. Accept other suggestions.


texts: Join only if one line line intersects buffer. Do not if two

texts:
I've got a road network (line features) and a trip of GPS-Points. I want to join attributes from the road network to every point. So for each point i create a buffer. If the buffer intersects with a line feature, the attribtes are written to the point. That is all no problem. However I do not want to write attributes to points if the buffer intersects with to or more lines. Any ideas how to do this?
I work in ArcGIS. Is there a way to do this without using python?


texts: How to hide Openlayers 4 Overlay while hover over controls

texts:
I use sth like

map.on('pointermove', function (evt) {
    if (evt.dragging /*or evt.target === control*/) {
        tooltipPopup.setPosition(undefined);
    }

    //get all polygon features and show info about them in popup
};

But when feature is behind control, popup overlaps control. It is still clickable but it doesn't look nice. So my intent is to prevent displaying popup when i hover over control similarly like above for dragging.


texts: QGIS 2.14 LTR Merge Raster repeatedly fails

texts:
I recently installed QGIS 2.14 LTR and I am seeing what it can do.  I have repeatedly tried to merge 2 rasters (USGS 1/3 arc second DEM's) in both an ESRI grid format and a GeoTIFF format and each time Merge fails.  Here are the diagnostics:

gdal_merge.bat -of GTiff -o C:/Users/Steve/Documents/QGIS/Projects/Garmin/Data/GDD/grdn40-41w106.tif C:\Users\Steve\Documents\QGIS\Projects\Garmin\Data\GDD\grdn40w106.tif C:\Users\Steve\Documents\QGIS\Projects\Garmin\Data\GDD\grdn41w106.tif

Traceback (most recent call last):
File "C:\PROGRA~2\QGIS2~1.14\bin\gdal_merge.py", line 540, in
sys.exit(main())
File "C:\PROGRA~2\QGIS2~1.14\bin\gdal_merge.py", line 526, in main
fi.copy_into( t_fh, band, band, nodata )
File "C:\PROGRA~2\QGIS2~1.14\bin\gdal_merge.py", line 270, in copy_into
nodata_arg )
File "C:\PROGRA~2\QGIS2~1.14\bin\gdal_merge.py", line 77, in raster_copy
m_band )
File "C:\PROGRA~2\QGIS2~1.14\bin\gdal_merge.py", line 128, in raster_copy_with_mask
data_dst = t_band.ReadAsArray( t_xoff, t_yoff, t_xsize, t_ysize )
File "C:\PROGRA~2\QGIS2~1.14\apps\Python27\lib\site-packages\osgeo\gdal.py", line 2320, in ReadAsArray
callback_data = callback_data)
File "C:\PROGRA~2\QGIS2~1.14\apps\Python27\lib\site-packages\osgeo\gdal_array.py", line 330, in BandReadAsArray
buf_obj = numpy.empty([buf_ysize,buf_xsize], dtype = typecode)
MemoryError

AND using ESRI grids:

gdal_merge.bat -of GTiff -o C:/Users/Steve/Documents/QGIS/Projects/Garmin/Data/GDD/grdn40-41w106.tif C:\Users\Steve\Documents\DataGIS\USGS\DEM_1_3as\n40w106\grdn40w106_13\w001001.adF  C:\Users\Steve\Documents\DataGIS\USGS\DEM_1_3as\n41w106\grdn41w106_13\w001001.adf

Traceback (most recent call last):
File "C:\PROGRA~2\QGIS2~1.14\bin\gdal_merge.py", line 540, in
sys.exit(main())
File "C:\PROGRA~2\QGIS2~1.14\bin\gdal_merge.py", line 526, in main
fi.copy_into( t_fh, band, band, nodata )
File "C:\PROGRA~2\QGIS2~1.14\bin\gdal_merge.py", line 270, in copy_into
nodata_arg )
File "C:\PROGRA~2\QGIS2~1.14\bin\gdal_merge.py", line 77, in raster_copy
m_band )
File "C:\PROGRA~2\QGIS2~1.14\bin\gdal_merge.py", line 128, in raster_copy_with_mask
data_dst = t_band.ReadAsArray( t_xoff, t_yoff, t_xsize, t_ysize )
File "C:\PROGRA~2\QGIS2~1.14\apps\Python27\lib\site-packages\osgeo\gdal.py", line 2320, in ReadAsArray
callback_data = callback_data)
File "C:\PROGRA~2\QGIS2~1.14\apps\Python27\lib\site-packages\osgeo\gdal_array.py", line 330, in BandReadAsArray
buf_obj = numpy.empty([buf_ysize,buf_xsize], dtype = typecode)
MemoryError

Other GDAL algorithms (hillshade, contour) that I tried work OK but QGIS does not seem able to merge 2 rasters.


texts: Creating database connection using ArcCatalog

texts:
I'm using ArcCatalog 10.2.2 which is installed on a Windows 2019 Server OS.  I'm having trouble creating a database connection.
Usually this is a very straight forward process.

Put in the server name under instance.
Database authentication.
When going to select the database it can see the different database on that server, so I know it can connect.

So it all looks fine right.

When I then double click the connect to browse it, usually on other Windows10 systems it takes maybe 30-60 seconds to establish the connection and download tables.

But this just sits there thinking... and thinking... and thinking.. its been over half an hour.  So there is something wrong, but I don't know where to start debugging it.

Any suggestions?


texts: How do I Identify and flip the lines from source to the destination solution in ArcMap?

texts:
I have got cable network data and I wanted to know if there is any way to flip all the lines in such a way that start (FROM) point should be the source which is direction of current.

Red labels are cable numbers. Circled blue lines need to be flipped, they have got 'arrow at the end' symbology.

Is there any solution to flip all of the lines in such a way that they are carrying electricity from source to the next cable.  Possible to find network analyst functions in arcmap or something in python


texts: using arcgis server connection file(.ags) inside python is not running in linux

texts:
When i am trying to execute the following script from my desktop it's executing successfully, but same is not executing/not throwing error in linux

              agsPath="/sharedFolder/arcgiscon(admin).ags"
               try:
                  print "stared"
                  for item in arcpy.ListDataStoreItems(agsPath,'FOLDER'):
                      pint item
               excpet:
                      print "problem in executuion"

In linux it is printing stared but after putty is not showing any thing. My arcgis version is 10.2

What could be the reason.


texts: Kernel density estimation across multiple UTM zones

texts:
I'm trying to compare the spatial distribution between two groups of fishing locations (for each of several different years) in the Bering Sea (lat's ~53 to 65N, long's ~ -178 to -158W. My preliminary explorations have used kernel density estimation (with an unconstrained plug-in bandwidth matrix) and then examined different overlap indices of their home ranges (quantified via utilization distribution overlap index and Bhattacharyya's Affinity) to compare similarities/differences (all performed in R).

However, I'm concerned about the interpretation of these data when the data may span more than one UTM zone. Most of my data are in UTM zone 3 with some data points seeping into zone 4. In more extreme cases, my data may span from zone 1 to 4.

Am I being really ignorant in thinking that I can compare across zone boundaries?

Since I am interested only in relative differences between calculations with the same distortion, will it not matter?

Should I simply make sure that I have my UTM origin set at the southwesterly most point of the data and I'm fine from there?

I've been reading Stack Exchange and other sites all morning and have seen several discussions about measuring across UTM zones, but I haven't been able to logically extend such responses to a method such as this one.


texts: How to exclude one specific polygon from a layer which have like 20 diffirent polygons?

texts:
I have a layer with like 20 polygons which I imported from Google Earth Pro.

And what I am trying to do is to exclude (isolate) that one big polygon which is in a way of all other smaller ones so they are not visible under big one.

I know that one way is to make it transparent, but while doing that all other polygons will become transparent aswell.

So what i need is like a tool or something which will help me to "extract" only one specific polygon from all those in a same layer, into another layer, and after that i can transform it into lines (Vectors) so there will only left boundarys of that specific polygon which can be seen on picture, or even make it transparent.


texts:
I have this set up: "POSTGIS="2.1.2 r12389" GEOS="3.3.8-CAPI-1.7.8" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.9.2, released 2012/10/08" LIBXML="2.7.8" LIBJSON="UNKNOWN" TOPOLOGY RASTER" on PostgreSQL-9.3

I installed CGAL 4.4 and SFCGAL

~# sfcgal-config --version
1.0.4

I would like to use 3D functions (ST_Extrude... etc) in my exiting databases, is it possible? I can't find proper information about this topic.


texts: How to access user credentials in database connection?

texts:
I am trying to automate a couple of things when a QGIS plugin is activated:

If Layer A is not in the TOC, add it. This requests username and password to access a PostGIS table using the inbuilt credentials dialog.
If Layer A now exists in the TOC, access the same PostGIS table to run a query using the same user credentials.

Apparently as I'm not storing the password, I can't access it to use in connecting to the database (using psycopg). Is there any way around this? I want to avoid having the user type credentials twice, and also I want to avoid storing the credentials permanently. I have tried to reimplement the QgsCredentialDialog class but that adds so much complexity (e.g. when a connection error is thrown).

What other ways can I add a layer using a non-stored username and password, and access the credentials later to query the same database? The easiest option I can think of is to access the inbuilt credentials dialog request, but I'm not sure how to do that.

EDIT

My implementation of this was to create a function which I call at the point my layer has been added. It then sets my plugin class attributes username and password correctly which I can use in my db query.

def get_credentials(self, lyr):
    connInfo = QgsDataSourceURI(lyr.dataProvider().dataSourceUri()).connectionInfo()
    success, username, password = QgsCredentials.instance().get(connInfo, None, None)
    if success:
        # QgsCredentials.instance().put(connInfo, username, password)
        self.username = username
        self.password = password


texts: Getting user input for QGIS FieldCalculator parameters

texts:
I'm trying to calculate and output a virtual field -- something that is a straightforward function of fields from a vector layer in QGIS 3. The FieldCalculator does the trick, but I need to manually specify parameters (e.g. coefficients in the function).
That's easy via FieldCalculator, but if I want others to use this calculation, what's the easiest way to prompt the user to input parameters?
A plugin seems like overkill for something this simple.


texts: Is there a way to link CID number (Fema Community ID) with the Census GEOID?

texts: I have a list of areas from the Fema Community Rating System that includes a large amount of communities (which are indeed Census places) and US Counties. I need to link each area to a Census GEOID (FIPS — Federal Information Processing Standards) to be able to map them and run some analysis. Fema publishes a "Data Capture Standards Technical Reference" (p. 4-5) that indicates that a table linking CIDs to GEOIDs exist, but I couldn't find it anywhere. Using some brute force work I can try to match the names using R, but it is not ideal ... I could find about 90% of roughly 2000 areas, and I'm not sure how good is this combination since there is a ton of places with a very similar name.


texts: Which is more performant GeoJSON or WMS/WFS using OpenLayers?

texts:
I'm developing a web application using PostGIS, Laravel, Vue.js and Openlayers, for the moment I generate the GeoJSON from my database using the following code :

 public function plosWithGeom()
    {

        return DB::table("plos")->selectRaw("ST_AsGeoJSON(geometry) as geom,*")->get();
    }

which is equivalent to

 select ST_AsGeoJSON(geometry) as geom,* from plos

using this approach serves me a result that's easy to manipulate using JavaScript, but I'm asking about performance compared to serve this data using an OGC service (WMS/WFS ) on MapServer

Update :

The server sends the data to the frontend like a text, in the frontend I parse this data using JavaScript and show it with OpenLayers. For the moment I'm working with a region, the server serves linestring and points with 530KB as size and I predict that will reach 10MB when working with all the country.


texts: Refering to 0° of longitude/latitude in a range

texts:
This is more of a question in diction...

In scientific writing, is there a preferable way to write a range, bounded by 0° longitude, in a way similar to "5° W to 10° E" if there is no cardinal direction for 0.

The two ways that seem most reasonable to me are "0 to 5 degrees longitude" or "0° longitude to 10° E".

Thanks!


texts: Touch 'tap' event not available on openlayers 2.11 map (update: click works fine)

texts:
According to Will mousedown events fire with Touch events on OpenLayers mobile?, openlayers 2.11 supports click-like behaviour for vector layers only.

Separate from that, regular pinching/zooming/panning works on the whole map. Weird thing, though, is that zooming on the ipad triggers a click event...

My problem: I don't have a vector layer, but I do want to handle click events on an ipad. I tried attaching 'touchstart' in addition to 'click' to my existing click handler, but that doesn't give any reaction.

So my question is: can I have an OpenLayers.Handler.Click handler on the whole map that reacts to touch taps?

Update: click events work anyway, see below.


texts:
Even after rotating a GeoTIFF file, the resulting file size do not change. however, the rotated file overlays on the original TIFF file. When I tried  gdal_translate with projwin switch an error occured:

  ERROR 1: The -projwin option was used, but the geotransform is rotated.  This configuration is not supported.

I am using GDAL 2.4.


texts: GDAL Python, projecting WGS84 geom to UTM to find area

texts:
I'm writing code with GDAL in Python. I'm trying to find the area of a geom in WGS84 with GDAL. I tried projecting it to UTM and then getting the area of the projection, but the answer isn't what I expected. I have created a simple script to demonstrate.
The answer from GDAL is 22658.299678236246. My expected area by using Postgis is something around 13798. Is there something wrong with my GDAL code?
#!/usr/bin/env python3

from osgeo import ogr

def main():
    text_geom = 'MULTIPOLYGON(((-22.5346102413589 21.7484412769131,-22.5346153453775 21.7479038269593,-22.53656903248 21.748114448725,-22.5365549834289 21.7488859855299,-22.5362457983049 21.7487836361578,-22.5359903673763 21.7486875138083,-22.5358014126894 21.748679989781,-22.5357257154141 21.7486479216644,-22.5355649068295 21.7485918834607,-22.5354704464861 21.7485921574616,-22.535215143558 21.7485283202296,-22.5350166478364 21.7484966071143,-22.5346102413589 21.7484412769131)))'
    geom = ogr.CreateGeometryFromWkt(text_geom)

    print('Starting geom:')
    print(geom.ExportToWkt())
    print()

    point = get_geom_point(geom)
    print('point: ' + str(point))

    logitude = point[0]
    latitude = point[1]

    # Find UTM Zone using a point from the geom.
    utm_zone = 1 + int((logitude + 180.0) / 6.0)

    print('utm_zone: ' + str(utm_zone))

    if latitude > 0:
        # North UTM zone
        projection_number = 32600 + utm_zone
    else:
        # South UTM zone
        projection_number = 32700 + utm_zone

    print('projection_number: ' + str(projection_number))

    wgs84 = ogr.osr.SpatialReference()
    wgs84.ImportFromEPSG(4326)

    utm = ogr.osr.SpatialReference()
    utm.ImportFromEPSG(projection_number)

    print()
    print('Reprojecting from:')
    print(str(wgs84))
    print()
    print('Reprojecting to:')
    print(str(utm))
    print()

    transform = ogr.osr.CoordinateTransformation(wgs84, utm)

    geom.Transform(transform)

    print('Reprojection:')
    print(geom.ExportToWkt())
    print()

    print('Area: ' + str(geom.GetArea()))
    print('expected around 13798')

    # => SELECT ST_AREA(ST_geomfromtext('MULTIPOLYGON(((-22.5346102413589 21.7484412769131,-22.5346153453775 21.7479038269593,-22.53656903248 21.748114448725,-22.5365549834289 21.7488859855299,-22.5362457983049 21.7487836361578,-22.5359903673763 21.7486875138083,-22.5358014126894 21.748679989781,-22.5357257154141 21.7486479216644,-22.5355649068295 21.7485918834607,-22.5354704464861 21.7485921574616,-22.535215143558 21.7485283202296,-22.5350166478364 21.7484966071143,-22.5346102413589 21.7484412769131)))')::geography);
    #      st_area
    # ------------------
    #  13798.4696836025
    # (1 row)
    #
    # =>

def get_geom_point(geom):
    if geom.GetPointCount() > 0:
        return geom.GetPoint(0)

    for i in range(geom.GetGeometryCount()):
        point = get_geom_point(geom.GetGeometryRef(i))
        if point is not None:
            return point

    return None

main()


texts: Creating p:commandlink in Popups in OpenLayers?

texts:
How do I create p:commandlink in Popups can be called to a bean method like actionListener="#{offLineMapsTestBean.fetchInstList(variable)}"?

I tried with this code:

map.on('click', function(evt) {
    var feature = map.forEachFeatureAtPixel(evt.pixel, function(feature) {
        return feature;
    });
    if (feature) {
        var coordinates = feature.getGeometry().getCoordinates();
        popup.setPosition(coordinates);
        $(element).popover({
                'placement' : 'top',
                'html' : true,
                'content' : 'Name:'+ feature.get('name')+ '\n'                                            + "
 <button onclick= '#offLineMapsTestBean.fetchInstList(feature.get('name'))}'>  </button>
"
        });
        $(element).popover('show');
    } else {
        $(element).popover('destroy');
    }
});


texts: how does arcgis server manage tile caches with an unlimited scale choice?

texts:
I am familiar with tile caching with geowebcache, but not with ArcGIS server. I imagine the basic concept is the same. First determine how many scales with be cached, set the image quality and format, and off you go. I would then usually design my own web client and set the scale choice to those set in the tile cache.

Now I am using ArcGIS server and do not have the ability to set the client scales.
So the client users can will basically hardly ever use the caches set to 1:1000, 1:2500, 1:5000, 1:10000,1:2500 etc...they will always be somewhere in between these set scales.

Is this how it will work, or am I missing an important ESRI concept?


texts: Plotting heatmap/locations in GeoServer & OpenLayers map from Lat/Lon data of PostGIS table

texts:
I am using Geoserver, Openlayers, PostGIS for Geospatial web implemenation.

I have an Indian map shape files for States, Districts, Talukas. I am able to generate a map with states, districts, taluka using WMS layers in Openlayers from Geoserver.
However, now on top of that layers I have to show state, district, taluka with some counts & locations which I am have in PostGIS table. The PostGIS table contains state, district, latitude and longitude of the location, as well as some counts which I need to show in the map generated in (1#).

How can I implement it?

I also have to provide a Drilldown menu on the map, using filters based on statewise, districtwise data.

How do I proceed further?


texts: QGIS Server Internal Server Error - OS 109) pipe has been ended. : mod_fcgid: get overlap result error

texts:
I have installed QGIS Server using the QGIS Server on Windows7 Step-by-step instructions and when I finished and tried doing http://localhost:83/qgis/qgis_mapserv.fcgi.exe?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities I get an error:
Internal Server Error
The server encountered an internal error or misconfiguration and was unable to complete your request.
Please contact the server administrator, admin@local.domain and inform them of the time the error occurred, and anything you might have done that may have caused the error.
More information about this error may be available in the server error log

The log shows:
 [Wed Dec 12 22:25:19 2012] [notice] Apache/2.2.14 (Win32) mod_fcgid/2.3.6 configured  resuming normal operations
[Wed Dec 12 22:25:19 2012] [notice] Server built: Sep 28 2009 22:41:08
[Wed Dec 12 22:25:20 2012] [notice] Parent: Created child process 3868
[Wed Dec 12 22:25:20 2012] [notice] Child 3868: Child process is running
[Wed Dec 12 22:25:20 2012] [notice] Child 3868: Acquired the start mutex.
[Wed Dec 12 22:25:20 2012] [notice] Child 3868: Starting 64 worker threads.
[Wed Dec 12 22:25:21 2012] [notice] Child 3868: Starting thread to listen on port 83.
[Wed Dec 12 22:26:08 2012] [warn] [client 127.0.0.1] (OS 109)The pipe has been ended.  : mod_fcgid: get overlap result error
[Wed Dec 12 22:26:08 2012] [error] [client 127.0.0.1] Premature end of script headers: qgis_mapserv.fcgi.exe
[Wed Dec 12 22:26:09 2012] [error] [client 127.0.0.1] File does not exist: C:/OSGeo4W/apache/htdocs/favicon.ico

What could be the problem ?


texts: How do I write out a mixed geometry geodataframe to a geopackage?

texts:
I have a geodataframe of several thousand rows, mixed geometry (Point, LineString, Polygon).  I would like to write this out to a geopackage with different layers per geometry type.  Here's what I've got so far (in jupyter nb):
import geopandas as gpd
%matplotlib inline

geojson_file = "C:/Temp/data.geojson"

gdf = gpd.read_file(geojson_file)

gdf.head() # yields good results
gdf.plot() # draws all the various geometries where they should be

gpkg = "C:/Temp/converted.gpkg"
gdf.to_file(gpkg, driver="GPKG")

Then when I load the layers...
import fiona
for ln in fiona.listlayers(gpkg):
    with fiona.open(gpkg, layer=ln) as dataset:
        print(ln, len(dataset))

it only gives me
converted 10000

Where am I going wrong?
EDIT: I'm on Windows 10, Python 3.8 (anaconda), and GeoPandas 0.8.1


texts: How to find the coordinate system of an existing Map in arcGIS using VB.NET

texts: I already have a georeferenced base map loaded in ArcMap. I want to find the spatial reference of this map and assign the same coordinate system to a polygon that I created. I can not get a start with the code as I am a beginner.


texts: Calculate Tiles X,Y, Z for current view or extent in OpenLayers

texts:
Im using OpenLayers with the openstreetmap tile server.

var urlTemplate = "http://c.tile.openstreetmap.org/{z}/{x}/{y}.png"

I would like to download tiles for the current view or an extent. So that I can view them later offline.

Therefore Im trying to calculate the z, x,y for the tiles

The coordinate Z should be avaible from the maps zoom

But how can I find the different x,y Coordinates for the view/extent?  Is there a OpenLayers funtion for that?


texts: Call ArcPy.arcpy.StripMapIndexFeatures_cartography for specific feature within layer?

texts:
I'm creating a custom python script, and having trouble calling the script on a specific route within a layer

routekeyfield = 'route_ID'
routename = 'I95 NB'

routeslayer = "C:\\Projects\\\\Routes2017.mdb\\Routes"
# table = arcpy.CreateTable_management("in_memory", "stripMaps")
#memoryFeature = "in_memory\\stripMaps"
memoryFeature = "c:\\path\\stripmaps.mdb\stripmapsstuff"

#query the routes layer for the specific route
fields = ["SHAPE", routekeyfield]
whereclause = routekeyfield + " = '" + routename + "'"
with arcpy.da.SearchCursor(routeslayer, fields) as cursor:
    for row in cursor:
        arcpy.StripMapIndexFeatures_cartography(in_features=row, out_feature_class=memoryFeature, use_page_unit="NO_USEPAGEUNIT", scale=str(currentextentscale), length_along_line=str(currentextentwidth) + " Meters", length_perpendicular_to_line=str(currentextentscale) + " Meters", page_orientation="HORIZONTAL", overlap_percentage="10", starting_page_number="1", direction_type="WE_NS")

I'm not exactly sure how to filter the layer for a specfic route. Note, I'm using this script from a standalone script. Anybody have some ideas?
`


texts: How to use Geoserver Authentication with Spring Security in WebApplication?

texts: I want to implement that Authentic login users of my website can only have access to Geoserver's specific layers.


texts: How to classify a raster using QGIS Plugin interacting with PostGIS?

texts:
I'm developing a QGIS plugin where users can enter some parameters (like distance, depth, area...) and by submiting the plugin interact with a PostGIS database and extract informations form vector and raster layers.
I want to classify a DEM raster by substituting the pixel values, for example:

If  200 < val < 210    Then   val = 205
If 211 < val < 220    Then   val = 215

This classification is based on the value entered by users using a TextBox or LineEdit
I'm using the query bellow and it's working but manually not automatically:
This query took 2 seconds for a DEM raster with (X:6000 , Y:5000) and (MIN(val):200,MAX(val):2400)

So each time I need to specify the gap.
Is there any way to allow to users to specify the interval they want to be applied for all values of raster?


texts: Create View error for pgRouting in Mapserver

texts:
I'm trying to implement pgrouting with MapServer. It works right now but I have to enter the number of the source and target nodes in my mapfile by myself. What I want to do, is being able to use a click on the map and select start and end points.
So the tutorial I'm following suggests that I create a VIEW in postgresql which is the following:

CREATE VIEW shortest_path AS
SELECT
 min(r.seq) AS seq,
 e.old_id AS id,
 e.street,
 e.district,
 sum(e.time) AS time,
 sum(e.distance) AS distance,
ST_Collect(e.geom) AS geom
FROM
  pgr_dijkstra(
'SELECT
id,
 source::INT4,
target::INT4,
%cost% AS cost,
FROM roads_noded', %source%, %target%, false) AS r,
 roads_noded AS e
WHERE
  r.id2 = e.id
GROUP BY
  e.old_id, e.street, e.district;

but I get the following error:

ERROR:  syntax error at or near "%"
 LINE 17:   FROM roads_noded', %source%, %target%, false) AS r,
                          ^
  ********** Error **********

 ERROR: syntax error at or near "%"
 SQL state: 42601
 Character: 302

In the beginning I assumed it because of the percentage symbol so I put %source% in quotes like '%source%'. But then I got another error:

ERROR: function pgr_dijkstra(unknown, unknown, unknown, boolean) is not unique
 SQL state: 42725
 Hint: Could not choose a best candidate function. You might need to add    explicit type casts.
Character: 192

So, I'm quite new to this. Am I going the right way? What is wrong here?


texts:
I have been trying this for a while without success. After installing the OpenLayers module in drupal 7.14 and trying to configure it i get the following error:

Cannot establish a connection to the server:

I have tried all permutations including downloading the openlayers library in the local machine, giving it a local path in a number of different ways.

The second error that is seen suggests that i am installing an incorrect version of the openlayers library when in fact i have been installing the recommended version i.e. version 2.11.

It has been a frustrating experience getting Openlayers, Geoserver , PostGIS and Drupal 7 working together on Windows 7. The number of combinations are many and there are many issues that are difficult to solve if you are not a webmaster.


texts: Are Google satellite images precise enough for georeferencing?

texts:
I need to georeference historical aerial images (from 1976). There is no metadata, e.g. reference system, available. I thought I could use Google satellite to set ground control points for georeferencing.

However, I am not sure if this is precise enough. Has anyone experience on the quality of georeferenced images via Google?


texts: Adding a search box to a Leaflet.js example

texts:
Consider the example in the following link for Folium: US Employment visualization per county.

The result looks like what you see below:

I would like to add a search box to it to be able to search by city, and have the visualization zoom in and/or center around it. How can I go about adding such search box to it? Can I add this functionality perhaps in the HTML file?


texts: Zonal statistics USFS tree cover by census block

texts:
I want average tree cover for every census block. I have a file of census blocks, and I have a USFS tree cover file from USFS Tree Canopy Cover Datasets that looks like this:

I'm trying to use zonal statistics but the files are too big- even though my computer is very powerful QGIS only uses <3GB of RAM for some reason.
Is there a way to speed this up?
I tried using R's exact_extractr but I got all NaNs.
rasta = raster('canopy.tif')
poly3 = st_read("cbg.shp")
foliage <- exact_extract(rasta, poly3, 'mean')


texts: Creating Select Then Zoom Tool using ArcPy?

texts:
I have been trying to create a tool to make a simple select then zoom to selected but I've been running into some issues. I've been using the code from this original posting (https://community.esri.com/thread/217697-select-layer-by-attribute-python-script-tool). It seems like his solution was never solved.

Here is the code that I have:

import arcpy

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
lyr = arcpy.mapping.ListLayers(mxd, "PARCELS", df)[0]
arcpy.env.overwriteOutput = True
arcpy.env.workspace = "C:/ArcMap/COGO_Template"

wrk = arcpy.env.workspace
arcpy.AddMessage("Workspace added")
FClass = arcpy.GetParameterAsText(0)
Field = arcpy.GetParameterAsText(1)
Feature = arcpy.GetParameterAsText(2)

FClassLayer = arcpy.MakeFeatureLayer_management(FClass, "FclassLayer")

where_clause = """{} = '{}'""".format(arcpy.AddFieldDelimiters(FClass, Field),
                                      Feature)  # I have added extra pair of single quotes for selecting strings
arcpy.AddMessage(where_clause)
arcpy.SelectLayerByAttribute_management(FClassLayer, "NEW_SELECTION", where_clause)
arcpy.CopyFeatures_management(FClass, "C:/ArcMap/COGO.gdb/SelectionSites")
df.zoomToSelectedFeatures()
arcpy.RefreshActiveView()
arcpy.AddMessage("All done!")

My parameters are as follows:

InputFeatureClass = Data Type: Feature Class, Required, Input, No Multivalue, Default: Database Connections\Connection to gis.sde\PARCELS

InputField = Data Type: Field, Required, Input, No Multivalue, Obtained from Input Feature Class

InputValue = DataType: String, Required, No Multivalue

The tool runs fine and shows no errors but unfortunately no parcel is selected and it doesn't zoom. Something to note as well is it takes about 2 minutes to run.

When I manually select by attribute in the attribute table (select, select by attribute, put in the parcel number, click ok), then run the tool, it zooms to the layer


texts: Exporting several images from same collection in one task using Google Earth Engine

texts: I need to export several images (around 60) using Google Earth Engine. They are from the same collection and for the same geometry, but have different dates. I already have a list of the dates I need, but no matter what I do, I cannot seem to export them in one task.


texts: Keeping color of shapefile constant in ArcGIS Desktop?

texts:
I am using ArcGIS Desktop 9.3 when I use a color for the shapefile it changes everytime when I close and open ArcMap.

How I can constantly keep the color of the shapefile that I work with?