texts:
I already use another library "shapefile", which was easy to import. However I just installed shapely 1.3.3 for windows (which still has an installer) and cannot import it. I tried it the same way as with the other library (import shapely) after the installation, but it only throws the error:

ImportError: No module named shapely

Edit: I run Python-2.7.9 and have GEOS installed. If I try "import shapely" I get the error "No module named shapely" and if I don't do it, I get this:

Point=(0,0)
print Point.geom_type

AttributeError: 'tuple' object has no attribute 'geom_type'

Can somebody resolve this? It should be a legit command. I do not understand what could be wrong.

Edit: It seems like shapely is not recognized by Python, eventhough it definitely is installed. I tried using another version of Shapely and got the same result. I really need a Shapely function to continue working on a work project but cannot figure out why it doesn't work. I also tried installing it with pip, which wasn't of any use as well.


texts: Find arc endpoint with bearing and radius counting with curvature of the earth?

texts:
So my problem is that I have start point [lat,long] for arch, bearing, radius and I need compute end point [lat,long] ? End point should be in 180 degrees angle from start point and I should be aware of curvature of the earth ... So let say I have start point [54.01,104.05], bearing is 170 and radius is 2 km ... Any ideas?

Thanks in advance.


texts:
I am using sharpmap to render my OPEN STREET MAPS (osm) maps. i need to show the direction of one way streets on my map.
How can i do that?!
After googling i found SharpMap.Rendering.Symbolizer and style.Line.DashPattern might help but i couldn't figure it how do it...!
i have the oneway column in my data, and i can filter my data for one way streets for styling, but i don't know how to style and draw a small arrow on my lines.

i want the arrows like this picture Street Directions


texts: QGIS2WEB white screen on export/preview?

texts: When I export my QGIS2WEB project it goes through, but when it loads up in my browser after the export is done all I see is a white background with the zoom/minus feature in the top left hand of the screen.


texts: ST_DISTANCE return not same distance as QGIS

texts:
I am searching distance between point and nearest road (OSM data). I wrote next SQL request:

SELECT osm_id, st_astext(geometry),  ST_DISTANCE(geometry::geography, ST_SetSRID(ST_POINT(37.72308, 55.47957), 4326), true) AS dist, roads.geometry::geography <-> ST_SetSRID(ST_MakePoint(37.72308, 55.47957),4326)  as op_dist
FROM roads ORDER BY ST_SetSRID(roads.geometry, 4326)::geography <-> ST_SetSRID(ST_MakePoint(37.72308, 55.47957),4326) ASC LIMIT 1;

result dist is: 27.830223659

I added this data to the QGIS and calculated distance with the meter instrument and it's show that the distance is 34.3 meters.

What can be wrong and why?

PostGIS 2.2


texts: Seeking Farm Resource Regions shapefiles?

texts:
Could someone direct me to the potential sources to find the Shapefiles of farm resource regions?

The map exactly looks like the one here https://www.usda.gov/media/blog/2014/04/22/celebrating-our-glorious-planet

but I have had a hard time finding the Shapefiles.


texts:
If I have my data an ArcGIS Server and I don't want to bring all the data, but I want to make a query first to select group some data together then bring only the result. How can I make that as the records number is huge and I don't want the records, I want only the result of grouping them together?

For example: If you have lands, each assigned with for a specific use (agricultural, industrial, etc..), the number of records is huge, you have a website that filters the lands with their location, use, or both, and the result is only the the percentage area of the filtered data w.r.t the total area. So you don't want to bring all the data you want only the result. How to make this?


texts: Finding closest points on another polygon

texts:
I have three layers:

polygons (having a unique id)
stops (having the polygon id it is laying on) (point)
possibilities (having the polygon id it is laying on) (point)

I now want to find the closest possibility-point of each stop-point, which lays on another polygon.

I could easily find the closest point overall, but how could I add the requirement that this point has a different polygon-id than the stop-point itself has?

Here is an example of what I have:

And what I would like to achieve:

I could also use street data from OSM (laying between the buildings) if that would help.

Preferably I would like to be able to integrate this process in a graphical model in QGIS 3.4. But if that's not possible, a different solution than processing tools (like plugins or virtual layers) would do it as well.


texts: When I use the function .toGeoJSON(); is there way to keep the layer leaflet_id?


texts: Calculating distance between two points using OpenLayers

texts:
I don't understand why I have two different results between
Openlayers distanceTo () and Movable Type Scripts

LatLng : (43, 3)
LatLng : (42, 4)
OL :  187.6
Movable Type Scripts : 138.1

var Geographic  = new OpenLayers.Projection("EPSG:4326");
var Mercator = new OpenLayers.Projection("EPSG:900913");
var point1 = new OpenLayers.Geometry.Point(3, 43).transform(Geographic, Mercator);
var point2 = new OpenLayers.Geometry.Point(4, 42).transform(Geographic, Mercator);
console.log(point1.distanceTo(point2)/1000);


texts: Adding vertex into polyline at interval using ArcPy

texts:
I am using ArcGIS 10.2 version and I have polyline features, for this I have to add vertex at some distance (10mtr) and label the z value of that vertex.
Is there any tool or script available?


texts: Files locked after running a Python script

texts:
I can run a Python script that uses qgsVectorLayer, qgsRasterLayer and many functions like processing:runalg("qgis.difference"......) and os.popen("gdalwarp -q -cutline.......)
but if I want to run it another time, I must close QGIS first and reload it (error: file locked because used by other programs...)
I want to close files handlers but I don't know which ones.
What can I do in this case?

Here is my some of my code


texts: Appending and Replacing some (but not all) Fields using ArcGIS Pro

texts:
I am using ArcGIS Pro 2.4.2.

My question is similar to Maintain Feature Data without Creating It.

I am looking for a way to update our proprietary parcel layer with a shapefile downloaded from our county source. A truncate and append would work, but there are some fields that are not maintained by the county that are parcel specific (e.g., variances).

Is there a way to replace the updated fields while retaining selected existing fields, or do I need to maintain variances as a separate feature class that I spatial join to the append and truncated parcel layer?


texts:
The following script fails, if I import the 'rasterstats' library.

import rasterstats
import ogr

geom1 = ogr.CreateGeometryFromWkt('POLYGON ((415924.806155724043492 4489457.402624096721411,416359.905146608769428 4489457.402624096721411,416359.905146608769428 4488983.206433195620775,415924.806155724043492 4488983.206433195620775,415924.806155724043492 4489457.402624096721411))')
geom2 = ogr.CreateGeometryFromWkt('POLYGON ((415924.806155724043492 4488989.295723957009614,416795.004137493495364 4489457.402624096721411,416784.428755077067763 4488509.010242293588817,415924.806155724043492 4488989.295723957009614))')

print geom1.Intersection(geom2)

The error I get is:

Assertion failed: (pts), function testInvariant, file ../../include/geos/geomgraph/Edge.h, line 88.

The line causing the problem is geom1.Intersection(geom2)

Unfortunately there is not more error code to it. If I don't import rasterstats, then I don't get the error.
 What could be the reason for that behavior?

Line 88 in Edge.h in GEOS is the following function:

void testInvariant() const {
    assert(pts);
    assert(pts->size() > 1);
}

What does this function test for?

GDAL Version
1.10.0100

Rasterstats Version
0.3.3

GEOS Version
3.4.2-3

Edits:

Originally I thought it is related to Flask, but I was able to isolate the problem, and figured it is not related to Flask.

OGR exceptions doesn't give any additional error.


texts:
I want to set the heap for Geonetwork and as the doc says, I need to add this line to the start-geonetwork.sh script :
java -Xms1g -Xmx1g -Xss2M -XX:MaxPermSize=128m -Dgeonetwork.dir=/app/geonetwork_data_dir -Dgeonetwork.lucene.dir=/ssd/geonetwork_lucene_dir

Yet I can't find where this file is located, knowing that Geonetwork is deployed under tomcat.
Can anyone help me find this file?


texts:
I have two polylines whose endpoints are near to each other, but do not overlap. I've put them into this shapefile. I want to use shapely's ops.snap() to snap their endpoints, but it's not working as expected.

from shapely import ops
import geopandas as gpd

wontsnap = gpd.read_file(path_to_wontsnap)
g1 = wontsnap.geometry.values[0]
g2 = wontsnap.geometry.values[1]

If I check the distance between the geometries, I get

g1.distance(g2)
Out[661]: 0.0006421029901831172

This agrees (roughly) with the measure tool in QGIS:

Now if I try to snap the two geometries with a tolerance that is greater than the distance between their endpoints:

snapped = ops.snap(g1, g2, 0.001)

The snapped object contains only the g1 geometry. I.e., they don't snap.

If I increase the tolerance by 10x:

snapped = ops.snap(g1, g2, 0.01)

The snapped object includes both geometries, snapped correctly.

If I trim the g1 and g2 geometries to the last and first 10 points, respectively, snapping works as expected.

I have tried various tolerances between 0.001 and 0.01--the geometries will snap at 0.003 but not 0.002. At this point I can only assume this is a bug in shapely--any ideas?

I have posted this as a bug report in the shapely repository. It appears to be a bug.


texts: Collect Values not working as expected in ArcGIS ModelBuilder

texts: I am trying to iterate through a folder of rasters, reclassify each raster, then combine all of the reclassified rasters using Cell Statistics (to get the overall SUM of the reclassified rasters).  I have read countless help forums and my Model should be working, but I am not able to connect the "Cell Statistics" block to the "Output Values" from the "Collect Values."  For some reason, the only variables available to be under the "Cell Statistics" input picker box are "Raster" and "Reclassified Rasters" - I would expect "Output Values" to be listed as well, as that is the link I am trying to make.


texts: Will WMTS clients accept redirect (code 307) responses for getTile operations?

texts:
We have all of our tiles available via AWS cloudfront. Cloudfront provides local caches of data, to speed up requests from geographically dispersed clients.

WMTS gettile operations typically return the tile (in say, a jpeg or png) straight to the client. However, instead of piping (via node) the tile straight back to the client, we want to redirect them to the cloudfront URL for that tile (which we has to sign, to confirm they can legitimately access the tile).

Does anyone know if WMTS clients will access this sort of redirect?


texts:
Has a problem. I have an tiff image (that have 4 layers). My tasks is to make small changes in pixel's color to make image better. In this case I use GDAL library. My source is:

GDALDataset  *poDataset;
GDALAllRegister();
poDataset = (GDALDataset *) GDALOpen(fileName.toStdString().c_str(), GA_ReadOnly);
if (poDataset == NULL) {
    QMessageBox::information(0, "error", "We have problems");
} else {
    QMessageBox::information(0, "Message", "All is ok");
}
int rasterCount = poDataset->GetRasterCount(); // Here is 4 raster images
GDALRasterBand *band = poDataset->GetRasterBand(1);
int width = band->GetXSize();
int height = band->GetYSize();

for (int i = 0; i < width; i++) {
    for (int j = 0; j < height; j++) {
        // cross all pixels
        // How to get pixel color here?
    }
}

So I dont know how to get pixel color in cycle. Can you give me advice pleasae?


texts: Moving a layer to a group layer by its index in QGIS 3.0

texts:
I try to use this code from QGIS 2.x version in QGIS 3.0:

li = iface.legendInterface()
li.moveLayer(demLayer, self.watershedGroupIndex

I found QgsLayerTreeLayer but can not use it. I want to move existing layer on layers panel by using gruop index, because I have layer and group index. However examples show me, adding map layer like that:

QgsProject.instance().addMapLayer(demLayer, False)
root = QgsProject.instance().layerTreeRoot()
group = root.findGroup(index.name())
group.insertChildNode(0, QgsLayerTreeLayer(myClone))

But I already loaded layer just want to move it to group by using its index.


texts: How do I set nodata with custom color for gdal_translate and gdal2tilespy?

texts:
I start with a bunch of 1-bit geotiff files that I had to convert to 8-bit RGB tiff files to use with GDAL. Then I tried to set borders to red (content is in black and white so I can't use those) and do

gdal_translate -a_nodata '200 0 0' -of vrt red.tiff nodata.vrt

gdalbuildvrt -srcnodata '200 0 0' -vrtnodata '200 0 0' -hidenodata -addalpha nodata.vrt trans.vrt

gdal2tiles.py -a "200,0,0" trans.vrt ...

But I am only able to hide nodata if it is in grayscale(ex. rgb(200,200,200)) I think the color information is just lost in the gdal_translate process. Is there any way to verify what happened and make gdal work with nodata that has custom colors?


texts: How to embed an OpenLayers map in a blog post?

texts:
I have a map in GeoServer which I need to publish in a blog post.

I'm pretty sure this is possible but I am not managing to make it work.

I'm posting the html source code from the openlayers layer preview I get from geoserver, but not getting the map...

Any help?

Thanks


texts: Using FME to compare multiple polygons (AreaOnAreaOverlayer)

texts:
I would like to use FME to compare multiple polygons from a geodatabase. Namely, I would like to know where the polygons overlap, export this overlap area as a shapefile, and have the attributes of the new shapefile contain all information from both/all polygons which overlap.
I plan to use this workbench frequently in the future, so that each time I receive a new polygon, I can run it through the workbench and see which existing polygons it overlaps with.
I have tried using the AreaOnAreaOverlayer with a Tester set to _overlap>1. When I run this, it does shows me the correct shape polygon where the overlap is. However, I can't figure out how to get the attribute table to say which polygons it is that are actually overlapping. For example, if I read in "Red" and "Blue", the output will say "Red" with two overlaps, but won't mention the name of the other polygon ("Blue").
Just wondering what I am doing wrong - I am familiar with ESRI products but have little experience using FME - I am essentially trying to recreate the Union tool.
Thanks in advance.


texts: Multiple groups for field calculator aggregates?

texts:
Does anyone know if it possible to use multiple groups in the QGIS field calculator?

I have a polygon layer that I want to check that the area calculations are being done correctly. The layer is derived from 3 layers intersected. A paddock polygon is intersected with a 2, 3, 5 and 8km buffer derived from points within the paddock. This intersect is then intersected with a land units layer (soil, vegetation etc).

What I would like to do is to sum the land unit area grouped by the land unit type and the buffered distance eg 8km buffer. But I can't quite figure out how to do it.

I have tried
sum("LUAREAKM2", WA_DIST" AND "LAND_UNIT") which returns zero.

I have also tried to create nested statements but they usually end up invalid.

This is the second time in many weeks that I have tried to figure this out, last week I was trying to calculate some statistics based on multiple groups and ended up resorting to exporting to a spreadsheet and using Excel.


texts: Esri Python script to edit multiple metadata elements failure

texts:
import os
import sys
import xml.etree.ElementTree as ET
import arcpy
from arcpy import metadata as md

# Script arguments...
Source_Metadata = arcpy.GetParameter(0)

# Local variables
#    new purpose text
newPurpose = "This is new text for an existing Purpose metadata element."
newCredits = "This is text for a new Credits metadata element."

# function to update purpose and credits in metadata
def update_metadata(root):
    num_elements = 0

    # modify purpose element's text
    # there is only supposed to be one purpose element in metadata
    # replace purpose element text if element exists
    # if element doesn't exist, do nothing
    purposeEls = root.findall(".//idPurp")
    for element in purposeEls:
        if element.text is not None:
            element.text = newPurpose
            num_elements += 1

    # add credits element to dataIdInfo parent, if the parent exists
    # ISO allows many dataIdInfo groups; ArcGIS supports only one, so get the 1st
    # ISO allows many idCredit elements, and many are supported in ArcGIS
    # append a new idCredit element with appropriate text, existing elements remain
    dataIdInfoEls = root.findall("./dataIdInfo[1]")
    for element in dataIdInfoEls:
        if element:
            newCreditEl = ET.SubElement(element, "idCredit")
            newCreditEl.text = newCredits
            num_elements += 1

    return num_elements

# get and save item metadata
for item in Source_Metadata:
    arcpy.AddMessage("Item: {0}".format(item))

    # get the item's metadata xml
    item_md = md.Metadata(item)
    metadata_xml_string = item_md.xml

    # create an ElementTree object from the metadata XML string
    root = ET.fromstring(metadata_xml_string)

    # call the update_metadata function to modify the item's metadata
    changes = update_metadata(root)

    if changes > 0:
        # get modified XML
        updated_xml_string = ET.tostring(root, encoding="utf-8")

        # import result back into metadata
        arcpy.AddMessage("Saving updated metadata with the item...")
        item_md.xml = updated_xml_string
        item_md.save()
    else:
        arcpy.AddMessage("No changes to save")

arcpy.AddMessage('Finished updating metadata for all source metadata items')

Trying to take this Esri-sourced script to update multiple metadata elements. If I run this in a standalone script without editing, it doesn't do anything, no errors. If I take the script and create a script tool in a GDB inside of ArcGIS Pro, it fails when I select the XML file as a parameter. What do I need to do to get this script to work on my metadata XML?
The script (meant for ArcGIS Pro) was sourced from:
https://pro.arcgis.com/en/pro-app/arcpy/metadata/migrating-from-arcmap-to-arcgis-pro.htm#ESRI_SECTION1_098BC26422304FAE96A19FB4B6528E07


texts: qgis plugin problem: Couldn't load plugin due to an error when calling its classFactory() method

texts:
I am not able to install any plugin.
Packages qtwebkit and qt5-qtwebkit are installed on my machine.

Below is an example for the error ocurring:

Couldn't load plugin 'go2streetview' due to an error when calling its classFactory() method
ImportError: cannot import name 'QtWebKit'
Traceback (most recent call last):
  File "/usr/lib64/python3.6/site-packages/qgis/utils.py", line 336, in startPlugin
    plugins[packageName] = package.classFactory(iface)
  File "/home/florian/.local/share/QGIS/QGIS3/profiles/default/python/plugins/go2streetview/__init__.py", line 34, in classFactory
    from .go2streetview import go2streetview
  File "/usr/lib64/python3.6/site-packages/qgis/utils.py", line 664, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/home/florian/.local/share/QGIS/QGIS3/profiles/default/python/plugins/go2streetview/go2streetview.py", line 23, in
    from PyQt5 import Qt, QtCore, QtWidgets, QtGui, QtWebKit, QtWebKitWidgets, QtXml, QtNetwork, uic
ImportError: cannot import name 'QtWebKit'

Python version: 3.6.5 (default, Mar 29 2018, 18:20:46) [GCC 8.0.1 20180317 (Red Hat 8.0.1-0.19)]
QGIS version: 3.0.3-Girona Girona, exported

Python Path:
/usr/share/qgis/python
~/.local/share/QGIS/QGIS3/profiles/default/python
~/.local/share/QGIS/QGIS3/profiles/default/python/plugins
/usr/share/qgis/python/plugins
/usr/lib64/python36.zip
/usr/lib64/python3.6
/usr/lib64/python3.6/lib-dynload
/usr/lib64/python3.6/site-packages
/usr/lib/python3.6/site-packages
~/.local/share/QGIS/QGIS3/profiles/default/python

What might be the problem here?


texts: Transform the Landsat NDVI value range the similar as MODIS13Q1's range?

texts: I am very interested that how can I transform the historical or overlapping  data(between MODIS13Q1 and Landsat NDVI)into the same value range, are there any regression method to achieve this goal?


texts: Coordinate system for accurately calculating areas of polygons that cross UTM Zones?

texts:
I have multiple polygons, all located offshore around the UK. I am trying to calculate the areas of these polygons in square kilometres, without splitting them by UTM zone.

Is it possible to do this accurately? I have been using World Mercator, but there are big distortions at this scale. If the polygons were limited to the east coast, say, I would use UTM31N. However, as they are all over the UK waters, I am unsure what to use.

Is there a coordinate system that covers this area? British National Grid does not extend far enough out to sea to be of use.


texts: How do I choose the geotransformation used at SpatialReference.org?

texts:
Each page at SpatialReference.org has a UI that updates coordinates as you pan/zoom around a small map control.  Behind the scenes, it is performing a geotransformation.  I looked, but don't see a way to choose which geotransformation to use.

Is there a way to choose the geotransformation, or at least display which one is being used?

I'd like to be able to choose from a pick-list similar to the way I can in Arcmap's geotransformation dialog.

.


texts: can I extract overlap of polygons within the same layer?

texts:
I am trying to extract the overlap of polygons that are within the same layer in 10. The layer came from a buffered street centerline and I'm trying to extract and delete the intersections from the layer.
I have tried 'intersect' but it's virtually impossible to tell which records are an overlap and which ones are not. We are trying Shape tools but so far the computer gave an 'overflow' error.


texts: Add all features to a popup using in Leaflet using WMS GetFeatureInfo

texts:
I'm using the function below to successfully add data to a popup in Leaflet.

Unfortunately, only the data from the last feature found is added. Is there any way to add the data from all features found to the popup.

function identify (e) {

var BBOX = map.getBounds().toBBoxString();
var WIDTH = map.getSize().x;
var HEIGHT = map.getSize().y;
var X = map.layerPointToContainerPoint(e.layerPoint).x;
var Y = map.layerPointToContainerPoint(e.layerPoint).y;
var url = 'http://maps.dartmoor.gov.uk/geoserver/wms/?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetFeatureInfo&layers=general:dnpa_planapps&BBOX='+BBOX+'&FEATURE_COUNT=5&info_format=application/json&propertyName=appcode,url&HEIGHT='+HEIGHT+'&WIDTH='+WIDTH+'&query_layers=general:dnpa_planapps&SRS=EPSG:4326&buffer=15&X='+X+'&Y='+Y;

$.ajax({
    url:url,
    datatype: "json",
    type: "GET",
    success:handleJson
    });
    function handleJson(data) {
    for (var i = 0; i < data.features.length; i++) {
    var feature = data.features[i];
    console.log(feature)
        L.popup()
        .setLatLng(e.latlng)
        .setContent(L.Util.template("
{appcode}

{url}
", feature.properties))
        .openOn(map);
        }
    }
}


texts: Classifying time periods in new field using Python Parser of ArcGIS Desktop?

texts:
I have a datetime field (Horario) in my feature class, showing only the time. Although the date doesn't appear, it seem to be stored automatically as 30 Dec 1899, I'm not sure why.

I want to use this field to populate another field (Periodo) , where I'll have different periods of the day: morning (6 am - 12pm), afternoon (12pm - 6pm), evening (6pm - 12am) and night (12am - 6am).

I'm trying to write an if-elif statement on the Field Calculator, but I'm always having an "invalid syntax" error, and I'm guessing that's because of the datetime expression I'm using. Can anyone see what's wrong?

def Reclass(horario):
    if (horario < date '1899-12-30 12:00:00' and horario >= date '1899-12-30 06:00:00'):
        return "Manhã (6-12h)"
    elif (horario < date '1899-12-30 18:00:00' and horario >= date '1899-12-30 12:00:00'):
        return "Tarde (12-18h)"
    elif (horario <= date '1899-12-30 23:59:00' and horario >= date '1899-12-30 18:00:00'):
        return "Noite (18-24h)"
    elif (horario < date '1899-12-30 06:00:00' and horario >= date '1899-12-30 00:00:00'):
        return "Madrugada (0-6h)"

Reclass(!Horário!)

I'm thinking the problem is I'm referring to the date as

date '1899-12-30 06:00:00'

How should I refer to datetime fields when using Python?


texts: How to select a range of rows with Select by Expression?

texts:
I can't find any similar questions as this. Suppose I have a layer with 20 features, and I only want the first ten.
so something like "col1" = 1:10 and then columns 1:10 are selected only.
How is this done?


texts: How to create a regular polygon grid in Openlayers?

texts:
I am trying to create a grid in my map and do some calculations.

I seen the following links:-

Create regular polygon grid using postgis

Create regular polygon grid inside a polygon using postgis

Here I understood, how to create a polygon grid using postgis spatial query.

Any one, please help me to do the same concept in Openlayers(any client side code).

My arguments are :

Bounds value. Example: new OpenLayers.Bounds(77.06173,28.50872,77.40918,28.74567)
Distance between the points in the grid. Example : 100

My actual requirements is

Create a grid on map.
Plot some POINTS.
Calculate the number of points in each grid.
Get the centroid(x,y) of the grid which have maximum number of points.

I just seen this code..

var points = new OpenLayers.Layer.PointGrid({
    isBaseLayer: true, dx: 15, dy: 15
});

var map = new OpenLayers.Map({
    div: "map",
    layers: [points],
    center: new OpenLayers.LonLat(0, 0),
    zoom: 2
});

Reference: Openlayers Point Grid

My Queries

How can I add polygon grid layer.
How can I set this layer in a particular bounds/extent.
How can I give "select" control to this layer.
How can I find out Number of points on a selected grid/polygon.

Note: I am using Openlayers 2.11.


texts:
I just now edited my Grass7Utils.py file. That didn't work. I have tried changing the path in my Processing toolbox.

QGIS 2.14.3 Mac OSX 10.11.5 I did not mess with my SIP despite running Capitan. This website claimed it was fixed https://github.com/OSGeo/homebrew-osgeo4mac/issues/118#issuecomment-200278686 (I got here from the website in #2 below)

What I have also done (sadly, I can't remember all of it):

started here- QGIS GRASS Missing Dependancy
compiled for macosx using homebrew (I got an error with the install --HEAD grass-71 command (no formulae found in taps). https://grasswiki.osgeo.org/wiki/Compiling_on_MacOSX_using_homebrew
I have tried leaving the path empty in the processing window. I also set the path as: /usr/local/Cellar/grass-70/7.0.4 in my processing window.
I installed the app from GRASS GIS for Mac website and then set my Processing path to: /Applications/GRASS-7.0.app/Contents/MacOS.

I am pulling out my hairs at this point. Almost ready to give up and just install on a VM. I'd much rather do this on the Mac though.


texts: I am trying to make a column in a table with data type geoJSON, I have been told that it is the same as the geometry data type, but when I insert a geojson in the field it give me an SQL error. I am using phpgeAdmin. I am confused.


texts: Store and observe changes of polygons through time?

texts:
What approach should be used to have the ability to store and observe changes of polygons?

For example if you have an area where many parcels of land have changed owners and dimensions (larger/smaller from acquisitions of nearby parcels or subdivisions) through the decades and you want to check the history of a particular polygon?

I was thinking about using "feature compare": http://resources.arcgis.com/en/help/main/10.2/index.html#//001700000004000000

But I am not sure how it works and if it's the ideal tool for this purpose.

Also I don't know how to store these past changes that the parcels would have in order to view how it was 15 years ago for example.


texts: Is there a way in Cartodb.js to include an autcomplete search bar?

texts:
I am using the code below to add a search bar to my mapping application, is there a way to include an autocomplete search bar? Something that would provide suggestions based on what the user is typing?

//add search bar to map
var v = cdb.vis.Overlay.create('search', map.viz, {})
v.show();
$('#map').append(v.render().el);
$('.cartodb-searchbox').find('input.text').attr('placeholder', 'Enter Location...');

I would like to use something that does not point back to my datasets but points to something like google's geocoder, bing, etc.


texts: Difference between FeatureLayer and FeatureSet in ArcGIS

texts: What is the difference between FeatureLayer and FeatureSet in ArcGIS?


texts: I'm new to postGIS, and am looking to use a spatial database with my Rails app (specifically, to retrieve points from the database that are within the bounding box of a Google Maps viewport). I posted a question on Stack Overflow with a few more details ( https://stackoverflow.com/questions/3555840/constucting-a-postgis-query-with-rails ), before coming across this site. I was following the tutorial at http://railsonpostgresql.com/2009/11/05/intro-to-postgis, but my implementation doesn't really seem to work - when zoomed in, the points don't appear as expected (but some appear when zoomed further out, although not those near the international date line). I'm using the spatial_adapter and GeoRuby gems at the moment - are there any good tutorials out there on how to begin integrating postGIS with Rails? From what I've read, it seems like using the geography type would be a good start for my project, but unfortunately I don't believe GeoRuby supports it. Any help/suggestions would be greatly appreciated!


texts: Get extent from layer in ArcGIS MapServer - ArcGIS API 3.28

texts:
I was wondering if it is possible to retrieve the extent of a layer from an ArcGIS map server? I have already tried using the following code:

var url = 'https://fragis.fra.dot.gov/fragis/rest/services/Blocks/MapServer';
var newLayer = new ArcGISDynamicMapServer(url);
newLayer.setVisibleLayers([3]);
newLayer.id = 'a name';
map.setExtent(newLayer.initialExtent);
map.addLayer(newLayer);

But it only returns the absolute max and min values for the extent.

{xmin: -180, ymin: -90, xmax: 180, ymax: 90, spatialReference: {wkid: 4326}}

I've tried using FullExtent and that returns null. I've tried using a FeatureLayer instead of a ArcGISDynamicMapServer with no luck.


texts: RTKGPS+ or RTKLIB configuration of rover and input streams for Stonex S8plus

texts:
I've tried to use RTKGPS+ and RTKLIB with external antenna and via NTRIP to communicate with Continuously Operating Reference Station (CORS) of my country (MAKPOS), but I couldn't make connection (i have user, password, port 9001, NTRIP mountpoint: iMAX-Auto, host: makpos.katastar.gov.mk) I have Stonex S8plus and i like to try to use smartphone+external antena via NTRIP to connect with MAKPOS and do Network RTK or DGPS.

Simply the configuration of this system is by these segments:

Referent base stations (14 stations in radius of 50-70km between each of them). Leica GRX 1200+ & Leica AR 25 … 9 stations Leica GRX 1200 GG Pro & Leica AT 504 GG … 2 stations Leica GX 1230 GG & Leica AX 1202 GG … 3 stations Control Center: Leica GNSS Spider (controlling and managing with network, RT Proxi Server and NTRIP Caster) Leica GNSS Spider QC (controling the quality of data and coordinates on network points) Leica GNSS Spider WEB (distribution of RINEX) Services in MAKPOS DGPS (accuracy 0.3-0.5 m), it uses RTCM 2.x and the transfer is by GPRS and NTRIP RTK (accuracy 0.02-0.04m), it uses RTCM 2.x and RTCM 3.x and the transfer is by GPRS and NTRIP Precise positioning with additional processing, (<0.01 m), RINEX, transfer via Internet, mail..

Does anyone have done RTK or DGPS with this app or with another app? Where can I find configuration file for Stonex S8plus?

I made configuration to sent NMEA via bluetooth in Stonex Assistant, at 1hz page 52 in user manual.

Link of RTKGPS+ https://play.google.com/store/apps/details?id=gpsplus.rtkgps&hl=en

Stonex S8plus user manual: http://www.stonexpositioning.com/jdownloads/Public/Manuals%20and%20tutorials/English/GNSS%20receivers/S8%20PLUS/s8_plus_user_manual.pdf


texts: Adressing files/directories in ogr2ogr using OSGeo4W shell?

texts:
I want to convert a .pbf-file (freiburg-regbez-latest.osm.pbf) into a SpatialLite Data Base using ogr2ogr. The file is located at W:\geodaten\UG
I open OSGeo4W.bat in the directory C:\Program Files\QGIS 2.18..
I enter following code:
C:\>ogr2ogr -f "SQLite" -dsco SPATIALITE=YES basemap_ug.db freiburg-regbez-latest.osm.pbf

However, this does not work; I get:

C:>ogr2ogr -f "SQLite" -dsco SPATIALITE=YES basemap_ug.db freiburg-regbez-latest.osm.pbf
ERROR 4: sqlite3_open(basemap_ug.db) failed: unable to open database file
ERROR 1: SQLite driver failed to create basemap_ug.db

Probably I need to specify the directory of the file. Unfortunately, I don't know how to do this.


texts: Use variables in PostgreSQL/PostGIS code

texts:
I'd like to create an SQL file, where the user can just change variables on the beginning , and use it in different configuration system. Typically for load layer from different PostGIS server with postgres_fdw.

I tried it with PL/pgSQL where I can declare variables before beginning processing, but I don't no how to use them after...

DO $$
DECLARE
    url character := '193.00.00.00';
    pt character := '5432';
    db character := 'bd_name';
    usr character := 'user_name';
    pwd character := 'password';
BEGIN
    CREATE EXTENSION IF NOT EXISTS postgres_fdw;

    CREATE SERVER foreign_srv
        FOREIGN DATA WRAPPER postgres_fdw
        OPTIONS (host url, port pt, dbname db);

    CREATE USER MAPPING FOR USER
        SERVER foreign_srv
        OPTIONS (user usr, password pwd);

END $$;

It obviously doesn't work like that.


texts:
I've been looking to use the Ordanance Survey's Meridian data (part of the their OpenData suite) for doing vehicle routing.

Whilst I can infer road class from the NUMBER property on each polyline, I cannot see anything that indicates if a polyline is unidirectional, or bidrectional. Am I just being unable to see this, or is this encoded somewhere else?


texts: How to connect to the ArcGIS Server through a firewall?

texts:
I know there were a lot of similar questions, but still can't find a correct description about that. What should be inbound and outbound rules on the server side? Which ports should I use? Will that be a local or a remote port in the rule's options? Is that important (or can I) rewrite the service.txt to not to be INCREMENT and to be FEATURE? Please write me a correct answer step-by-step.

I don't understand, because you can see some of the drawings, for example:
hardcoding port number in Arcgis License Manager

says the License Server Manager is with the ArcGIS Server, not with the ArcGIS Desktop. But I can't install the License Manager to be with the ArcGIS Server. So what am I doing wrong? On my PC there's the ArcGIS Desktop and the License Server Manager. The ArcGIS Server is on an other computer. I can connect to that without the server's firewall, but that should be active...


texts: Move Raster image in raster shape fill in QGIS

texts:
I am trying to make a simple flag map of the world. I use raster fill and import an image of the flag, which is a PNG file, and that fills up the polygon of the country.
While I can increase and decrease the size of the image, I am not able to move it.
In this example I want the blue centre wheel of the indian flag to be in the centre of the country.

When I offset the image, the whole country seems to move but the flag position remains the same

When I use draw effect to translate the image even a bit, the whole image disappears
Is there any way to move the raster image?
I use QGIS 2.18


texts: Using Weighted sum and Grouping Analysis?

texts:
I have a collection of points representing buildings, and each has a number value. I need to look at building within 50 miles of each other, where the total number values of the building in that cluster is less than 200.

I'm having trouble figuring out how to put this together.

Is there a way to create clusters based on where the value of points within 50 miles is less than 200?

I feel like a should create a buffer, and than than calculate wieghted sum. What I would really like to do is create clusters (maybe with grouping analysis?), and and then select the point that has a highest wieght from there.

Could someone point me in the right direction getting started?


texts: Visual Studio's CMake with vcpkg: Error gdal is not found

texts:
I have been trying to build a cmake c++ project. More specifically I am trying to use the gdal library in this project. In the CMakeLists.txt it says find_library(GDAL gdal) after doing some research i found, that visual studio can open cmake files by default as mention in this thread: https://docs.microsoft.com/en-us/cpp/build/cmake-projects-in-visual-studio?view=vs-2019.
Moreover, visual studio should also automatically include the gdal library once i have set it up with vcpkg correctly. I've already downloaded the x64-windows version of the library (vcpkg install gdal:x64-windows) in order to build for the right architecture and made it available via vcpkg integrate install on a user-wide scope.
After some trial and error, everything works fine now, the toolchain gets included accordingly and the library is found automatically, resulting in a configuration like that:

However, when trying to include the header files (or anything else; see code snippet), visual studio does not seem to link the library correctly as it will result in the error message: cannot open source file "gdal/gdal.h".
#include <gdal/ogrsf_frmts.h>
#include <gdal/gdal.h>
#include <gdal>

Where should I further investigate?


texts: Given a GeoTIFF that has a map, extra textual descriptions of rock, and cross-sections all in one .tif file, is there an automatic way to crop the GeoTIFF to just the map? A programmatic solution in Python or one that uses GDAL would be preferable.


texts: Labels appear and disappear for no apparent reason using QGIS?

texts:
I'm going through the Alaska sample data tutorial for QGIS 2.2.0. Some of the lake labels appear and disappear at different scales when zooming. At 1:314,460 Becharof is visible, but Lower Ugashik doesn't appear until 1:157,230. It seems like others have had similar problems.

I have played with the rendering and placement dialogs to no avail. Even if I don't have the "magic combination" of parameters set, shouldn't they all behave the same?

For my needs, this may well render (no pun intended) the program useless. I need it to print labels when I need them!


texts: How can I connect seven schools worldwide to create and share QGIS projects?

texts:
I am looking at a possible programme to connect seven schools ( all English capable ) located around the world to set up and share QGIS data and projects about their own local area among the group.

What type of open source software ( OSGEO ???) is available to do this and where can I find documentation on how it all goes together?  Is there a cloud server system that I can use?

Has this been done before using QGIS?


texts: Generating MapQuest Style tiles for offline use

texts:
I need to create about a 9-level deep SLIPPY map of the world from OSM data (for offline use). So I cannot use the MapQuest, CloudMade, etc. servers). The projection must be EPSG4326 (LL84) which is another reason I cannot use them. The viewer is Leaflet which I have working fine with EPSG4326 tiles (currently I am diced the NASA Blue Marble to 9 levels).

I like the look of the mapquest tiles the best: http://otile{s}.mqcdn.com/tiles/1.0.0/osm/{z}/{x}/{y}.jpg

but I could deal with others. The main reason I like MapQuest is the worldwide english labeling (with foreign subtitles) and the look of the topographic terrain (useful in regions with no roads).

Tools that I have at my disposal include FME Professional (my favorite), QGIS, PostGIS 2.0, and GDAL. I am working on a large Windows 7SP1 box.

I cannot use cloud services like like mentioned in the following article at OSM.org:

http://wiki.openstreetmap.org/wiki/Deploying_your_own_Slippy_Map

See Switch2Osm.org for more information on serving your own tiles.
Tiledrawer is a free, quick-start, out-of-the-box renderer and server using the standard OSM toolchain for Amazon EC2, other virtual servers, and fresh installs of Ubuntu on any hardware.


texts: Loop in arcpy.describe with single print if all rasters have same properties?

texts:
I have a bunch of rasters with different file formats in a folder (most of them are GRID).  Using describe of arcpy, I want my script to inform the following: when all rasters have the same format, inform once the file format e.g print "All the rasters has the XXX format", but when there are different formats, list the formats of each raster. I'm very new to scripting so by now the code I wrote print all the files with their format.

How can I do that?

The code block is:

for dpath, dnames, fnames in os.walk(lu_folder):
    arcpy.env.workspace = dpath
    print "\nThe format of LU rasters are:"
    for raster in arcpy.ListRasters('*'):
        desc = arcpy.Describe(raster)
        print  "{0}: {1}".format(raster, desc.format)


texts: Problem in accessing data form Landsat5 in Google Earth Engine

texts:
I have written a code to calculate NDWI from LANDSAT5 and to creat a time series chart.
the problem is that I get very weird results- it says there are only 2 images  between 1990-2011 and for that reason I get wrong graph.
 what can be the reason? is this data is blocked for use or there is something wrong in my code?

this is the code I have written:

/**
 * Function to mask clouds based on the pixel_qa band of Landsat SR data.
 * @param {ee.Image} image Input Landsat SR image
 * @return {ee.Image} Cloudmasked Landsat image
 */
var cloudMaskL457 = function(image) {
  var qa = image.select('pixel_qa');
  // If the cloud bit (5) is set and the cloud confidence (7) is high
  // or the cloud shadow bit is set (3), then it's a bad pixel.
  var cloud = qa.bitwiseAnd(1 << 5)
                  .and(qa.bitwiseAnd(1 << 7))
                  .or(qa.bitwiseAnd(1 << 3));
  // Remove edge pixels that don't occur in all bands
  var mask2 = image.mask().reduce(ee.Reducer.min());
  return image.updateMask(cloud.not()).updateMask(mask2);
};

var dataset = ee.ImageCollection('LANDSAT/LT05/C01/T1_SR')
                  .filterDate('1990-01-01', '2011-12-31')
                  .filterBounds(geometry)
                  .map(cloudMaskL457);

var visParams = {
  bands: ['B3', 'B2', 'B1'],
  min: 0,
  max: 3000,
  gamma: 1.4,
};

//clip the dataset according to the geometry
var clippedCol=dataset.map(function(im){
  return im.clip(geometry);
});

// Get the number of images.
var count = clippedCol.size();
print('Count: ',count);
//print(clippedCol);
print(dataset,'dataset');

//function to calculate NDWI in LANDSAT5
var addNDWI = function(image) {
  var NDWI = image.normalizedDifference(['B4', 'B5'])
  .rename('NDWI')
  .copyProperties(image,['system:time_start']);
  return image.addBands(NDWI);

};

//NDWI to the clipped image collection
var withNDWI = clippedCol.map(addNDWI).select('NDWI');

var NDWIcolor = {
  min: -1,
  max:1,
  palette: ['0000ff', '00ffff', 'ffff00', 'ff0000', 'ffffff'],
};

var listOfImages =(withNDWI.toList(withNDWI.size()));
var randomImage=ee.Image(listOfImages.get(0));
Map.centerObject(geometry,12);
Map.addLayer(randomImage,NDWIcolor,'Random_Image_test');

print(ui.Chart.image.series(withNDWI, geometry, ee.Reducer.mean(), 30));


texts: Can't display vector layer using L.vectorGrid.protobuf

texts:
I use mapbox/tippecanoe  to create .pbf tiles.
But when I try to display them using L.vectorGrid.protobuf nothing is rendered, although I see that proper .pbf files are requested and downloaded by browser.

I suspect that missing options parameter is the root of the problem but I am not sure:

var layer = L.vectorGrid.protobuf('http://my_static_files_server/{z}/{x}/{y}.pbf')
map.addLayer(layer)

Is it necessary to set options if I want just display vector layer represented by single polygon?

How can I set common options?


texts: Running ZTV (ZVI) on ArcGIS Desktop for 3m high object?

texts:
What tool in the spatial toolbox would you use to run a ZTV (ZVI)? or is there a better way of carrying this type of query out using ArcGIS 10?

Also what other information would I need to do this type of query?


texts: Visualising opening hours of a building with rule-based symbology [QGIS]

texts:
I have the shapefile of a neighbourhood of buildings. Within the attribute table of this shapefile I have the opening and closing hours of each building, written in decimal format (e.g. 10:30am is 10.5, 5:00pm is 17). Opening and closing times are in their own real number fields each.

How can I create a rule-based symbology that shows those buildings that are open, from say 7am-9pm, if I have only the opening and closing time? How can I ask QGIS to compute the complete range of hours a given building is open for, and based on the values returned from that, return either a positive or negative for the purposes of symbology?


texts: Loading JSON with image link and extents into ArcGIS Pro

texts:
One of the ways to get data from ArcGIS REST API in image format is this.
https://services.twdb.texas.gov/arcgis/rest/services/Base/BaseLayerQueryService/MapServer/export?bbox=-11360770.018835%2C3942451.651995%2C-11357251.989618%2C3945703.856100&bboxSR=&layers=show%3A8&layerDefs=&size=847%2C783&imageSR=&format=png&transparent=true&dpi=&time=&layerTimeOptions=&dynamicLayers=&gdbVersion=&mapScale=&rotation=&datumTransformations=&layerParameterValues=&mapRangeValues=&layerRangeValues=&f=json

The return is a JSON like this
{
 "href": "https://services.twdb.texas.gov/arcgis/rest/directories/arcgisoutput/Base/BaseLayerQueryService_MapServer/_ags_mape38a2c53b769401a86d7d8f1751fc8a2.png",
 "width": 847,
 "height": 783,
 "extent": {
  "xmin": -11360770.018835019,
  "ymin": 3942451.651995,
  "xmax": -11357251.989617988,
  "ymax": 3945703.8561,
  "spatialReference": {
   "wkid": 102100,
   "latestWkid": 3857
  }
 },
 "scale": 15698.301990971015
}

If I save this as a file, is there a way to load this in ArcGIS Pro?


texts: Leaflet: how to use a custom marker on a geojson layer?

texts:
I'm trying to put on a Leaflet map a geojson and all works fine until I use the default blu marker.

Now I'd like to use a custom marker (a little .png icon) and I've changed my code in the follow

 var my_json;
 $.getJSON('../Dati/my-geojson.geojson', function(data) {
           my_json = L.geoJson(data, {
            pointToLayer: function(feature, latlng) {
                var smallIcon = L.Icon({
                    options: {
                        iconSize: [27, 27],
                        iconAnchor: [13, 27],
                        popupAnchor:  [1, -24],
                        iconUrl: 'icone/chapel-2.png'
                    }
                });
                return L.marker(latlng, {icon: smallIcon});
            },
           onEachFeature: function (feature, layer) {
                   layer.bindPopup(feature.properties.ATT1 + '<br />'
                                                 + feature.properties.ATT2);
           }
         });
 my_json.addTo(markers.addTo(map));
 TOC.addOverlay(my_json, "My layer name in TOC");
 map.removeLayer(my_json);
 });

the error that I can see in Firebug is

TypeError: this.options.icon is undefined
var anchor = L.point(this.options.icon.options.popupAnchor || [0, 0]);

something is going wrong but I don't know how to fix it.


texts:
Are there any known bugs with the getCover argument in the rasterize function?
(R package 'raster'). or has anyone got a fix for the following:

I have a series of landcover polygons, here is a small sample of them in shp format, that i want to rasterize based on coverage at 25m resolution. There are 5 types of landcover.

I would like to perform a rasterization that shows the percentage of each cell that is covered, per Land Cover type (i get 5 resulting rasters that I can stack), using the getCover option in 'rasterize' in R.

Here is some code, im very sorry it is only a download dataset;

require(raster)
require(rgdal)

setwd("C:\\your_wd")

# read in shapefiles, extract .prj
fields <- readOGR(".", "sample")
projGB <- projection(sample)

# make a 25m raster (made specifically for above data)
rOS <- raster(nrows=16, ncols=27, xmn=254110.156238, xmx=254785.156238, ymn=251530.203086, ymx=251930.203086, crs=projGB)

# rasterize the polygons per LandCover code in a loop, using getCover

cats <- unique(fields$LAND_CAT)
st <- stack()

for (i in cats){
 x <- fields[fields$LAND_CAT == i,]
 y <- rasterize(x, rOS,getCover=TRUE)
 writeRaster(y, filename=paste("fields_",i,"_25m"),  format="GTiff", overwrite=TRUE)
 st <- stack(st, y)
}

Unfortunately, the produced GeoTiffs do not show all coverage stats across the polygons. Only some polygons appear to have been used, with no detectable pattern. The biggest category 'GP' is the most problematic.

I have:

plotted all polygons as a group, per Land Cover and individually to check they all plot (they do)
used fortify to check for holes etc, there are none

require(ggplot2)
fields$group <- 1
fields.f <- fortify(fields,region="group")
rasterized the whole dataset disregarding landcover type and not using getCover, just to check - all is fine.
plotted individual polygons and rasterized each one using getCover - no problem
used 2 other small samples from the big dataset but the problem persists
used another shapefile from another project, problem persists

I wonder why getCover in rasterize does not work and if anyone has had this problem?

R version 3.0.3

If and when I get this fixed, I will use the raster stack to then determine which land category has the biggest value in every cell to make a raster that is based on majority cover instead of centroid. Doing this in R seems long winded to me.


texts: Trying to get Spiderfier to work with MarkerClusterer in Google Maps

texts:
I have just over 100 points to display on a map and had it working fine.  Quite a few of these points are displayed at the exact same coordinate position so sit on top of each other. To make it so these are able to be viewed I added Spiderfier which also worked fine. However I wanted to have it so that these were represented with a number so I tried to add MarkerClusterer. Now while MarkerClusterer is working it is overwriting spiderfier and now they won't split when clicked.

How can I get them to work together?

Here is my code:

<script type="text/javascript">

  google.load("visualization", "1", {packages:["corechart"]});

<!--Make Map-->
window.addEventListener('load',function initialize() {
    <!--Map Options-->
    var mapDiv = document.getElementById('map-canvas');
    var MapOptions = {
    center: new google.maps.LatLng(-31.950804, 115.860485),
    zoom:10
    };

    <!--Load Map Window-->
    var map = new google.maps.Map(mapDiv,MapOptions);

    var oms = new OverlappingMarkerSpiderfier(map);

    <!--Query Fusion Table "Teams"-->
    var query = "SELECT 'TEAM NAME',SUBURB,Lat,Lng,LogoName FROM 15XW2l-mxElk3qejK2lMgV2kukDh-FMOdjl6IRnwO";
    query = encodeURIComponent(query);
    var gvizQuery = new google.visualization.Query(
        'http://www.google.com/fusiontables/gvizdata?tq=' + query);

    <!--Create Marker Function-->
    var createMarker = function(teamname,suburb,coordinates,icon) {

        var marker = new google.maps.Marker({
            <!--Marker Options-->
            map: map,
            position: coordinates,
            icon: icon,
            title: teamname

      });

       oms.addMarker(marker);

      <!--Listener for clicking point "logo" -->"
      google.maps.event.addListener(marker, 'click', function(event) {
        infoWindow.setPosition(coordinate);
        infoWindow.setContent(teamname);
        infoWindow.open(map);
      });

      return marker

    };

    <!--get number of rows-->
    gvizQuery.send(function(response) {
        var numRows = response.getDataTable().getNumberOfRows();

        var markers = [];

        <!--Loop through rows-->
        for (var i = 0; i < numRows; i++) {

            <!--Get row values-->
            var teamname = response.getDataTable().getValue(i, 0);
            var suburb = response.getDataTable().getValue(i, 1);
            var lat = response.getDataTable().getValue(i, 2);
            var lng = response.getDataTable().getValue(i, 3);
            var logoname = response.getDataTable().getValue(i, 4);

            var coordinates = new google.maps.LatLng(lat, lng);
            var icon;
            icon = new google.maps.MarkerImage(
                'images/'+logoname+'.png',
                null,
                null,
                new google.maps.Point(0, 0),
                new google.maps.Size(40, 40)
                );
            <!--Feed row values to createMarker-->
            var mark = new createMarker(teamname,suburb,coordinates,icon);

            markers.push(mark);

        };

    var markerCluster = new MarkerClusterer(map, markers,
        {imagePath: 'https://developers.google.com/maps/documentation/javascript/examples/markerclusterer/m'});
    });

});

</script>


texts:
How can I integrate RGB image and nDSM data provided in the dataset (ISPRS Vaihingen) that I want to use in a multimodal deep learning based semantic segmentation method? Can anybody provide a Python code to do so.
The below code provides a way to make semantic segmentation only from RGB images.
https://github.com/nshaud/DeepNetsForEO/blob/master/SegNet_PyTorch_v2.ipynb
But how the data preprocessing has to be done to make the fusion of RGB and nDSM data both to make the segmentation?


texts: How to close the current mxd document through arcmap add-in

texts:
First, I open a mxd document named A in arcmap, then i use my add-in to generate another  mxd document named B,

then I use

 ArcMap.Application.OpenDocument(arcmapMapFileName);

to open the document, but the document does not display normally, it has both A and B's layers.
I want to only open B, so before open B, I have to close A.

How can I achieve this through an arcmap add-in in c#?


texts: Influence of the scale factor on the projection

texts:
I'm currently reading documentation about map projections to understand the source code of the Proj4 project.

The scale factor is named in a variety of sources I read. This sources explained its definition and its value for some projections.

In the source code of Proj4, for the mercator projection (sphere and ellipse case), the scale factor influences the coordinates on the projection :

//P->k0 is the scale factor
xy.x = P->k0 * lp.lam;
xy.y = - P->k0 * log(pj_tsfn(lp.phi, sin(lp.phi), P->e));

Why and how I should use the scale factor during the computation of the projection ? Is there any valuable resources on the web ?

This question in asked in the sense of projection computation. I can find the formula for the inverse and forward projection as well as the scale factor in a various of resources, but no one explains how I should use both in a algorithm. You have the definition of the projection and the definition of the scale factor, but it's not clearly written that I should multiply or divide the result by the scale factor.

Is it a general rule : if I find the formula for any projection with the related scale factor, should I, in all cases, always divide or multiply the results by the scale factor ?


texts: Using ArcGIS Python API to consume standalone ArcGIS REST API

texts:
I have access to a ArcGIS REST API service running on a standalone, unfederated ArcGIS Server. (Not AGOL or ArcGIS Enteprise) ie  (https://xxxx.xxxxx.com/arcgis/rest/services/xxxxxx/xxxxMap/MapServer)
Is there any way to consume the REST API services like MapServers, ImageServers using ArcGIS's Python API?
Also, I am working in a secured environment with no internet access. Using MapImageLayer from ArcGIS's Python API forces a GIS object as a parameter by default. And that GIS object calls out to www.arcgis.com for anonymous authentication.


texts: QGIS an PDF templates: getComposerItemById is changed with the QGis version?

texts:
In a QGis python plugin i have this code to fill in the fields of a PDF template (.qpt):

    myComposition = QgsComposition(myMapRenderer)
    myFile = path + "template.qpt"
    myTemplateFile = file(myFile, 'rt')
    myTemplateContent = myTemplateFile.read()
    myTemplateFile.close()
    myDocument = QDomDocument()
    myDocument.setContent(myTemplateContent)
    myComposition.loadFromTemplate(myDocument)

    item = myComposition.getComposerItemById("id_field")
    item.setText("my value")

working with QGis 2.18.7 item is a qgis._core.QgsComposerLabel and the method setText() is correct.

working with QGis 2.18.13 item is a qgis._core.QgsComposerItem and the method setText() not exist.

Is getComposerItemById changed with the QGIS version?


texts: ArcGIS: ModelBuilder output is not added to the display

texts:
At 10.7.1, I have a simple model that is supposed to add the final FGDB feature class to the ArcMap display:

I also have the correct Geoprocessing option set:

I run the tool in ArcMap by double-clicking it from its toolbox, not from within the editor window, as explained in this SE GIS post.

Sadly, the feature class is not added to the display.  However, the model really does create the feature class, and I can manually add it to the display.

What am I missing?


texts:
Wondering if there is a way to read and access realtime JSON data within QGIS. I am aware of how to work with realtime GeoJSON, however I would like to work with JSON data, which does have lat/long coordinates, like in the link provided below. The ultimate goal is to access realtime data within QGIS. Is there any such tool like what is available to work with CSV data to convert to geo, but that will allow for realtime updates?

https://tor.publicbikesystem.net/ube/gbfs/v1/en/station_information


texts: CRS is UTM 51, the coordinates are shown in meters, show same coordinates as lat long as well?

texts:
We are working with a CRS of UTM 51N, the Philippines. The coordinates under cursor are shown as metric values. Clients however like so see the Latitude Longitude values when they move the cursor.

Is there a way to display the cursor position in Latitude Longitude as well?

Thank you Frederik


texts: Create point between two points with given distance or length

texts:
What I have:

Coordinates of Point A and Point C and therefore vector and azimuth from A to C

What I want:

Coordinates of Point B (between A and C) with a distance from A to B of 35.2 meters.
Suppose to work in every coordinate system

Some code:

import shapely.geometry as g
line = g.LineString([(0, 0), (0, 100)])
print line.interpolate(35.2).wkt

..works in UTM 32N but not in geographic coordinates?!


texts: Selected output format is not supported with this function? Wikimapia API

texts:
I have been trying with limited luck to instigate the &format=kml with the place.getbyarea function. I have been trying to isolate any additional paramater that might be conflicting with the exporting process. This has produced the same error every time.

http://api.wikimapia.org/?key=PUT_YOUR_KEY&function=place.getbyarea&coordsby=latlon&lon_min=35.37735&lat_min=33.055148&lon_max=36.333847&lat_max=34.694027&format=&pack=&language=en&page=1&count=50

Is there an issue with any of this input?


texts: QGIS Retrieve field used in data-defined override

texts:
I have a point layer in QGIS with an SVG marker. I set the rotation via a data-defined override to the value of one field in the layer.

How do I retrieve that field in Python? Not the numeric value for a given feature, but the field (ie the expression) used?


texts: Configuring lidR catalog to save files with same filename?

texts:
I'm trying to learn how to use the options of catalog() within the lidR package in R. I'd like to save directly the processed files, for example, using grid_terrain() function over a catalog and save files keeping the filename of the original LAZ/LAS file.

As you can see in the guide of the package, catalog has the option of save the file using things like {XBOTTOM}_{ID}:

# Internal engine will not return results into R. Instead it will write results in files.
opt_output_files(ctg) <- "/path/to/folder/templated_filename_{XBOTTOM}_{ID}

I'd like to save the files using the same filename, but, I don't know how to configure that part with {} in the opt_output_files() option. I tried several things like, for instance: opt_output_files(cat) <- paste0(output,"/{data$filename}") but, it doesn't work.

lasdir <- "C:/lazfiles"
output <- "C:/output"

cat <- catalog(lasdir)
lasfiles <- cat@data$filename #with this you can see the filenames
opt_progress(cat) <- TRUE
opt_output_files(cat) <- paste0(output,"/{data$filename}")
opt_cores(cat) <- 3
opt_chunk_buffer(cat) <- 20

#function that I want to use over the catalog files
mdt <- grid_terrain(cat, res = 5, algorithm = knnidw(k = 5, p = 2))


texts: Can "Getting to know ArcGIS Desktop" for version 9 be used as a tutorial for version 10?

texts:
I am a GIS basic/novice user.

Our company has ArcGIS 10 desktop, and since it is multi-user license, I was allowed to install on my computer. I was given a copy of "Getting to know ArcGIS desktop" but it is for ArcGIS 9.

Are the two versions so different that I should just get myself the latest copy for 10? ESRI sells it for about 70 dollars so I'd like to not spend the money, but I'd also like to learn the system as well as possible.


texts: PostGIS ST_Collect/ST_Dump not working as expected?

texts:
I wanted to calculate the number of points that fall within some distance of a point inside a PostgreSQL/PostGIS window function, and I ran into a problem I don't understand.  The core of the problem is as follows:

This simple query finds the number of points within 2500 m of the point with ID 863040 (the point geometries are in UTM Zone 48, SRID 32648):

WITH onepoint AS (SELECT id, geom FROM mytable WHERE id = 863040),
     allpoints AS (SELECT geom AS points FROM mytable)
SELECT op.id,
 sum(ST_DWithin(op.geom, allp.points,2500)::int)
FROM onepoint op, allpoints allp

but if I try the same query, this time collecting the geometries of "allpoints" and then re-dumping them, the query fails:

WITH onepoint AS (SELECT id, geom FROM mytable WHERE id = 863040),
     allpoints AS (SELECT ST_Collect(geom) AS points FROM mytable)

SELECT op.id,
 sum(ST_DWithin(op.geom, (ST_Dump(allp.points)).geom,2500)::int)
FROM onepoint op, allpoints allp

> ERROR: set-valued function called in context that cannot accept a set

I needed to collect and dump to be able to use ST_DWithin inside the window function, but my plans have all come to naught because of this error.


texts: GDAL_TRANSLATE Assumes Picture is Orientated North?

texts:
I used gdal_translate to change a jpeg to a geotiff. Attempting to use gdalwarp to create a mosaic I found the images skewed and this lead me to the observation in the title of the question. Below I've posted the line I used to translate and gdalinfo on the new image. Notice that the upper left and lower left longitudes are the same (-76.4303686) and the upper right and upper left latitudes are the same (38.1480727).

Is this a flaw/bug or am I doing something wrong?

C:\Users\path\to'Desktop>gdal_translate -a_ullr -76.43036800559511 38.14807271147
146 -76.43183999600483 38.148231283935985 IMG028.jpg img28.tif
Input file size is 400, 267
0...10...20...30...40...50...60...70...80...90...100 - done.

C:\Users\path\to\Desktop>gdalinfo img28.tif
Driver: GTiff/GeoTIFF
Files: img28.tif
Size is 400, 267
Coordinate System is `'
Origin = (-76.430368005595113,38.148072711471457)
Pixel Size = (-0.000003679976024,0.000000593904362)
Metadata:
  EXIF_BitsPerSample=8 8 8
  EXIF_ColorSpace=65535
  EXIF_DateTime=2013:11:02 15:13:42
  EXIF_ExifVersion=0221
  EXIF_Orientation=1
  EXIF_PhotometricInterpretation=2
  EXIF_PixelXDimension=400
  EXIF_PixelYDimension=267
  EXIF_ResolutionUnit=2
  EXIF_SamplesPerPixel=3
  EXIF_Software=Adobe Photoshop CC (Windows)
  EXIF_XResolution=(72)
  EXIF_YResolution=(72)
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( -76.4303680,  38.1480727)
Lower Left  ( -76.4303680,  38.1482313)
Upper Right ( -76.4318400,  38.1480727)
Lower Right ( -76.4318400,  38.1482313)
Center      ( -76.4311040,  38.1481520)
Band 1 Block=400x6 Type=Byte, ColorInterp=Red
Band 2 Block=400x6 Type=Byte, ColorInterp=Green
Band 3 Block=400x6 Type=Byte, ColorInterp=Blue

I'm new to GIS and really programming in general.

The pictures were taken from a high altitude and the camera was pointed directly below. The picture is effectively a rectangle.
I do not have 'world files.' That being said, I do have information that was read from a GPS unit. The information is: latitude, longitude, altitude, heading (degrees clockwise from North). These co-ordinates are assumed to be at the exact center of the picture. The hard-coded ulx uly, lrx, lry values are a result of trigonometry. They were intended to represent the latitude and longitude of the upper-left and lower-right corners of the picture (which is what I understood GDAL wanted). Thus it might totally be possible that the upper-left corner of the picture is actually more southward than the lower right (i.e. the picture is upside down). I thought GDAL automatically reads the size of the image and assigns proper bounds for the upper right and lower left corners.

I read over the translate docs again and there seem to be no parameters with which I can tell GDAL to read data about orientation using a 'world file.'

What methods are available to factor in the rotation of the image?

Can someone provide easy-to-read links on how to use gdal_translate properly?

In particular, multiple methods of factoring in image rotation would be preferred. I do not want to create world files. It would be good if I could somehow factor in rotation using command line flags or something. Additionally, I'm unclear about what a projection is. I know 'projection' is like assigning parameters for what a pixel means in meters and what latitude and longitude means. But something like this is not helpful:

-a_srs srs_def:
Override the projection for the output file. The srs_def may be any of the usual GDAL/OGR forms, complete WKT, PROJ.4, EPSG:n or a file containing the WKT.

I will continue to work on projection because I have a few leads but this 'world file' business and image rotation I do not even know where to start.


texts:
I'm using gdal_rasterize function to create a raster from my PostGIS database which contains measurements over the entire earth of a gas concentration. (almost 600'000 points over the polygon below)
my goal is to create a raster with resolution 700x700 over a rectangular polygon.

'POLYGON ((-11.0000000000000000 34.0000000000000000, 42.0000000000000000 34.0000000000000000, 42.0000000000000000 72.0000000000000000, -11.0000000000000000 72.0000000000000000, -11.0000000000000000 34.0000000000000000))'),4326))

my problem is that setting:

-ts 700.0 700.0

-te 34 -11 72 42

does not return a correct raster over the polygon but a smaller square with
latitude in (34,42)  and longitude in (33,42)

I have tried to set the normal -te -180 -90 180 90 but then i get the correct
polygon but not with the resolution i want (since the whole world is 700x700).

this is the full command i'm using

gdal_rasterize -l public.measurements -a no2_tropo -ts 700.0 700.0 -a_nodata 0.0 -te -180.0 -90.0 180.0 90.0 -ot Float32 -of GTiff "PG:dbname='testdb' host=localhost port=5432 user='postgres' password='admin' sslmode=disable" /home/dave/Desktop/fullPoints4k.tif^C


texts: KML - Network Link Control 'Update' Usage

texts:
I am having trouble getting a NetworkLinkControl KML command to work.

The drive.google.com URL matches Network Link to Slot 2 and I am trying to use NetworkLinkControl to delete kmlslot2.

I know the targetHref tag works, because I have successfully been able to append a placemark under the Network Link to Slot 2 directory, but have been unable to get the delete function to work.

Below is the KML text, along with an included picture of the KML that I am trying to modify.  I am using Google Earth Pro 7.1.8 on macOS for reference.

**I just realized I was confusing name and object id, but even after adding an object id ex: <NetworkLink id="kmlslot2"> it still isn't working.

Any idea why this isn't working?

<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://earth.google.com/kml/2.2">
  <NetworkLinkControl>
    <Update>
      <targetHref>https://drive.google.com/uc?export=download&amp;id=**</targetHref>
      <Delete>
        <NetworkLink targetId="kmlslot2"/>
      </Delete>
    </Update>
  </NetworkLinkControl>
</kml>


texts:
How to load custom style (mapbox GL JSON) when running a docker container of tileserver-gl?
I have my own MBTILES vector tiles made with tippecanoe and Ive made a stylesheet.
Want to load that JSON file to use with the mbtiles.
I need to rasterize (request PNG) these tiles for a mobile client that does not support vector tiles.


texts: Installing Lecos plugin for QGis 2.0.1 on linux fedora 20

texts:
I try to install Lecos plugin, which requires python imaging library (pil), but the installation fails with the following: "you need to have the image library pil installed"

I installed the python-pillow library instead of the pil library: pil is almost dead and I can't find it in fedora repositories so i suppose it has to be build from source. Pillow is the "new pil" library and using pillow the command "from PIL import ModuleX" does work from the qgis python console, but Lecos plugin installation still fails.

Is there a way to get Lecos recognize pillow library, or do I have to install the outdated pil library ?

thanks for your help


texts: QGIS 3x combining auxiliuary label rotation with the expresion

texts:
How can I combine these 2 features in the expression box?
   "auxiliary_storage_labeling_labeldistance"

and
 CASE
 WHEN "fid" = 31 Then 34.0
 WHEN "fid" = 54 Then -62.0
 WHEN "fid" = 35 Then -55.0
 WHEN "fid" = 36 Then -55.0
 WHEN "fid" = 29 Then 28.0
 WHEN "fid" = 27 Then -62.0
 WHEN "fid" = 33 Then 26.0
 WHEN "fid" = 16 Then -67.0
 WHEN "fid" = 6 Then 25.0
 WHEN "fid" = 93 Then 40.0
 end

When I merge them 2, the console says, that expression is invalid.
The basic problem is, that if I start to use the auxiliary rotation, I am losing all the expressions at the same time.


texts: MapServer runtime execution: Cannot evaluate expression, no item index defined error

texts:
I'm trying MapServer's run-time substitution to dynamically set my Mapfile's CLASS EXPRESSION values.

The expressions take the form of:

CLASS
    EXPRESSION (([my_attr_col] >= num_1) && ([my_attr_col] < num_2))
    ...
END

When manually setting the expression values (not using run-time substitution) it works. (And only seems to do so when I use no quotes at all.)

When trying to pass the expression in dynamically with run-time substitution, then I add a VALIDATION block:

VALIDATION
    'my_attr'    '^[\s\S]*$'
END

and I modify the CLASS EXPRESSIONs along these lines:

CLASS
    EXPRESSION '%my_attr%'
    ...
END

and I then append and send the information via the query URL.

However, I get this error:

msEvalExpression(): General error message. Cannot evaluate expression, no item index defined.

I've tried many permutations of quotes and regex filtering but I can't seem to figure out what would be accepted by MapServer.


texts: Scrolling features in QGIS following field?

texts:
I have a lot of features with labels associated. Because some of this features are wrong placed, I have to scroll them all. Usually I work with the attribute table open and select every single feature and zoom at it.
This is a quite annoying process but it's the only way that I found for now. Anybody knows if there is a tool/plugin that permits to scroll all the features in a layer following a field (like the primary key field)?

The plugin ItemBrowser does something similar for what I'm looking for but it needs a selection of the items. Selection means that if I have to move only one feature I have to deselect all the others, move the feature, and then reselect all (not so efficient and could cause a lot of errors)


texts:
I'm currently making a proportional symbol map for ArcMap. Safe to say I'm very bad at building queries. The data I'm using has "0" values in the attribute table. The proportional symbol tab is telling me to "exclude" 0 using the Data Exclusion Properties tab.

However, I have no idea what query I would need to build to make sure the "0" is excluded for both the value field and normalisation field. What would an exclusion query look like?


texts: GeoPortal - "Not an XML response received from the host " when adding a resource using URL protocol

texts:
I have a local directory I want to connect to. It has been configured in IIS to list the contents. If I visit it in my browser, it displays the contents.

When adding it as a resource, it passes the connection test fine. However I get a red text error saying "Not an XML response received from the host" when I use the URL protocol (Which is the correct one I take it?)

I tried using the WAF protocol, it adds the folder as a resource, however it shows up as a single entry, and the metadata looks like it was generated on it's own.

See screenshot:

The metadata it shows:

<rdf:RDF><rdf:Description rdf:about="http://localhost/MUNICIPAL/">
<dc:identifier>http://localhost/MUNICIPAL/</dc:identifier><dc:title>MUNICIPAL</dc:title>
<dct:references>http://localhost/MUNICIPAL/</dct:references></rdf:Description></rdf:RDF>

My question: How can I add a folder full of XML files to GeoPortal?

I am approving it after I add it.

EDIT: When I click preview, it shows this:

EDIT: Checked the gpt.xml log file. Seen an interesting error with web harvester:

Sep 8, 2014 2:53:45 PM com.esri.gpt.framework.context.ApplicationConfigurationLoader loadHarvesterConfiguration
INFO: [SYNCHRONIZER] Missing "webharvester.autoSelectFrequency" parameter. Default 3,600,000 will be used instead.
Sep 8, 2014 2:53:45 PM com.esri.gpt.framework.context.ApplicationConfigurationLoader loadHarvesterConfiguration
INFO: [SYNCHRONIZER] Missing "webharvester.watchDogFrequency" parameter. Default 60,000 will be used instead.
Sep 8, 2014 2:53:45 PM com.esri.gpt.framework.context.ApplicationConfigurationLoader loadHarvesterConfiguration
INFO: [SYNCHRONIZER] Missing "webharvester.baseContextPath" parameter. Harvest notification messages will be sent without information about harvest report.
Sep 8, 2014 2:53:45 PM com.esri.gpt.framework.context.ApplicationConfigurationLoader loadHarvesterConfiguration
INFO: [SYNCHRONIZER] Missing "webharvester.resource.autoApprove" parameter. Default false will be used instead.


texts: Convert NAM Grid (Lambert Conformal Conic) to MERRA2 Grid (lat/lon) using pyproj

texts:
I am attempting to convert NAM grid data into a simple lat/lon grid that my current MERRA2 data is on (NAM Grid Info). The NAM data is currently in lat/lon but on a Lambert Conformal Projection. I've tried a few other posts with regards to converting coordinate using pyproj (SE Post 1,SE Post 2), but with no luck. Below is some sample data and what I have tried.

from pyproj import Proj, transform
lat = np.array([12.19, 12.2199, 12.2497, 12.2795, 12.3092])
lon = np.array([-133.459, -133.354, -133.249, -133.144, -133.038])

nam = Proj('+proj=lcc +lat_1=12.190, +lon_0 = -133.459 +datum=NAD83 +no_defs')
merra2 = Proj('+proj=longlat +datum=WGS84 +no_defs')

x, y = nam(lon, lat, inverse=False)
x2, y2 = transform(nam, merra2, x, y)

And x2, y2 give me the same result as lat/lon. Am I missing something? I should note I am fairly new to geographic transformations. Let me know if you need any other information.


texts:
In QGIS on a new layer I'm trying to calculate the area of certain landtypes individually from a Grounds Maintenance layer. Previously I have used the following expression in Cadcorp SIS:

CalcInterior(Overlay('Grounds Maintenance',{Landtype$='Grass'}),'_area$',1)

Which works fine in Cadcorp,but I don't know how to write the same expression to work in QGIS.


texts: How to create a ArcGIS Online Organizations account?

texts:
Maybe this question is silly and I'm overseeing something. My question is: how to create a ArcGIS Online Organizations account?

We have to two ArcGIS for Desktop Basic licenses of the type Single Use and should have therewith access to an ArcGIS Online Organizational Account (if I believe the information of our dealer).

I found the option to Request Permissions but I don't get it from who and which permissions.

EDIT: I got a nice mail with the content:

  Your Permissions Request Has Been Received

  [...]

  The request will be reviewed by the administrator of the organization. You will receive an e-mail once your request has been processed.

The content is confusing as I should be already the administrator.


texts: Storing non-orthorectified images in GeoTIFF

texts:
I'm reading GeoTIFF and GDAL documentations and I'm trying to understand what type of transformations are supported.

GeoTIFF documentation says you can use ModelTransformationTag to define a transformation matrix, but later adds:

  For Baseline GeoTIFF, the model space is always 2-D, and so the matrix will have the more limited form

And then goes on to describe a matrix that only supports 2D affine transformations.

Later in the document it seems to imply that non-affine projections are only supported via use of tiepoints, and still would not work with Baseline GeoTIFF images.

  Case 5: The raster data cannot be fit into the model space with a
  simple affine transformation (rubber-sheeting required).
  Use only the ModelTiepoint tag, and specify as many
  tiepoints as your application requires. Note, however, that
  this is not a Baseline GeoTIFF implementation, and should
  not be used for interchange; it is recommended that the image be
  geometrically rectified first, and put into a standard projected
  coordinate system.

This seems to be confirmed in GDAL documentation, which describes SetGeoTransform, which is used to "set the affine transformation coefficients".

Does that mean there is no way of describing non-orthorectified images using GeoTIFF?

Is there some non-baseline GeoTIFF that supports non-affine (i.e. projective transformations)? If so, what software supports it?

Is there any other image format that natively supports projective transformations?


texts: How to Install GDAL with ECW Support on Windows

texts:
What do I download from http://www.gisinternals.com/stable.php for Win7-32 to cut 1:1M map sheets out of Australia.ecw file in GDA94/Geographic and save them as geotiff in WGS84/Geographic? I would like to automate the process to cut all of the 1:1M map sheets then reproject and convert them to geotiffs.

(I am thinking gdalinfo, gdaltranslate, gdalwarp but I am not clear about the ecw support and don't understand the gisinternals release numbers.)


texts: Specifying EPSG transformation method in GeoTools?

texts:
I am a software developer who knows a bit about Gis/Geodesy, but I am not an expert so I might misuse terminology somewhat.

I have been tasked to extend an app in Java and have a need to transform data from one Coord System to another. We started using GeoTools which lets us transform using EPSG codes.

    private MathTransform getMathTransformForEpsgCode(String epsgCode) throws FactoryException {
    MathTransform sourceToTargetConversion = null;
    try {
        // Get the result SourceCRS
        CoordinateReferenceSystem sourceCRS = CRS_FACTORY.createCoordinateReferenceSystem(epsgCode);
        // Create and apply Transform
        sourceToTargetConversion = CRS.findMathTransform(sourceCRS, TARGET_CRS);
    } catch (NoSuchAuthorityCodeException e) {
        log.info("No CRS found for code, assuming it is a transform", e);
        throw e;
    }

    return sourceToTargetConversion;
}

Which was then used something like the following

JTS.transform(coordinate, coordinate, returnedMathtransform);

HOWEVER, what we want to do is actually specify the math transformation method used (Note: I think this is akin I think to this question "How to specify transformation method when using ogr2ogr to reproject geometry?".

Is this possible with GeoTools?

Ideally specifying the math transform using an ESPG code (i.e. http://epsg.io/27700-5339)
Alternativly specifying 3 or 7 parameters

So I am a wee bit out of my depth here (but learning fast). I believe the above can be accomplished using a (licensed) ESRI java api. Can it be achieved easily in GeoTools or any other open source libraries, again ideally usable as a java api and not the command line.


texts:
Where can I find in depth documentation on modifying arcgis stylesheets? I have found some of the basics here :
Changing tool dialog appearance

I want to change the height on one element, the only element in my script. I have a filter value list of 82 counties and I would like to see more than 9 when the tool opens. Usually ESRI documentation is thorough unless it is a topic they expect everyone to know. Is this the case with xsl stylesheets? I think I could figure it out if I just knew what the element was called. does anyone have a suggestion?


texts: Low-pass filtering raster with a mask - matrix

texts:
I have to do low-pass filtering with a mask - matrix:

1 2 1
2 4 2
1 2 1

This is my code

import arcpy
import numpy as np
arcpy.env.workspace = r"C:\Users\ASUS\Documents\ArcGIS\Default.gdb"

raster = "Stadiony_Clip_Clip1"

nr = str("_05")
arcpy.CheckOutExtension("Spatial")
R = arcpy.Raster(raster)
lowL = arcpy.Point(R.extent.XMin, R.extent.YMin)
cellS = R.meanCellWidth

Rnp = arcpy.RasterToNumPyArray(R, nodata_to_value = -9999)
#A = [[1, 2, 1], [2, 4, 2], [1, 2, 1]]
Rnp = (1*Rnp[:-2,:-2] + 2*Rnp[1:-1,:-2] + 1*Rnp[2:,:-2] + 2*Rnp[:-2,1:-1] +
4*Rnp[1:-1,1:-1] + 2*Rnp[2:,1:-1] + 1*Rnp[:-2,2:] + 2*Rnp[1:-1,2:] + 1*Rnp[2:,2:])/16.

R = arcpy.NumPyArrayToRaster(Rnp, lowL, cellS, value_to_nodata = -9999)
R.save(raster + nr)
arcpy.CheckInExtension("Spatial")

print "KONIEC"

is my solution correct?

The output file gets a lock extension instead of tif, I think it's an error. How can I change this?


texts: Previewing File Geodatabase in QGIS Browser - Is it Possible?

texts:
Can I Preview/look at file geodatabase attributes in the QGIS browser if I have installed QGIS using OSGEO4W installer and FileGDB API? Currently I can add/edit file geodatabase feature classes to QGIS desktop, but the browser does not allow me to see preview/attributes/metadata. All I get is the Param tab previewing all the indexes, tables, yadda, yadda.

Thanks,
Mike


texts: I need some clarification on oracle sdo_gtype 4402. I understand that sdo_gtype is codded as DLTT. D: Dimension L : linear reference T T: TT identifies the geometry type. It's clear to me that first 4 indicated that the feature is 4-D. I am not sure about the second "4" and linear reference term. I need to create a spatial table using sdo_gtype 4402 that's why seeking some clarification on it.


texts: I need my Javascript browser application to convert/export KML documents into ESRI shapefiles. Is there any JS library that can do it? I tried googling, but no answers so far (at least, not in Javascript).


texts: Setting "REAL" work/scratch workspace in ESRI/ERDAS etc

texts:
I thought the issue was only in ARCGIS but it seems that ERDAS does the same thing. I have my C with 1.7TB and my D with 10TB and my scratch set to d:\scratch but still ERDAS (and ESRI) fills my profile with it's temp data...why? see the image below (in this example I have tried to set the temp space to c:\scratch)

In ARCGIS
 I have tried C:\Program Files (x86)\ArcGIS\Desktop10.0\Utilities\AdvancedArcMapSettings.exe--> run as admin-->System Paths-->last option "Output Temp file path" and set to D:\scratch.

I have tried "solutions" from http://forums.arcgis.com/threads/3300-ArcGIS-10-Default-Geodatabase

How can we resolve this?


texts: GeoServer + ImagePyramid: Request is empty where data exists

texts:
I have set up an image pyramid based on 5 levels of tiling. The data is served using the ImagePyramid plugin for GeoServer. When I do a layer preview the data is visualized just like I expect them to be. But when I zoom in a couple of levels, I get a white tile in return.

The GeoServer log tells me the following:

  2012-05-14 08:07:07,495 DEBUG [gce.imagemosaic] - Reading mosaic from file:/e:/GeoServer/data_dir/data/pyramid/2/2.shp
  2012-05-14 08:07:07,495 DEBUG [gce.imagemosaic] - Highest res 0.6400000000000005 0.6400000000000005
  2012-05-14 08:07:07,495 DEBUG [gce.imagemosaic] - The computed CropBoundingBox ReferencedEnvelope[534687.5 : 535156.25, 6225625.0 : 6226093.75] Doesn't intersect the coverage BoundingBox > ReferencedEnvelope[-1.0 : 0.0, -1.0 : 0.0] resulting in an empty request
  2012-05-14 08:07:07,495 DEBUG [gce.imagemosaic] - RequestedBBox empty or null
  2012-05-14 08:07:07,495 DEBUG [gce.imagemosaic] - Request is empty: RasterLayerRequest description:
      RequestedBBox=ReferencedEnvelope[534687.5 : 535156.25, 6225625.0 : 6226093.75]
                                   RequestedRasterArea=java.awt.Rectangle[x=0,y=0,width=768,height=768]
      RequestedGridToWorld=PARAM_MT["Affine",
          PARAMETER["num_row", 3],
          PARAMETER["num_col", 3],
          PARAMETER["elt_0_0", 0.6103515625],
          PARAMETER["elt_0_2", 534687.8051757812],
          PARAMETER["elt_1_1", -0.6103515625],
          PARAMETER["elt_1_2", 6226093.444824219]]
      ReadType=JAI_IMAGEREAD
  2012-05-14 08:07:07,495 DEBUG [gce.imagemosaic] - The response is empty. ==> returning a null GridCoverage

If I zoom in a couple of levels further, I get data again.

The extent of the data is [530000, 6220000, 540000, 6230000] and the projection EPSG:25832, which the GeoServer admin interface agree on.

Then what could be the reason for GeoServer to look inside the BBOX [-1, 0, -1, 0], when my data is inside the requested map?


texts:
I wrote simple code to check the xlwt module in Qgis. But I receive an error.

from xlwt import Workbook

wb = Workbook()
sheet1 = wb.add_sheet('Sheet 1')

sheet1.write(0,0,'Hello!')

wb.save("test12.xls")

ERROR:

Traceback (most recent call last):
  File "C:\PROGRA~1\QGIS3~1.0\apps\Python36\lib\code.py", line 91, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "<string>", line 10, in <module>
  File "C:\PROGRA~1\QGIS3~1.0\apps\Python36\lib\site-packages\xlwt\Workbook.py", line 710, in save
    doc.save(filename_or_stream, self.get_biff_data())
  File "C:\PROGRA~1\QGIS3~1.0\apps\Python36\lib\site-packages\xlwt\CompoundDoc.py", line 262, in save
    f = open(file_name_or_filelike_obj, 'w+b')
PermissionError: [Errno 13] Permission denied: 'test12.xls'

I need to solve this problem so that I can then take specific data from different attribute tables and write them to an Excel file.


texts: How to sum raster image values with vector overlay

texts:
I'm new to this and haven't been able to figure out how to do this on my own.

I have two files: a raster image of land use data and a vector-based shapefile of my territories.

The raster image has values like so: 11 = water, 21 = developed - open, 22 = developed - low, 23 = developed - medium, 41 = forest, 71 = grassland, 82 = cultivated crop, etc. Overall, I've got about 30 distinct values. The resolution of this image is pretty big (1m squares, I think) and the filesize is about 16 gigs.

I was to add a set of columns to my shapefile, one new column for each of the possible land use values. The column names would be like so: value_11, value_21, vaule_22, value_23, etc. Then I'd like to sum up the number of pixels in each shape that correspond to each value. So, if there are 1200 pixels that have a value of 11, then column value_11 would have 1200 in it.

I've used qgis to get the vector and raster files to load at the same time and they're georeffed correctly. It displays beautifully.

I'm not sure if I'm making this clear, but I'd like a count for every combination of each vector-based polygon and the unique values from the raster image within that polygon. If I need to do this in multiple passes, that's fine.

Using qgis, I was able to get average/median/mode/min/max of raster values within each polygon, but that's useless to me. I need the count of each distinct land use value. Eventually, I want to load the final shapefile into SQL Server 2008 and be able to write queries against it. I'd like to take a series of points, find out which polygon they are in, and determine how many water, developed, forest, etc. pixels are in that polygon.

Regarding software, I've played with several free/open source programs (such as QGIS, SAGA, etc), but I just don't know enough about what I'm trying to do to get anything to work. I can download whatever I might need, but I don't have a budget, so if ARCGIS or something expensive is the only way to do this, then I'm out of luck. Once I get the shapefile created, I can get it back into SQLServer and I can write the queries there to do what I want. I just need to get build the shapefile with the land use data added to it.


texts:
In my postgres database i have a set of points(lat, long) and i created the_geom index for that. Now i have to find the centroid of these points.

I Google a lot and i find that for closed objects(polygon etc) it possible to find the centroid. I am not sure how to find the centroid of points using postgis.

Also i found this:

SELECT ST_AsText(ST_Centroid('MULTIPOINT ( -1 0, -1 2, -1 3, -1 4, -1 7, 0 1, 0 3, 1 1, 2 0, 6 0, 7 8, 9 8, 10 6 )'));

But i didn't get any thing from this.

EDIT:

My requirement is:

Consider i have 5 places latlong, and i have to find which is center of this 5 places.

Help me to solve this. Thanks in advance.


texts:
According to release notes version 5.1.0 is capable of generating concave shapes. However, I cannot find any information on how to use this new feature.

Has anyone had any luck with this?


texts: Can I use the Google Elevation API for academic research (licensing)?

texts:
I would like to use the Google Elevation API to generate about 5 elevation maps for master's degree research on aircraft collision avoidance. The data would be meshed and used in simulation to calculate trajectory intersections of aircraft with the terrain.

I know the terms and conditions state that the data must be displayed in conjunction with a Google Map. If I publish my thesis with the meshed elevation data alongside a Google map of the area - is it allowed? What else would I need to do to be able to use the data for this purpose? The amount of data I need is within the usage limits.


texts: Changing symbol and color of vector Layer using QGIS 3 C++ APIs

texts:
I am trying to change the default symbol for the point/line/ polygon vector layer using QGIS3 C++ APIs.
I was able to render a vector layer with default symbols on QgsMapCanvas. When I try to change the default symbol of the point vector layer, somehow it is not changing. When I am trying to change color, it does not change color as expected, but the symbol color changes in a different color.
Here is the piece of code.
layer->startEditing();
qDebug()<< "is editable " << layer->isEditable();  // returns true
QgsFeatureRenderer * layerRenderer= layer->renderer();
QgsSingleSymbolRenderer *mSingleRenderer = QgsSingleSymbolRenderer::convertFromRenderer(layerRenderer);
QgsMarkerSymbol* symbol = new  QgsMarkerSymbol();
QgsStringMap mp;
mp[QString("name")]= QString("square");
mp[QString("color")]=QString("yellow");
symbol->createSimple(mp);
symbol->setColor(QColor(255,255,255));
mSingleRenderer->setSymbol(symbol);
layer->triggerRepaint();
mLayerTreeView->refreshLayerSymbology(layer->id());
layer->commitChanges();

variable-name layer defines QgsVectorLayer.


texts: Asking QGIS to provide upgrade for particular distro?

texts:
I effortlessly installed QGIS 2.4 and upgraded to 2.6, for my Linux Mint 16, because I just had to put http://qgis.org/debian in my repositories list and note that my distro was saucy.

2.8 has been out for some time now, and QGIS nicely updated to the new version for several different debian distros, but...not saucy.

I cannot figure out how to e-mail whatever volunteers do this and ask for saucy to be added.  I've attempted claiming to be other distros - they even updated for "precise" which is much older. But any older or newer and the upgrade won't go, it can't find various dependencies from other repositories.
Do I have to upgrade my Linux itself to get this QGIS upgrade to work?  It's only about a year old.


texts: Export.image.toDrive in TIFF format in GEE takes forever

texts:
I'm working with high resolution water data (30 m, Globas Surface Water Explorer). I need to ensure that all calculations are done in a certain projection (EPSG:4326) and after all calculations I want to extract to my Drive results that are aggregated to 1km, ~1.8 km and 5.5 km.
I can see results, and extraction is working without errors but it's been already more than 10 days and it keeps running. I have space in my GoogleDrive.
I tried following extraction:

var template_Global_01km = ee.Image("OpenLandMap/CLM/CLM_LST_MOD11A2-DAYNIGHT_M/v01").select('jan');
var globe_export = ee.Geometry.Rectangle([-179.9999999999999999, -90, 180, 90], 'EPSG:4326', false);
var dimensions_Global_01km = '43200x21600';
var maxPixels_numberD = 1e10; // max pixels for region reduction and download
var fileFormat_name = 'GeoTIFF'; // file format
Export.image.toDrive({
  image: lsm_Global_01km,
  description: 'lsm_Global_01km',
  folder: folder_name,
  dimensions: dimensions_Global_01km,
  region: globe_export,
  crs: template_Global_01km.projection().crs().getInfo(),
  maxPixels: maxPixels_numberD,
  skipEmptyTiles: false,
  fileFormat: fileFormat_name
});

I also tried to exract to Asset as thought it would be faster but it's the same - no errors and runs forever.
Script link: https://code.earthengine.google.com/?scriptPath=users%2Fmargaritachoulga%2FSTATIC_MAPS%3AStaticFields_water
Process Id 1: 66RXL5VFRB5O6UQF3JGE7X5G
Process Id 2: K65C7CW3QOYNSR5XGX6VC67U


texts: Select polygons which do not contain points

texts:
I am able to use the intersect tool to find polygons which overlap with points however I would like to find polygons which do not overlap with any points.

I think the Symmetrical Difference Tool would be an appropriate tool to do so, but I do not have the pro license to use it. I cannot use the Union tool either as it does not accept points as input.


texts: Radar Shadow for Sentinel 1 images in GEE?

texts: Is there a way to mask out areas of radar shadow in GEE? This should be produced during pre-processing I believe but I can't find a layer to add that would let me mask out the shadowed areas for Sentinel 1 in GEE?


texts: Individual houses in OS Zoomstack using QGIS - splitting terraced homes

texts: I am mapping the UK in QGIS using an OS zoomstack geopackage. I would like to look at the area and perimeter of properties in the UK but this is not possible with terrace and semidetached houses as there is no distinction of the number of properties are within that block. Does anyone know of a layer I can overlay to solve this? Is there a shapefile that could be added that will show each individual home?


texts: How to query a PostGIS database through Geoserver and OpenLayers

texts:
I've built a web application using OpenLayers, which is connected to my PostGIS database through Geoserver. Now that I'm able to successfully display all of the layers, I'd like to add querying/interactivity capabilities to the map. For instance, I'd like to create a dropdown menu of checkboxes so the user could "Choose X who have done Y", another dropdown to "Choose X who haven't done Y", etc. I'm having trouble figuring out how exactly to add in this functionality. Does anybody know any good resources/examples for how I might go about doing this?

Thanks for your help.


texts: gdal_grid with CSV file and VRT "No point geometry found on layer"

texts:
Given the following .vrt file which I manually created:

<OGRVRTDataSource>
    <OGRVRTLayer name="DGM5_BE">
        <SrcDataSource>DGM5_BE.txt</SrcDataSource>
        <SrcLayer>DGM5_BE</SrcLayer>
        <LayerSRS>EPSG:25833</LayerSRS>
        <GeometryType>wkbPoint</GeometryType>
        <GeometryField encoding="PointFromColumns" x="x" y="y" z="d"/>
    </OGRVRTLayer>
</OGRVRTDataSource>

and the associated DGM5_BE.txt - here are the first lines:

x;y;d
397200;5837250;59,15
397205;5837250;59,04
397200;5837245;59,03
397205;5837245;58,94
397195;5837240;58,93
397200;5837240;58,9
397205;5837240;58,86
397210;5837240;58,9
397190;5837235;58,98

I run the following command to convert the DGM5:

$ gdal_grid -of GTiff -ot Float32 -l DGM5_BE DGM5_BE.vrt DGM5_BE.tif

This fails with the following error message:

  ERROR 1: Failed to open datasource `DGM5_BE.txt'.
  No point geometry found on layer DGM5_BE, skipping.


texts: I have parcel data for a small town and I am trying to determine the impervious area of each parcel.  I have the impervious surfaces delineated into different shapes. How do I best subtract these impervious areas from the parcels?


texts: CartoDB data not displaying on Google basemap

texts:
I'm new to CartoDB but hoping to use it to display various types of real estate data on my website.  That said I've been playing around with some of the sample code (Google Maps V3 integration example - https://github.com/CartoDB/cartodb.js/blob/develop/examples/gmaps.html) and cannot get this data to display within my website (see test page - https://metropolisiq.com/MapTest).  I've tried the same code within Dreamweaver and it works fine.  Not sure what the issue is, but perhaps it relates to the differences between Dreamweaver and my website:

My site is encrypted (HTTPS)
My site was built using the Joomla CMS.

How do I solve this?