texts: Is ArcMap 10.0 Query Layer passing invalid SQL to Oracle 11gR2?

texts:
In two different environments I'm having trouble using ArcMap 10.0 Query Layers to view SDO_GEOMETRY (point) data in Oracle 11.2. This suggests I'm doing something wrong...

Points are in Oracle stored under SDO_GEOMETRY, 2001 type (2D point), 4326 SRID (WGS84 lat/lon). The table / column has a spatial index and is entered in user_sdo_geom_metadata.

I can successfully query the table in a way that uses its spatial index (SDO_WITHIN_DISTANCE), so I know this is functional.

When I add the table to ArcMap as a query layer the query validates fine, however nothing is displayed when I zoom to the layer. I have verified that the zoom extent is correct. If I try to show labels I see a drawing error:

Underlying DBMS error[ORA-29902: error in executing ODCIIndexStart() routine
ORA-22060: argument [2] is an invalid or uninitialized number
ORA-06512: at "MDSYS.SDO_INDEX_METHOD_10I", line 333
]

Not understanding how this involves me I decided to trace the database.

ArcMap makes the following query to Oracle:

select LOCATION_GEOMETRY
from (select * from TEST.BC_POINT_TEST_1) a
where sdo_filter(LOCATION_GEOMETRY,sdo_geometry( to_blob( :i1 ), :i2),'querytype=window') = 'TRUE'

I filled in the blanks to create the following query:

select LOCATION_GEOMETRY
from (select * from TEST.BC_POINT_TEST_1) a
where sdo_filter(LOCATION_GEOMETRY, SDO_GEOMETRY(2003, 4326, NULL, SDO_ELEM_INFO_ARRAY(1, 1003, 3), SDO_ORDINATE_ARRAY(53.276, 23.232, 53.587, 23.451)),'querytype=window') = 'TRUE'

Which worked fine. However, if I replace 4326 with NULL it throws a similar error (but not exactly the same) to what I'm seeing in ArcMap. This makes me think possibly ArcMap isn't passing the correct SRID for Oracle to know what to do. I set the data frame's coordinate system to WGS84 but this didn't make a difference.

The error reported through ArcMap references argument [2]. If this is actually the second bind parameter in the traced query that would be the SDO_ORDINATE_ARRAY, but as this is an array it would make no sense to complain that this is an invalid or uninitialized number.

If anyone has prior experience with this issue or suggestions on where to take it I would be most appreciative!

Edit: so far nothing has helped solve this problem, but I recently had an opportunity to test with 10.1 and everything worked as expected. This suggests there's nothing wrong with the Oracle instance.


texts: How to create shortest route and multiple routes service using openlayers?

texts:
I created a digital database of many important features and locations of my city. and I want to create a web application in which shortest distance between two locations can be calculated and displayed. Is it possible with openlayers as client scripting. My shapefiles are in postgis.
What I am expecting is :

The database contains all roads, point features etc.

When the user types a location, the required location name can be listed out from the database and then user can similarly select another location

Then my technique is to find the nearest road points for the two locations (may be using buffer method) and then

finding the shortest and possible routes between the two road points along the road vector features.
But I don't know how to write the code and algorithms.


texts: Changing user attribute format in dynamic translation using FME?

texts:
I'm a newbie in FME, and I've been asked to do a massive translation from MapInfo TAB to ESRI File Geodatabase. There are a lot of files so I don't know the schema.

The merge feature reader and dynamic writer look to work well. But for some features, there is an attribute that the writer cannot write, normally logic or date attribute.

I'm trying to change those attributes format to text.

Is there anyway to do this in FME?


texts: QGIS2.0 on Ubuntu 12.04 cannot see the hints when hover over

texts:
I am running QGIS2.0 on Ubuntu server 12.04. After installing some libraries for Grass and Orfeo toolbox (I don't remember exactly which ones), I cannot see the hints when i hover over a menu bar or a layer. I can see only the colored frame of the hint without the text. Has anyone experienced something similar?
Thanks!


texts: Copy attribute data from one polygon shapefile to a multiple point shapefile

texts:
I have two shapefiles:

Polygons of parks in Africa. Each park (polygon) has a row of data regarding it's size, name, date of establishment etc.
Random points in parks in Africa. Here, I generated 500 random points in each park.Therefore this attribute table has 500 (random points) x 224 parks = 96096 rows.

However, in my random points shapefile there is no indication of which points are in which park. I would like to know how I can get the park's name from the polygon shapefile to match each point ID that is in that park? I see that my random points shapefile has a column called CID which goes from 0 - 224 (Indicating the parks that the points are in), but I did not specify this and I would like park name to be there instead (even though I know it will be repeated 500 times for each matching row).

I am using ArcMap 10.2


texts: How should I make use of railway routes from Google Maps API?

texts: I have a specific requirement wherein I have to display or highlight only railway stations in the map rather than the whole map. How can this be achieved? This is the code I tried:


texts: How to fix Error:  could not access file "$libdir/librouting"?

texts:
I am trying to run the sql file of pgrouting-> routing_core.sql in mac. However it displays the following error

psql:routing_core.sql:32: Error:  could not access file "$libdir/librouting": No such file or directory

How can I remove this?


texts: Understanding ArcGIS field calculator behaviour with geometries?

texts:
In Getting Properties of Next Record Using Calculate Field Tool in ModelBuilder? the OP asks for suggestion on calculation of distances between pairs of points in a feature class.

I thought of alternative solution using field calculator with a very basic logic:

for first record memorise point geometry from Shape field into global variable pG and return negative distance
for other records calculate distance between current Shape and one in memory, assign pG = current Shape, return distance.

Expression:

i=0
def FirstOrOthers(shp):
 global pG,i
 if i==0:  i, dist, pG = 1, -1, shp
 else: dist, pG=shp.distanceTo(pG),shp
 return dist
#------
FirstOrOthers( !Shape! )

Output shown in field DISTANCE:

When I changed expression slightly, to memorise first point of point geometry:

i=0
def FirstOrOthers(shp):
 global pG,i
 if i==0:  i, dist, pG = 1, -1, shp.firstPoint
 else: dist, pG=shp.distanceTo(pG),shp.firstPoint
 return dist
#------
FirstOrOthers( !Shape! )

It worked as expected:

It means that field calculator refuse to compute distance between 2 point geometries, but does it for point geometry and a point. This is risky, because point does not have spatial reference defined. This is why in scripts I always use 2 point geometries to compute distance between 2 points.

Any thoughts on that odd behaviour of field calculator?

I tested using ArcGIS Desktop 10.3.


texts:
I am trying to overlay points on a polygon, but identicalCRS complains the reference system is not the same. It actually is the same:

> polygon@proj4string
CRS arguments:
 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
> points@proj4string
CRS arguments:
 +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
>
> identicalCRS(polygon@proj4string, points@proj4string)
Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘proj4string’ for signature ‘"CRS"’
> identicalCRS(polygon, points)
[1] FALSE

Bug or feature? Any workaround?


texts: Zonal Statistics not running in QGIS 3.12

texts:
I am working on Nighttime light data from NOAA. The study area is India at state level. After masking the raster, when I am trying to run the zonal statistics with raster and vector polygon layer (screen shot attached), the small window of zonal statistics just flickers and does not run. Though the log section says that 'Algorithm zonal statistics finished' but it does not run. It does not add any new columns in the vector polygon attribute table. The attribute table is as it is as it was before. I have been trying to run this for last two days, but it does not help.
What could be the possible reason for it? I am currently using QGIS version 3.12.


texts:
I can't add the TauDEM Tools.tbx to ArcMap 10.6.1. Windows 10 64bit.

  .... taudem tools.tbx is either not designed to run on Windows or it
  contains and err. ...... error status  0x000012f


texts: Is it possible to temporarily suspend snapping with a shortcut Key in Qgis, like its possible holding down space in Arcmap?


texts: Mapping multiple datasets using Google Fusion Tables - possible?

texts:
I've used ArcGIS for a few years but only recently discovered Google Fusion Tables. 2 very simple questions but not ones I've been able to find any information on online - hoping someone in this community might know:

Is it possible to plot more than 1 kml file at a time and have all map properly? (e.g., is it possible to recreate how multiple layers would be
displayed in ArcGIS)
Is it possible to plot different shapefile types at the same time? e.g., plot a polygon, then plot polylines superimposed on top of the polygon? or do Fusion Tables only support 1 kml type?

Trial-and-error attempts at getting both 1 & 2 to work suggest it doesn't work, but I was hoping to get a more concrete and documented "yes" or "no" - or whether there is some formatting trick that can be used?


texts: Not scaling and/or rotating when exporting map using ArcObjects?

texts:
I am writing a console application using ArcGIS (engine or desktop). I'm trying to rotate and scale a map based on a polygon, and I'm not having any luck at all. With some tinkering, I can usually do one or the other but not both.

Anybody have any suggestions? I feel as though I've tried every variation in terms of rotating the map and adjusting the scale either using the Imap interface, or the IActiveView interface.

    Private Sub Zoom2AndExportJpgs(mapDocument As IMapDocument, imagePath As String)
    'Get the relevant strip maps for the route
    Dim map As IMap = mapDocument.Map(0) '(assumed that there is only one map in the mxd)
    Dim pageLayout As IPageLayout = mapDocument.PageLayout
    Dim pageAV As IActiveView = pageLayout
    Dim mapAV As IActiveView = map
    Try
        'get a list of strip map objects... EF object containing a list of the 4 corners of a polygon
        Dim stripmaps As List(Of MRLAStripMapExtent) = getStripMaps()
        Dim filename As String = ""

        'get all the associated stripmaps for the route  for each scale
        If Not stripmaps Is Nothing Then
            For i = 0 To stripmaps.Count - 1
                Dim pgon As IPolygon = geometryHelper.CreatePGON(stripmaps(i).RectLowerLeftX, stripmaps(i).RectLowerLeftY, _
                                                                                stripmaps(i).RectLowerRightX, stripmaps(i).RectLowerRightY, _
                                                                                stripmaps(i).RectUpperLeftX, stripmaps(i).RectUpperLeftY, _
                                                                                stripmaps(i).RectUpperRightX, stripmaps(i).RectUpperRightY)
                filename = imagePath & stripmaps(i).ZoomScale & stripmaps(i).PageNumber & ".jpg"
                CreateJpgFile(filename, mapAV, pageAV, pgon, stripmaps(i).ANGLE, stripmaps(i).ZoomScale)

            Next i
        End If

    Catch ex As Exception
        publishError.PublishError(publishError.ErrorSeverity.majorStopRouteProcessing, ex)
    End Try

End Sub

Private Sub CreateJpgFile(ByVal sFilename As String, ByVal pAV As ESRI.ArcGIS.Carto.IActiveView, _
                          ByVal pPageActiveView As ESRI.ArcGIS.Carto.IActiveView, _
                          ByVal pgon As IPolygon,
                          rotationAngle As Double, _
                          mapScale As Double)

    Dim pExport As ESRI.ArcGIS.Output.IExport
    pExport = New ESRI.ArcGIS.Output.ExportJPEG

    pAV.Extent = pgon.Envelope
    pAV.ScreenDisplay.DisplayTransformation.ScaleRatio = mapScale
    pAV.ScreenDisplay.DisplayTransformation.Rotation = rotationAngle '0
    pAV.Refresh()

    'Export Image....
    pExport.ExportFileName = sFilename
    pExport.Resolution = 96
    Dim exportRECT As ESRI.ArcGIS.esriSystem.tagRECT
    exportRECT = pPageActiveView.ExportFrame
    Dim pPixelBoundsEnv As ESRI.ArcGIS.Geometry.IEnvelope
    pPixelBoundsEnv = New ESRI.ArcGIS.Geometry.Envelope
    pPixelBoundsEnv.PutCoords(exportRECT.left, exportRECT.bottom, exportRECT.right, exportRECT.top)
    pExport.PixelBounds = pPixelBoundsEnv
    Dim hDC As Integer
    hDC = pExport.StartExporting
    pPageActiveView.Output(hDC, pExport.Resolution, exportRECT, Nothing, Nothing)

    'Finish writing the export file and cleanup any intermediate files.
    pExport.FinishExporting()
    pExport.Cleanup()
End Sub


texts: Python in ArcGIS, CalculateField_management, updating with single value

texts:
I have a float field which I want to simply update all values with 1. I can easily do this within the field calculator but I have started using python today, what is the syntax to simply update a field with the value 1. I'm stuck on the expression bit of the script:

      arcpy.CalculateField_management("ROCHDALE_TEST_NET","FRMEAS",??????


texts: Error while seeding via MapProxy from GeoServer


texts: QGIS Layout: Filter Nth Row in Attribute Table

texts:
I have an attribute table in a Qgis Layout and only want to see the first line.
Can I do that?
It's a workaround for overlapping labels.
I want an attribute table with the 1. Line, then a table with the 2. line... and so on
Only features intersecting the atlas are shown.
Most of the time it's not more than 1 to 10 at the same time.
I tried various things with arrays in the filter expression but can not find a solution.
In the "Filter with" expression I probably need something like "the first result of the array that is intersecting the atlas" or something like that!


texts: What is the relationship between Autodesk and the MapGuide Open Source project?

texts:
In the year 2006 Autodesk started MapGuide as an open source project (MGOS). At the same time they began to develop the commercial version of MapGuide (MapGuide Enterprise MGE /Autodesk Infrastructure Map Server AIMS).

With respect to the software development: What is the relationship between Autodesk and the MapGuide Open Source project today (and in the future?)? How does Autodesk support the MapGuide Open Source community today? Are developers from Autodesk writing code for MGOS? How benefits MGOS from new features of AIMS?

Regards Jens


texts:
I am using geotools, my goal is to print locations (of grid) in normal way.

How to  prepare style in a correct way?

 public Style createLineStyle2() {
        StyleBuilder sb = new StyleBuilder();
        FilterFactory ff
                = CommonFactoryFinder.getFilterFactory(GeoTools.getDefaultHints());

        Font font = new StyleBuilder().createFont("Verdana", 12.0D);
        Font fonts[] = {font};

        AnchorPoint anchorPoint = sb.createAnchorPoint(0.5, 0.5);

        PointPlacement pointPlacement = sb.createPointPlacement(anchorPoint, null, sb.literalExpression(0));
        LinePlacement linePlacement = sb.createLinePlacement(0.0);
        TextSymbolizer text = sb.createTextSymbolizer(null, fonts, sb.createHalo(),
                sb.attributeExpression("Location"), pointPlacement, "Location");

        Style style = sb.createStyle();
        LineSymbolizer lineSymbolizer = sb.createLineSymbolizer(Color.red, 1.0);
        Rule rules[] = {sb.createRule(text), sb.createRule(lineSymbolizer)};
        FeatureTypeStyle fts = sb.createFeatureTypeStyle("Feature", rules);
        style.featureTypeStyles().add(fts);

        return style;
    }


texts: Shapely/Pyproj converts all polygon coordinates to inf

texts:
After some libraries got updated, I can't convert a polygon between CRS in shapely/pyproj. Everything seems to be done according to the documentation.

Pyproj version 2.4.2.post1, Shapely version 1.7.0

Here's the test code:

from shapely.ops import transform as transform_
from shapely.wkt import loads
import pyproj
from functools import partial
print(pyproj.__version__)
obj = loads('POLYGON ((101.6862738336452 3.1597141052896, 101.6862738336452 3.176801618144858, 101.6706701928565 3.176801618144858, 101.6706701928565 3.1597141052896, 101.6862738336452 3.1597141052896))')
crs_from = pyproj.CRS.from_string('epsg:4326')
print(repr(crs_from), crs_from.__class__)

crs_to = pyproj.CRS.from_string('epsg:3857')
print(repr(crs_to), crs_to.__class__)

print(obj.__class__)
print(obj.wkt)
print(transform_(partial(pyproj.transform, crs_from, crs_to), obj).wkt)

Output:

$ python3 /tmp/crs.py
2.4.2.post1
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich
 <class 'pyproj.crs.CRS'>
<Projected CRS: EPSG:3857>
Name: WGS 84 / Pseudo-Mercator
Axis Info [cartesian]:
- X[east]: Easting (metre)
- Y[north]: Northing (metre)
Area of Use:
- name: World - 85°S to 85°N
- bounds: (-180.0, -85.06, 180.0, 85.06)
Coordinate Operation:
- name: Popular Visualisation Pseudo-Mercator
- method: Popular Visualisation Pseudo Mercator
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich
 <class 'pyproj.crs.CRS'>
<class 'shapely.geometry.polygon.Polygon'>
POLYGON ((101.6862738336452 3.1597141052896, 101.6862738336452 3.176801618144858, 101.6706701928565 3.176801618144858, 101.6706701928565 3.1597141052896, 101.6862738336452 3.1597141052896))
POLYGON ((inf inf, inf inf, inf inf, inf inf, inf inf))

After changing the code to use pyproj.Transformer... still getting the same results.


texts:
I'm running QGIS 2.8 Wien on Windows 8.1. I got this error:

  The program can't start because qgis_core.dll is missing from your
  computer. Try reinstalling the program to fix this problem.

I got the same error last week, uninstalled and reinstalled QGIS and it worked. But there must be a better way to fix the error than reinstalling QGIS every time it bombs out.


texts: OSGeo4W updating is failing - OpenCL.dll was not found

texts:
Windows 10. I've not used my 3 installations of QGIS for 6 weeks maybe. I updated them via OSGeo4w but I get the error "The code execution cannot proceed because OpenCL.dll was not found. Reinstalling the program may fix the problem".

OSGeo4w gives me the error

Package: qgis-common qgis-common.bat exit code -1073741515.
Package: qgis-ltr-common qgis-ltr-common.bat exit code -1073741515

Attempting to open any of the QGIS installations gives me the error

Reinstalling the programmes produces the same errors. Can anyone help please?


texts: Add spatial csv-layer in QGIS3.10 in a standalone script

texts:
I want to add a csv-layer by PyQGIS in a standalone script. This is my code:
uri = 'file:///C:/Users/sml/Desktop/point.csv?delimiter=,&crs=epsg:4326&xField=x&yField=y'
vlayer = QgsVectorLayer(uri, "point", "delimitedtext")
print(vlayer.isValid())
QgsProject.instance().addMapLayer(vlayer)

The same code works fine in the QGIS Python console, but it can't work in a standalone script. I don't know why.


texts:
I have written a script to go through each layer in MapInfo (from 1 to 168), turn layer 1 on, style it, label it, print, turn layer off, then the same for layer 2,3,4 etc. The printing works, however I think it's going a bit too fast for MapInfo/the printer and will crash every 25-30 pages, to fix this I am removing my printed layers and tweaking my variables.

I've had a google and from what I can see there is no wait function in MapBasic, is there something else I can do? E.g having it perform a pointless task just to kill a bit of time between pages so the printer does not get overloaded?


texts: Watershed other than stream ( between two points)

texts: Is there any method available to demarcate watershed for a location other than stream? I mean, if I want to find out the watershed for a place, where no stream available but, sheet of water can flow (like a causeway / Irish bridge). Thanks....


texts:
I am trying to do the pagination using my own geometry filter. However, I am always getting this: KeyError: 'id'. I can run the example/pagination.py successfully without any error.

Here is my pagination.py:

import os
import requests

# our demo filter that filters by geometry, date and cloud cover
from CCIL_filters import geometry_filter, champaign_county

#item_types can be the following:
#"PSScene4Band"
#"PSScene3Band"
#"REScene"
#"REOrthoTile"
#"Sentinel2L1C"
#"PSOrthoTile"
#"Landsat8L1G"

print(geometry_filter)
# Search API request object
search_endpoint_request = {
    "item_types": ["PSScene4Band"],
    "filter": geometry_filter
}

################################
session = requests.Session()
session.auth = (os.environ['PLANET_API_KEY'], '')

# Create a Saved Search
saved_search = \
    session.post(
      'https://api.planet.com/data/v1/searches/',
      json=search_endpoint_request)

# after you create a search, save the id. This is what is needed
# to execute the search.
saved_search_id = saved_search.json()["id"]

# What we want to do with each page of search results
# in this case, just print out each id
def handle_page(page):
    for item in page["features"]:
        print item["id"]

# How to Paginate:
# 1) Request a page of search results
# 2) do something with the page of results
# 3) if there is more data, recurse and call this method on the next page.
def fetch_page(search_url):
    page = session.get(search_url).json()
    handle_page(page)
    next_url = page["_links"].get("_next")
    if next_url:
        fetch_page(next_url)

first_page = \
    ("https://api.planet.com/data/v1/searches/{}" +
    "/results?_page_size={}").format(saved_search_id, 5)

# kick off the pagination
fetch_page(first_page)

Could someone give me some help on this?


texts: Watermarking TIFF with pillow in Python breaks georeferencing

texts:
I am trying to automate the watermarking of thousands of georeferenced TIFF files with a Python script. The files have related tfwx world files and aux.xml files. When I run the code the spatial referencing seems to be corrupted. The cell size of the output TIFF seems to be 200x larger. I can then run arcpy.Rescale_management to recover their proper size, but they are still incorrectly placed. I am a noob, so I apologize if this is coming across as vague or ignorant. I have tried several workarounds without luck and I just want to make sure I'm not missing something easy.

Here is a stripped down version of the code I am using:

from PIL import Image, ImageDraw, ImageFont
import os, sys, arcpy
from arcpy import env

env.workspace = r'D:\mbrush\scratch\tools\CalTrans\testing'
workspace = env.workspace
env.overwriteOutput = True

# copywrite logo and name
text = u'\u00a9' + ' Company Name'

# iterate through tiffs and watermark
for file in os.listdir(workspace):
    if file.endswith('.tif'):

        # open image
        im = Image.open(workspace + '\\' + file)

        # open image for drawing
        drawing = ImageDraw.Draw(im)

        # determine height and width of image
        width = im.size[0]
        height = im.size[1]

        # determine where to place the text
        posx = width - (500)
        posy = height - (80)
        pos = (posx, posy)

        # set font size
        fontsize = 33

        # define outfile to overwrite original
        img_out = workspace + '\\' + file

        # choose font
        font = ImageFont.truetype(r'C:/Windows/Fonts/Arial/arial.ttf', fontsize)

        # create text shadow
        shadow = (pos[0]-1, pos[1]-1)
        shadow2 = (pos[0]-2, pos[1]-2)
        shadow3 = (pos[0]+1, pos[1]+1)
        shadow4 = (pos[0]+2, pos[1]+2)
        drawing.text(shadow, text, fill=(200,200,200), font=font)
        drawing.text(shadow2, text, fill=(255,255,255), font=font)
        drawing.text(shadow3, text, fill=(200,200,200), font=font)
        drawing.text(shadow4, text, fill=(255,255,255), font=font)

        # draw watermark
        drawing.text(pos, text, fill=(3, 8, 12), font=font)

        # save image
        im.save(img_out)


texts: Calculating length of road over different kinds of land cover using ArcGIS Desktop?

texts:
I have road network (shapefile polyline) and landcover (polygon).

How can I calculate the length of road that cross different kinds of landcover？

According to the answer, I write the following code.

6 kinds of landcover and 35 roads

# Import arcpy module
import arcpy

for j in range(1, 7):
    land_cover = "Land_" + str(j)
    for i in range(1, 36):

        # Local variables:
        in_land = "U:\\4_work_space\\LandCover_shp\\" + land_cover + "_Dissolve.shp"
        in_line = "U:\\4_work_space\\Line\\Line_" + str(i) + ".shp"
        out_shp = "U:\\4_work_space\\output\\" + land_cover + str(i) + ".shp"
        out_txt = "U:\\4_work_space\\output\\" + land_cover + str(i) + ".txt"

        print(i)

        # Process: Intersect
        arcpy.Intersect_analysis([in_land, in_line], out_shp, "ALL", "", "INPUT")

        # Process: Add Geometry Attributes
        arcpy.AddGeometryAttributes_management(out_shp, "LENGTH_GEODESIC", "METERS", "", "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]")

        arcpy.ExportXYv_stats(out_shp, "LENGTH_GEO", "SPACE", out_txt, "NO_FIELD_NAMES")


texts: In R, How to read in GeoJSON file, modify, then export back as new GeoJSON file

texts:
Familiar with R, but not GIS.
I have a GeoJSON file with tens of polygons, each of which has thousands of coordinates, and I want to read it into R, perform some validation of the structure, fix where needed, and export it back to file.
I've used FROM_GeoJson (package geojsonR) to read in the file. It represents it as a list of lists. I can easily manipulate this structure, but there isn't, for example, a corresponding TO_GeoJson. What is the easiest way to export this to file?
I read this question: How can I export a GeoJSON file from R?, and it shows how to generate a GeoJSON from a data frame. But I would have to convert my list of lists to a data frame. I'm unclear how to specify the coordinates when there is more than one point, however.
I'm hoping for a straightforward way that even a novice GIS/R person would easily know so I don't have to spend several hours learning how to do it myself.


texts: Exporting data by month/year using ModelBuilder?

texts:
Using ArcGIS Desktop 10.6.1. I have a large dataset that I need to view sequentially by month AND year. Instead of using Select by Attribute for each month/year combo (i.e. February 2017-May 2019 for 28 total feature classes).

Can this be done in ModelBuilder?

Data is currently in a file geodatabase and I'm looking to export each date/month as a feature class.


texts: How to limit the search for facilities to only 1 unit with PG Routing? ... Closest Facility in PostGIS?

texts:
I'm working in PGROUTING - POSTGIS
I have my road network and two sets of points -firestation and the entire cityblocks from town (in diferent tables, 1 for origin, and 1 for destination).

The plan is to calculate the optimal routes from the fire station to the blocks. Later it would filter for maximum time. I have already achieved this first point just explained using algorithm pgr_dijkstra (many-to-many).

DO
$$DECLARE
mzn_o integer[] = (SELECT ARRAY (select origen from origin_table));
mzn_d integer[] = (SELECT ARRAY (select destino from destination_table));

BEGIN
create table my_routes_table as
SELECT seq::integer, path_seq::integer, start_vid::integer, end_vid::integer, concat(start_vid::integer, end_vid::integer) as id, node::integer,edge::integer, a.cost::double precision, agg_cost::double precision, b.class_id::integer, b.the_geomx::geometry
    FROM pgr_dijkstra(
    'SELECT id, source, target, traveltime as cost, traveltime_rev as reverse_cost FROM my_network_table', mzn_o, mzn_d1,directed := false) a left join my_network_table b on a.edge=b.id;
END$$;

I would like to improve this operation, because in the above code you get all possible routes from all sources to all destinations. I would like that when a block finds its nearest fire station, it stops searching. In other words, each station is assigned only the closest station and not all of them. In this way, achieve a higher calculation performance and with less cost, since the combinations would be significantly limited. That's the idea, something similar or identical to what is achieved with the Network Analyst Closest Facility tool in ArcMap, where it seems that you can enter the maximum number of facilities to calculate the route.https://desktop.arcgis.com/en/arcmap/latest/extensions/network-analyst/closest-facility.htm

Maybe a picture for understanding:

My purpose is not to insert a filter after main operation, but to integrate that idea so that it results in an improvement in the performance of the operation that causes a considerable decrease in the amount of operation and can work with large networks and improve the capacity of the computer.

I think I have already tried everything within my reach and my technical level, and my eyes are bleeding from looking at possibilities on the internet. In Postgis or PG Routing I have not been able to find anything related to Closest Facilities (if to calculate with service area and pgr_drivingtime algorithm), something that limited the calculation to find only the closest route.
I thought that with the pgr_astar algorithm it could be achieved since its concept is to find the best route and reduce possibilities but it also returns the same as with Dijkstra in my tests. I have also tried to make a cost matrix and try filtering with the minimum cost min (cost), I thought that this way I would filter out the final nodes not needed, as explained in this example https://www.qgistutorials.com/en/docs/3/origin_destination_matrix.html. But it's not working for me.

SELECT * FROM pgr_dijkstraCost(
    'select id, source, target, cost from my_network_table',
    (SELECT ARRAY (select source from origin_table)), (SELECT ARRAY (select target from destination_table)));

select source, min(cost) as shortest_path from costs_table group by source;


texts:
I've implemented a openlayers bbox strategy on a vector layer. I have my server return a nice jsonp object wrapped in the callback given by the request.

The problem is that during normal use of the map I get al lott of errors like:

TypeError: OpenLayers.Protocol.Script.registry.c18 is not a function

where c18 is the number of the request. The map keeps working, but these show up as actual errors in Firebug. This happens very often just playing around for a minute results in a dozen or so of these errors. In the 'NET' tab I can see that it seems to happen every time the map interface issues two requests for data.

GET poi?bbox=-0.597037...cript.registry.c18, 200 OK, localhost, 101.8KB, 127.0.0.1:80, 250ms

GET poi?bbox=-6.831778...cript.registry.c19, 200 OK, localhost, 102.3KB, 127.0.0.1:80, 421ms

I guess the callback corresponding to the first one of the set of two gets deleted by the second request, so when the first bit of data returns, you get the error.
So, why does it fire two requests? and how to fix this. My guess is that it has something to do with resizing the map and zooming at the same time and both events result in a request??

Here is how I configured the layer. It's the only layer with a BBOX strategy:

var in_options = {
    'internalProjection': config.projection,
    'externalProjection': config.displayProjection
};

// add a layer to the map
this.layer = new OpenLayers.Layer.Vector("POI",{
    strategies: [
        new OpenLayers.Strategy.BBOX({
            resFactor: 2,
            active: false,
            autoActivate: false,
            ratio: 1.3
        }),
        new OpenLayers.Strategy.Cluster({
            distance: this.clusterDistance,
            threshold: this.clusterThreshold
        })
    ],
    protocol: new OpenLayers.Protocol.Script({
        url: 'http://localhost/tests/poi',
        callbackKey: 'callback',
        format: new OpenLayers.Format.GeoJSON(in_options),
    }),
    styleMap: new OpenLayers.StyleMap({
        "default": style,
    })
});
this.layer.strategies[0].calculateBounds = this.myCalculateBounds;
this.layer.strategies[0].activate();

Edit
I've looked at the layer properties, there are two listeners registered. One for the BBOX strategy and one for the cluster. None of these invoke each other. So I looked at the source of the BBOX strategy. The TriggerRead Method of the OpenLayers.Strategy.BBOX class has one argument: options. If you console log this one, it's just the event object send by the layers eventHandler. It's not an opion object associated with the BBOX instance. That also explains why the setting noAbort wasn't working.
I changed this line:

if (this.response && !(options && options.noAbort === true)) {

into this:

if (this.response && !(this.noAbort && this.noAbort === true)) {

and ofcourse set the noAbort: true in the layer strategy

 new OpenLayers.Strategy.BBOX({
        resFactor: 2,
        active: false,
        autoActivate: false,
        ratio: 1.3,
        noAbort: true,
    })

Did I discover a bug in OpenLayers?


texts: How to download satellite image from planet.com using API key

texts:
I would like to know what to do (step by step in simple words) and how to get data from www.planet.com

There is tutorial
https://www.planet.com/docs/api-quickstart-examples/
but anyways I have problem..

So I have API key and Trial 14 days access to the platform. I tried to use Python command line to get satellite image that I chosed, but actually I am not good at Python at all and I don't know how to build request.

export PL_API_KEY=mykey (does not work either)


texts: Null values coerced to zero (0) when exporting to csv or txt

texts:
I extracted values from several rasters, and in some cases, the value is null. However, when I export the attribute table to a csv or a txt file (from ArcMap), the null values are converted to 0. Is there any way to prevent this from happening?
These zeros are a real downer, as throw off my analysis,  which is done in R.


texts: How to invoke a script in the currently running QGIS instance from the command line?

texts:
From the Windows Command Prompt (or Batchfile) I want to find the running QGIS instance  with an open project. Afterwards a script (eg. zoom to something) should be automatically started.
I've tried it with qgis.bat --code

That works fine, but always starts a new QGIS instance. - I don't want that!

In future I want to register a Windows Protocol for something like link2qgis:layer="layer".featid="123".action=zoomto

When clicking on that url in some third party application, the QGIS instance should be maximized and the map canvas should be zoomed to the feat.

Maybe someone has already done something similar and can help me with that.


texts: Obtain coordinates and corresponding pixel values from GeoTiff using python gdal and save them as numpy array

texts:
How can I obtain projected coordinates as well as the actual pixel values at those coordinates from a GeoTiff file and then save them into a numpy array? I have arsenci020l.tif file, and its coordinates are in meters. Below is the abridged output of gdalinfo I ran on it.

~$ gdalinfo arsenci020l.tif
Driver: GTiff/GeoTIFF
Files: arsenci020l.tif
       arsenci020l.tfw
Size is 10366, 7273
Coordinate System is:
PROJCS["Lambert Azimuthal Equal Area projection with arbitrary plane grid; projection center 100.0 degrees W, 45.0 degrees N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Lambert_Azimuthal_Equal_Area"],
    PARAMETER["latitude_of_center",45],
    PARAMETER["longitude_of_center",-100],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-6086629.000000000000000,4488761.000000000000000)
Pixel Size = (1000.000000000000000,-1000.000000000000000)
...

There was a similar question here about obtaining lat/long coordinates from tiff (Obtain Latitude and Longitude from a GeoTIFF File) and the answer showed how to obtain only top left x and y pixel coordinates. I need to obtain ALL projected pixel coordinates as well as get the pixel values and save them in a numpy array. How can I do it?


texts:
Some Python software I have written generates raster files via GDAL. I would like the software to create an XML metadata file for each raster so that I can include details about how the raster was made, who made it, and so-on.

I understand that ISO 19115 is the international geospatial metadata standard. I really have no idea where to start with generating a metadata file, though. Can GDAL do this from scratch? Or would it be better to work with a dedicated XML library in Python?


texts: How to create a profile with profile points at discrete altitude intervalls, e.g. each 10 m

texts:
In both SAGA and QGIS theres a profile-from-line-function that does a nice job of creating profiles. But there's no way to control the altitude spacing.

Any ideas for workflows that would allow me to achieve a profile (river as line-shape on a DEM-layer), with individual profile-points situated only at discrete altitude-intervalls of 10 m ? Of course profile points will need to contain information on segment-length and altidude...

Any help appreciated!


texts: Openlayers 3 set icon pointerEvents to none for click through

texts:
I've read through the OL3 API documentation without finding how or where to disable click on an icon so the map underneath gets the click instead.

Perhaps it's not possible?
Or does it go in ol.Feature, ol.style.Style or ol.style.Icon somehow?

code example...

xFeature = new ol.Feature({
    geometry: new ol.geom.Point([0,0])
    });
xStyle = new ol.style.Style({
    image: new ol.style.Icon(({
        anchor: [40,40],
        anchorXUnits: 'pixels',
        anchorYUnits: 'pixels',
        src: 'image.png'
        }))
    });
xFeature.setStyle(xStyle);

Update 1: erilem has helped me understand that clicks are ignored by default. Further in my code I am using forEachFeatureAtPixel which is obviously intercepting the click.

map.on('click', function(evt){
    var feature = map.forEachFeatureAtPixel(evt.pixel,
        function(feature, layer){
            return feature;
            });
    if(feature){
        /* do whatever */
        }
    });

My Solution add uuid to icon feature & ignore that in click iteration

xFeature = new ol.Feature({
    geometry: new ol.geom.Point([0,0]),
    uuid: 1
    });
xStyle = new ol.style.Style({
    image: new ol.style.Icon(({
        anchor: [40,40],
        anchorXUnits: 'pixels',
        anchorYUnits: 'pixels',
        src: 'image.png'
        }))
    });
xFeature.setStyle(xStyle);

map.on('click', function(evt){
    var feature = map.forEachFeatureAtPixel(evt.pixel,
        function(feature, layer){
            if( feature.get('uuid') != 1 ){
                    return feature;
                    }
            });
    if(feature){
        /* do whatever */
        }
    });


texts: ArcGIS Calculate Near by looping through and matching .shp names from two different directories

texts: I have several thousand shapefiles in one directory (direct1) and several thousand other shapefiles in another directory (direct2).  The files in direct2 will have a match in direct1 and will have the same file name.  I need to do a near analysis by calculating the the nearest point from direct2 to the nearest point in direct1 for files that share the same name.  I found this example that is pretty much what I want to do, but it is for ArcGIS 9.3.  I am working with ArcGIS 10.3 and looking for some type of updated code and or something that will perform faster!


texts:
I'm a bit of GIS rookie but I seem to have found my way around most of the intermediate functions of QGIS that I need for my work over the last year or so. What I want is really simple, but multiple searches have returned nothing. I'm using a base map from quick map services and I want to be able to manually trace over all of the water courses in a certain area in order to highlight them as well as annotate a number of  other features such as location of assets etc. I know this can be done in the print composer but I want to be able to add it as a layer so that I can produce multiple prints of different "zoom frames".

Is there a tool that I'm missing that allows you to manually annotate and add as a layer; essentially paint in GIS 2.14.2?

I'd also be interested to know if there is a layer already out there for mapped watercourses in the uk.


texts: Linear Referencing Complex Intersections/Road Segments

texts:
What is the proper way to handle complex/branching intersections in linear referencing? I have a road shapefile I was going to linear reference, but I have some segments that are very complex. I have tried ignoring their complexity and locating them on my route like they were a simple line, but I ended up with problems since the tool doesn't distinguish between the parallel segments and assigns them to segments seemingly at random. I have tried researching solutions, but am stuck. See below my failed attempt at locating features along the route.

The points' intended locations are shown below spread across the road segments at my points of interest.

Do I have to simplify this segment of the route somehow or are there other steps/tool settings that will enable the proper location of the points on my existing route?


texts: ogr2og2 fails with memory allocation error post v2.0 for MSSQLSpatial

texts:
There are a few of us encountering "Memory allocation" errors when using ogr2ogr to import spatial (shp/gpx) data into a MSSQLSpatial databse. According to this post: Writing shapefile to SQL Server database with ogr2ogr Error
 the problem is due to using the latest version of the gdal libraries. The "answer" is to use a previous version, but I find this unacceptable for a number of reasons not the least of which is that an OSGeo4W installer will not enable you to install the 1.x version required and therefore it seems that its no longer possible to install older versions.

Im using gdal 1.11.2 with native SQL 11 driver successfully; however, 2.1.0 fails with same SQL driver.

Has anybody successfully had this work? What are the parameter/driver changes required to get this to work?


texts: How to extract a non contiguous multilinestring portion between 2 points?

texts:
I have a non contiguous multilinestring.
I try to cut these multilinestring with points (see schema below)
I tried something like the code below, but the function St_LineSubString and ST_LineLocatePoint only work with linestring and do not work not with multilinestring...
How can I cut these non contiguous multilinestring?

DECLARE
    r1 record; -- Roads table (1 multilinestring per road))
    r2 record; -- Points table
    p1 geometry; -- point
    p2 geometry;  --point
    res varchar;
BEGIN
    FOR r1 in SELECT * FROM roads ORDER BY road_name ASC
    LOOP
        p1=ST_StartPoint(r1.geom);
        FOR r2 IN SELECT * FROM points WHERE road_name = r1.road_name ORDER BY "PRD" ASC
        LOOP
            p2=r2.geom;
            SELECT ST_AsText(St_LineSubString(r1.geom,ST_LineLocatePoint(r1.geom,p1.geom),ST_LineLocatePoint(r1.geom,p2.geom))) INTO res;
            RAISE NOTICE '%',res;
            p1=p2;
        END LOOP;
    END LOOP;
END;


texts: Creating polygon grid in ArcGIS model using ET Geowizard?

texts:
When bring in the vector grid extent tool from ET Geowizards surface tools into a model An error pops up:

  ERROR 00824 The tool is not licensed

Outside of the model the tool works fine. What I need is a grid of 100 x 100 m polygons which covers the extent of another layer.


texts: How to create shaded relief from xyz file?

texts: I have a DEM (xyz format) which I want to apply to a "flat" map texture to create a shaded relief, preferably with software that is available for free.


texts:
I want to be able to build a measure in ArcGIS (ArcMap).

(1) I have some shapefiles of sovereign states by centuries. These are polygons with the name of the state as ID. (2) Moreover, I have a shapefile of current European regions at NUTS 2 level.

I would like to get a measure of the average number of States that existed in a given NUTS 2 region, over centuries. I have tried to build such a measure for days now.

My bet was to build this measure by using the tool Union. Take two maps of sovereign states for 1000 and 1100, for example, unionize them.
Union will give me a new attribute table with the polygons non-intersecting at t=1000, t=1100 and those intersecting in between.

What I want is a measure that will tell me how many States did rule in a given polygon.
Project the two sovereign states shapefiles for 1000 and 1100. Some polygons will likely maintain their borders, hence the measure should give 1 for the latter.
In other cases, the polygon in the first period will be crossed by another polygon. This will give three polygons: two that belong to the same states as before and another one (the intersection) which was belonging to state A in period t= 1000 and now to B in period t=1100. The measure should give 2 for this polygon.
Then I was thinking about using something like Spatial Join to be able to get the average count of polygons in a NUTS 2 region.

However the very first problem is that the attribute table that I get with GIS after using the Union tool as polygons with the State related at t=1000 and t= 11000 but I don't ever know how I can convert these different names to count (and a number).


texts: Highlight all features associated with selected feature via relationship in ArcObjects C#

texts:
I am trying to highlight those features which are associated (via Relationship) to the selected feature in ArcObjects C#.
public void IRelationship_Example(IFeature feature)
{
    //assumes that only one RelationshipClass exists for the Origin feature class
    IEnumRelationshipClass enumRelClass = feature.Class.get_RelationshipClasses(esriRelRole.esriRelRoleAny);
    IRelationshipClass relClass = enumRelClass.Next();
    //if a feature with no Relationships established has been selected, exit
    if (relClass == null)
    {
        return;
    }
    ESRI.ArcGIS.esriSystem.ISet relSet = relClass.GetObjectsRelatedToObject((IObject)feature);
    relSet.Reset();
    //If an Attributed Relationship does not exist, exit
    if (relClass.IsAttributed != true)
    {
        return;
    }
    IFeature destinationFeature = (IFeature)relSet.Next();
    while (destinationFeature != null)
    {
        IRelationship relationship = relClass.GetRelationship((IObject)feature, (IObject)destinationFeature);
        IRow row = (IRow)relationship;
        object attributeValue;
        if (row.get_Value(0) == null)
        {
            attributeValue = "0";
        }
        else
        {
            attributeValue = row.get_Value(0);
        }
        MessageBox.Show("Destination OID: {0}  Origin OID:  {1}  Attribute value:  {2}", relationship.DestinationObject.get_Value(0), relationship.OriginObject.get_Value(0), attributeValue);
        destinationFeature = (IFeature)relSet.Next();
    }
}

I am newbie in ArcGIS Products.

As shown in above image, once, I select the feature, I can able to highlight all associated features.
For selecting feature I am using following Code:
newselection = true;  //set newselection to true on every mouse click.
HighlightClickedFeature(Layer, 22, activeView);

private void HighlightClickedFeature(IFeatureLayer layer, int featureOid, IActiveView activeView)
{
        IFeatureClass featureClass = layer.FeatureClass;
        IFeature feature = featureClass.GetFeature(featureOid);
        IGeometry geometry = feature.Shape;

        ISpatialFilter spatialFilter = new SpatialFilter();
        spatialFilter.Geometry = geometry;
        spatialFilter.GeometryField = featureClass.ShapeFieldName;
        spatialFilter.SpatialRel = esriSpatialRelEnum.esriSpatialRelIntersects;
        spatialFilter.WhereClause = "OBJECTID = " + featureOid;

        IFeatureSelection featureSelection = layer as IFeatureSelection;
        if (featureSelection != null)
        {
            if (newselection == true)
            {
                //featureSelection.Add(feature);
                //featureSelection.CombinationMethod = esriSelectionResultEnum.esriSelectionResultAdd;
                //featureSelection.SelectionSet.Combine(selectionSet, esriSetOperation.esriSetUnion, selectionSet2);
                featureSelection.SelectFeatures(spatialFilter, esriSelectionResultEnum.esriSelectionResultNew, false);
                newselection = false;
            }
            else
            {
                featureSelection.SelectFeatures(spatialFilter, esriSelectionResultEnum.esriSelectionResultAdd, false);
            }
            activeView.Refresh();
        }
    }


texts:
I am following this question with the first answer Dissolve/aggregate polygons with ogr2ogr or gpc? where I am trying to dissolve my shapefile down to 1 row by the geometry field

ogr2ogr C:\Users\name\Desktop\output.shp C:\Users\name\Desktop\All_constraints_v5.shp -dialect sqlite -sql "SELECT ST_Union(geometry) FROM input"

ERROR 1: In ExecuteSQL(): sqlite3_prepare(SELECT ST_Union(geometry) FROM input):
 no such table: input

right so then I did

C:\Program Files\GDAL>ogrinfo -so C:\Users\rizagha\Desktop\All_constraints_v5.shp
INFO: Open of `C:\Users\rizagha\Desktop\All_constraints_v5.shp'
      using driver `ESRI Shapefile' successful.
1: All_constraints_v5 (Polygon)

So my question is

"SELECT ST_Union(geometry) FROM input" does input need to be my shp full path variable? what is the exact syntax
I cannot find out what the geometry column is called.

I ran this:

C:\Program Files\GDAL>ogrinfo -geomfield -fields=YES C:\Users\rizagha\Desktop\All_constraints_v5.shp
INFO: Open of `C:\Users\rizagha\Desktop\All_constraints_v5.shp'
      using driver `ESRI Shapefile' successful.
1: All_constraints_v5 (Polygon)

but it doesn't tell me what column name the geom is


texts: TextFormating.NewLine not honored in Field Maps

texts:
I have an arcade expression that removes null fields from my pop up in AGOL Map Viewer Beta. I need the same pop up settings to also be honored in Field Maps. When using the TextFormatting.NewLine or "\n" all my fields get concated into one line in Field Maps. Here are the expressions I've tried using:
var returnstring = TextFormatting.NewLine

    for (var i in $feature) {
        if (!isEmpty($feature[i])){
          returnstring = returnstring + i + ": " + $feature[i] + TextFormatting.NewLine
      }
    }
    return returnstring

And this expression:
var result = "";
if (!IsEmpty($feature.ABUNDANCE)) {
    result += "\n" + "Abundance: " + $feature.ABUNDANCE;
}
if (!IsEmpty($feature.OCCSTATUS)) {
    result += "\n" + "Occurence Status: " + $feature.OCCSTATUS;
}
if (!IsEmpty($feature.COMNAME)) {
    result += "\n" + "Common Name: " + $feature.COMNAME;
}

return result;

In field maps the popup window looks like
`Abundance: 2 Occurence Status: 2 Common
Name: Example`

Anyone know of a workaround? I'm trying to just hide null values from my pop up and not have blank spaces between rows.
Desired results:
`Abundance: 2
 Occurence Status: 2
 Common Name: Example`


texts: Adding raster catalog to MXD using arcpy.mapping?

texts:
Using the bit of code below adds the raster catalog footprint polygon layer.

mxd = arcpy.mapping.MapDocument("Current")
df = arcpy.mapping.ListDataFrames(mxd, "*")[0]
addOrtho2011 = arcpy.mapping.Layer(r"P:\Mapping\Data\Orthos.gdb\Orthos_2011")
arcpy.mapping.AddLayer(df, addOrtho2011, "TOP")

I can't figure out what else I have to do to add the actual ortho part of the raster catalog.  Any suggestions?


texts:
I would like to to iterate through polygon features in a shapefile and spatially join them with points. I am able to do this with model builder, as well as when the polygons are in separate shapefiles (using ListFeatureClasses) but I would like to know how to do this with many polygons in a single shapefile using Search Cursor.

import os
import sys
import arcpy
from arcpy import env

points = r'C:\Temp2\points.shp'
polygons = r'C:\Temp2\polygons.shp'

output = r'C:\Temp2\output'

bCursor = arcpy.da.SearchCursor(polygons, '*')

for row in bCursor:
   print row
   arcpy.SpatialJoin_analysis(row, points,output+'%i%.shp',"JOIN_ONE_TO_ONE","KEEP_ALL","", 'INTERSECT', '', '')

print 'Done'

Error (when run in WinPython Spyder):

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\WinPython-32bit-2.7.10.3\python-2.7.10\lib\site-packages\spyderlib\widgets\externalshell\sitecustomize.py", line 790, in runfile
    execfile(filename, namespace)
  File "C:\WinPython-32bit-2.7.10.3\python-2.7.10\lib\site-packages\spyderlib\widgets\externalshell\sitecustomize.py", line 77, in execfile
    exec(compile(scripttext, filename, 'exec'), glob, loc)
  File "C:/GISDATA/python/batchSpatialJoin.py", line 26, in <module>
    arcpy.SpatialJoin_analysis(row, points,output+'%i%.shp',"JOIN_ONE_TO_ONE","KEEP_ALL","", 'INTERSECT', '', '')
  File "C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\analysis.py", line 481, in SpatialJoin
    raise e
RuntimeError: Object: Error in executing tool

Error (when run in ArcMap Python Window):

Runtime error
Traceback (most recent call last):
  File "<string>", line 15, in <module>
  File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\analysis.py", line 481, in SpatialJoin
    raise e
RuntimeError: Object: Error in executing tool

I have tried this on different computers, ArcMap 10.1 and 10.4, with the same errors


texts: How to determine the right strategy for WFS requests with Getotools?

texts:
I found that the strategy is selected in the determineCorrectStrategy() method in org.geotools.data.wfs.internal.WFSClient.

I don't know how reliable the suggestions are for all the strategies but e.g. for mapserver I found multiple services that do not contain "mapserver" or "map=" but something like "mapserv" and are therefore not found if the strategy is set to "auto" (default).

Are there more reliable indicators to choose the right strategy?


texts: OpenLayers 3, problem with GeoJSON file?

texts:
I have two layers:

var layer = new ol.layer.Tile({source: new ol.source.OSM()});

and

var layer_g=new ol.layer.Vector({
    projection : 'EPSG:4326',
    style: new ol.style.Style(
          {
              fill: new ol.style.Fill({color: [0, 255, 0, 0.2]}),
              stroke: new ol.style.Stroke({color: 'black', opacity: 0.1})

          }
            ),
    source: new ol.source.GeoJSON({
    url: 'file.geojson',
    projection : 'EPSG:3857'
    })
  });

Where the file.geojson is:

{
  "type": "LineString",
  "crs": {
    "type": "name",
    "properties": {
      "name": "urn:ogc:def:crs:OGC:1.3:CRS84"
    }
  },
  "coordinates":
  [
[21.919146, 51.962142]
.....
,[21.955264, 51.955445]
,[21.955894, 51.955762]
,[21.956169, 51.955924]
,[21.957933, 51.956962]
  ],  "properties": {
    "distance": "9.783024",
    "description": "To enable simple instructions add: 'instructions=1' as parameter to the URL",
    "traveltime": "653"
    }
}

I add these two layers to map:

var map = new ol.Map({
  layers: [layer,w_gminy],
  target: document.getElementById('map'),
  view: view
});

But only OSM layer is visible.

Is the problem in a structure of the GeoJSON file?


texts: Trying to display some data in OpenLayers - it's not erroring, just not loading :(

texts:
I'm trying to display some KML data on a map, using OpenLayers.

Nothing complex - but it's just not rendering.

Here's my javascript code

var map;

function init(element, kml) {
    map = new OpenLayers.Map(element, {
        projection: 'EPSG:3857',
        layers: [
            new OpenLayers.Layer.Google(
                'Google Streets', // the default
                { numZoomLevels: 20 }
            ),
            new OpenLayers.Layer.Vector('Kml data')
        ],
        center: new OpenLayers.LonLat(153.132, -26.987)
            // Google.v3 uses web mercator as projection, so we have to
            // transform our coordinates
            .transform('EPSG:4326', 'EPSG:3857'),
        zoom: 10
    });

    // Grab the features.
    var features = convertKmlToFeatures(kml);

    // Display it.
    map.layers[1].addFeatures(features);
}

function convertKmlToFeatures(kml) {
    var kmlFormat = new OpenLayers.Format.KML();
    return kmlFormat.read(kml);
}

Basically,

Create a map with two layers.
Convert the Kml to features.
Display the features on the map.

I just don't see what I've done wrong :(

Oh - on a side note, the KML data is fine, because when I tried it to visualise it on this website, it worked.


texts: Script tools will not run and error is blank?

texts:
When I try and run any script tools from ArcToolbox, the red X icon shows up next to the input features line, but there is no error when I hover over the X. If I proceed with the tool and click OK, I get a blank error message with the yellow exclamation point error (see screen capture below). I have tried multiple script tools and am having the same error with each. An example tool that is not working is "Generate Points Along Lines". The input features are shapefiles with projected coordinate systems.

I am running ArcMap 10.4.1 (Basic) on Windows 10 Pro.

I have no experience with Python, but based on some reading I have done online, I am wondering if this error is Python related.


texts:
Taking a long shot but does anyone have an idea of whether GDAL or OGR has a function to determine a shapefile or layer geometry type using python? I specifically want the output to print "POINT" for example.

I understand that the geometry can be printed to well-known text format and have the code to complete that, but specifically require the string format.


texts: UNION two LEFT JOINs as Virtual Layer in QGIS

texts:
In QGIS 2.18, there are two layers, e.g. "Locations" and "Connections"

"Locations" has values

Id | Dimension
---------------
 1 |    4
 2 |    8
 3 |    2

"Connections" maintains attributes

Origin | Destination | Value | Distance_KM
-------------------------------------------
   1   |      2      |  500  |     30
   1   |      3      |  100  |     20
   2   |      1      |  100  |     10
   2   |      3      |  300  |     10
   3   |      1      |  100  |     40

I want to create a Virtual Layer with the following Attribute Table. Where "In" correspond to "Destination" from "Connections" and "Out" to "Origin" accordingly.

Id | Dimension | In_Value | In_Count | In_Dist | Out_Value | Out_Count | Out_Dist
----------------------------------------------------------------------------------
 1 |     4     |    200   |     2    |    50   |    600    |    2      |     50
 2 |     8     |    500   |     1    |    30   |    400    |    2      |     20
 3 |     2     |    400   |     2    |    30   |    100    |    1      |     40

I can achieve the result that I strive for separately with two queries.

  Query 1

SELECT C.Destination, SUM(C.Value) AS In_Value, COUNT(C.Destination) AS In_Count, SUM(C.Distance_KM) AS In_Dist
FROM Connections AS C
GROUP BY C.Destination

  Query 2

SELECT C.Origin, SUM(C.Value) AS Out_Value, COUNT(C.Origin) AS Out_Count, SUM(C.Origine_KM) AS Out_Dist
FROM Connections AS C
GROUP BY C.Origin

Nevertheless, there should be only one query that solves my issue, is not it? I tried this but no success. What can be a problem?

SELECT L.Id AS Id, L.Dimension AS Dimension, C.In_Value, C.In_Count, C.In_Dist, C.Out_Value, C.Out_Count, C.Out_Dist
FROM   Locations AS L
LEFT JOIN (
           SELECT C.Destination,
                  SUM(C.Value) AS In_Value,
                  COUNT(C.Destination) AS In_Count,
                  SUM(C.Distance_KM) AS In_Dist
           FROM Connections AS C
           GROUP BY C.Destination
        ) ON L.Id = C.Destination
LEFT JOIN (
           SELECT C.Origin,
                  SUM(C.Value) AS Out_Value,
                  COUNT(C.Origin) AS Out_Count,
                  SUM(C.Origine_KM) AS Out_Dist
           FROM Connections AS C
           GROUP BY C.Origin
        ) ON L.Id = C.Origin

Moreover, I would like to know if I am eligible to add a second LEFT JOIN ON to the query with already existing LEFT JOIN ON.

References:

Two SQL LEFT JOINS produce incorrect result
Multiple left joins on multiple tables in one query


texts:
I am trying to make a map using r.sunmask.datetime in QGIS GRASS. Now i have to following problem. The Netherlands is largely under sealevel, and I am starting to believe that files, which heights are below 0 are crashing.

What can be the reason of this constant crashing?

Remember, I also have parts of maps which are above sealevel which just work correctly.


texts: Converting point kriging surface to small polygon data set?

texts:
I have a dataset showing soil samples as points and the nutrient contents at each of these sample locations.

I have created an interpolation surface that represents the distribution of nutrients across my study area.

Within my study area are a number of smaller blocks or polygons. I want to either create new polygons or update an attribute of the existing polygons to represent the mean results from the interpolated surface generated from the point samples. The polygon blocks are equivalent to the vegetative block shown in the image below

A result similar to those generated from areal interpolation would be great but the aeral tool can only handle aeral interpolated.

Currently using ArcGIS 10


texts: Geometry generator and calling an expression from a custom function

texts:
I am using geometry generator symbol layer type and some expressions that makes the result a line that is shorter then the original line.

The expression works and looks like this:

make_line(line_interpolate_point($geometry,10),line_interpolate_point($geometry,35))

However I would like to make my own function that makes it easier to do this operation. Inside the function I am trying to use the make_line expression but this results in an error in the expressions panel when using the function.

My function looks like this:

from qgis.core import *
from qgis.gui import *

@qgsfunction(args='auto', group='Custom')
def smallerline(geometry,value1, feature, parent):
tmpline = make_line(line_interpolate_point(geometry,value1),line_interpolate_point(geometry,value1))

return tmpline

I'm using the function in the expression tab but there's an eval error that says "Eval Error:
global name 'make_line' is not defined"


texts:
I am trying to enable audit logs in geoserver. Apparently there should be a monitoring/ directory within the data directory. It should then contain the file monitor.properties.

I cannot find this directory in any of the later versions. If i try to manually create, restart then load some WFS layers no log files are being generated.

Anyone able to give advice?

Its very sparsely documented. I was going from this example https://geoserver.geo-solutions.it/edu/en/adv_gsconfig/auditing.html and http://docs.geoserver.org/latest/en/user/extensions/monitoring/audit.html

I've checked versions 2.10, 2.11 and 2.13. None of these versions had a monitoring/ directory.

Here is the content of the monitor.properties file I created:

audit.enabled=true
audit.path=[PATH TO FOLDER]
audit.roll_limit=20

EDIT: ok... Ignore this... I don't have the Monitoring extension installed. I didn't realise I needed it!


texts:
I am dealing with landsat ETM+ images after 2003. While filling the line drop outs usiung QGIS Fill Nodata,  I keep getting this error, ERROR 1: TIFFOpen:7948_0fill_y_work.tif: Permission denied

What could be the Problem


texts:
I am trying to open a dialog box which shows up only raster datasets so that I can select one and create Hill shade based on it. But I am facing a run time error for Dialog box when I click on the addin button. Could any one please help me out regarding this

Protected Overrides Sub OnClick()
    '
    '  TODO: Sample code showing how to access button host
    '
    My.ArcMap.Application.CurrentTool = Nothing

    Dim pMxDoc As IMxDocument = My.ArcMap.Document
    Dim pMap As IMap = pMxDoc.FocusMap
    Dim bObjectSelected As Boolean
    Dim pGxDialog As IGxDialog = New GxDialog
    'only show raster datasets
    Dim pRasterFilter As IGxObjectFilter = New GxFilterRasterDatasets
    'Layers and document declaration
    Dim pGxObjects As IEnumGxObject = Nothing
    Dim pGxDataset As IGxDataset
    Dim pRasterLayer As IRasterLayer

    Try
        With pGxDialog
            .AllowMultiSelect = False
            .StartingLocation = "C:\Users\rxm160030\Downloads\Lab12\Data"
            .Title = "Select Raster datasets to Add to Map"
            .ObjectFilter = pRasterFilter
            .DoModalOpen(My.ArcMap.Application.hWnd, pGxObjects)
        End With

        'If no file selected, exit
        If bObjectSelected = Nothing Then Exit Sub

        'Get first dataset
        pGxDataset = pGxObjects.Next
        'no dataset selected
        If pGxDataset Is Nothing Then Exit Sub
        'iterate through remaining datasets
        Do Until pGxDataset Is Nothing
            pRasterLayer = New RasterLayer
            'assign dataset to rasterLayer
            Dim pInputRaster As IRaster = pRasterLayer.Raster
            'assign dataset name to rasterLayer
            pRasterLayer.Name = pGxDataset.Dataset.Name
            ' Do something with the featureLayer
            Dim choice1 As String = InputBox("Enter Azimuth")
            Dim choice2 As String = InputBox("Enter Altitude")

            Dim pSurfaceOp As ISurfaceOp = New RasterSurfaceOp
            Dim pEnv As IRasterAnalysisEnvironment = pSurfaceOp
            Dim pWSF As IWorkspaceFactory = New RasterWorkspaceFactory
            Dim pWS As IWorkspace = pWSF.OpenFromFile("C:\Users\rxm160030\Downloads\Lab12\Data", 0)
            pEnv.OutWorkspace = pWS
            Dim pOutput As IGeoDataset = pSurfaceOp.HillShade(pInputRaster, choice1, choice2, True)

            Dim pOutLayer As IRasterLayer = New RasterLayer
            pOutLayer.CreateFromRaster(pOutput)
            pOutLayer.Name = "HillShade"
            pMap.AddLayer(pRasterLayer)
            pMap.AddLayer(pOutLayer)
            ' Go to next dataset
            pGxDataset = pGxObjects.Next
        Loop
    Catch ex As Exception
        MsgBox("BtnAdd_Click Exception: " & ex.Message)
    Finally
        ' Release ArcObjects objects
        pGxObjects = Nothing
        pRasterFilter = Nothing
        pGxDialog = Nothing
        pGxDataset = Nothing
        pRasterLayer = Nothing
        ' Call garbage collection
        GC.Collect()
        GC.WaitForPendingFinalizers()
    End Try
    pMxDoc.ActiveView.PartialRefresh(esriViewGeography, Nothing, Nothing)
End Sub


texts: How to calculate a value based on another field under arcpy.da.UpdateCursor and cursor.updateRow?

texts:
I am trying to do a calculations under arcpy.da.UpdateCursor. I want to calculate the values for the records that have FIPS=06037. I want to take the values from "kg" column and populate the new value by multiply by 3 in "kg_new" column.

I think I am not using the right syntax under

row_field_Name[0] = "!kg!"*3

and for the next FIPS in the shapefile , i.e. FIPS=06308,

I want to take the "kg"value and multiply by 5.

    # Set the workspace environment
arcpy.env.workspace = r"E:\Test_Errors_onE\myPath"
arcpy.env.overwriteOutput = True

infc = "myFile.shp"
#Add a new filed named
field1 = ["kg"]

field_Name = "kg_new"
field_Nametype = "DOUBLE"
field_Namename = arcpy.ValidateFieldName(field_Name)
#arcpy.AddField_management(infc, field_Name, field_Nametype, "", "", "")

cursor = arcpy.da.UpdateCursor(infc, field_Name, "\"FIPS\" = 06307")
for row_field_Name in cursor:
    row_field_Name[0] = "!kg!"*3
    cursor.updateRow(row_field_Name)

cursor = arcpy.da.UpdateCursor(infc, field_Name, "\"FIPS\" = 06308"")
for row_field_Name in cursor:
    row_field_Name[0] = "!kg!"*5
    cursor.updateRow(row_field_Name)


texts:
While coding to obtain individual layers' attributes (field) list, I am stymied by the code.

I can identify a specific layer:
MyLayer = QgsMapLayerRegistry.instance().mapLayersByName("TestPoint")
which returns on "print MyLayer":
qgis._core.QgsVectorLayer object at 0x000000000FCC6730
from which I infer I am correctly using MyLayer as a Layer object/pointer (or whatnot...)

However, when I run this:
fields = Mylayer.pendingAllAttributesList()
I force this error:
AttributeError: 'list' object has no attribute 'pendingAllAttributesList'

Following this statement, eventually I would run a "For field in fields...:" indent block, to iterate through the fields, but I can't even get past my error.

I get same/similar errors for other attribute functions within the API, using versions provided in code samples, and best I can tell I am duplicating relevant python code within the samples found.  Somehow, I am missing something important about the way the Attributes work, or the coding, or the list or object.

I am using the recently downloaded QGIS 2.6.0 package for windows.


texts:
What I'd like to try and do is intersect some point buffers with census tracts, maintain their ID's and the tracts attributes.

I found a sample from http://www.macwright.org/2012/10/31/gis-with-python-shapely-fiona.html but can't exactly substitute the cascaded union for an intersect. I'm also having trouble getting the intersection function to import from shapely.ops.

from shapely.geometry import Point, mapping, shape
from fiona import collection
#from shapely.ops import intersection
import shapely.ops

bufSHP = 'data/h1_buf.shp'
intSHP = 'data/h1_buf_int_ct.shp'
ctSHP  = 'data/nyct2010.shp'
with collection(bufSHP, "r") as input:
    schema = input.schema.copy()
    with collection(intSHP, "w", "ESRI Shapefile", schema) as output:
        shapes = []
        for f in input:
            shapes.append(shape(f['geometry']))
        merged = shapes.intersection(ctSHP)
        output.write({
            'properties': {
                'uid': point['properties']['uid']
                },
            'geometry': mapping(merged)
            })

I've also got the code and sample shapefiles up here https://github.com/nygeog/questions/tree/master/shapely_intersect

UPDATE:

So for some reason using Shapely, the larger buffers (1 km) I lose some of the features (census tracts) in the intersect.

Also, since the Shapely responder said the OGR command is good to use I started looking for into that. Anyway, I was able to test cd'ing to the directory and it works great. However, I've been trying to experiment with it so I can loop through the many sized buffers and dif. census years. However, why is it that I cannot get the following script to work. Can I not substitute paths for the folder 'data' in the ogr2ogr response below?

import os
wp = '/Volumes/Echo/GIS/projects/naas/tasks/201411_geoprocessing/data/processing/'
folder = wp+'test'
theGDALcmd = 'ogr2ogr -sql "SELECT ST_Intersection(A.geometry, B.geometry) AS geometry, A.*, B.* FROM nyct2010 A, hr1000m B WHERE ST_Intersects(A.geometry, B.geometry)" -dialect SQLITE '+folder+' '+folder+' -nln hr1000m_int_ct10_ogr_pytest'
os.system(theGDALcmd)

The error I get is

sh: ogr2ogr: command not found

However, when I go to the terminal, I have access to ogr2ogr.

UPDATE2:

Mike T include some code so I can run this ogr2ogr command in python.


texts:
"I would like to calculate the Euclidian distance from a residence to a public service of interest e.g. the distance from a residence to a green space (as the crow flies) in QGIS. After browsing the Internet it seems that I need to rasterize my vector before I can use the Proximity tool provided by GDAL.
When I try to make the conversion from vector to raster my thought was to take the shp.file which contains the polygons. The function "rasterize" tells me to chose an "output file for rasterized vectors (raster)" but I have no idea what to chose here.
Should I create a raster for this? And how do I this?
I'm using QGIS 1.8 Mac and GDAL 1.9.
Thanks in advance."

Ok, it seems like I have not done my homework good enough. Let me start over: I have a vector layer with points from which I want to calculate the distance to vector layers with lines or polygons. Still I want to calculate the Euclidian distance (as the crow flies) from the points in the first layer to the lines or polygons in the other layers: The distance from a residence (point) to a stream (line) or forest (polygon). Additionally, I would like to put the distances in the attribute table of the vector layer with points in a new column. Sorry for the inconvenience.


texts: Capturing shape of several polygons in one polygon

texts:
I have a dataset of horizontal polygons, which have a vertical distance of 175 meters to each other. These polygons mark certain areas. I would like to have one polygon that covers all the area and doesn't have these blank stripes in it. Basically I want to sort of buffer the polygons, but only vertically, so that they connect to one big polygon.
I am working on ArcGIS 10.3. with a Standard License, but have also QGIS. I have tried to Buffer but I only managed to Buffer in all directions, so that the shape was lost.
Does someone know how to do this?
Here is a sample of the data:


texts: import kml in qgis : kml.kml is not a valid or recognized data source

texts:
When I import kml.kml (generated from ArcGIS) into Google Earth, there is no problem.

However, when I try importing this file as vector layer in QGIS 1.8 on XUbuntu 12.04, I get this error :

kml.kml is not a valid or recognized data source

I need help to resolve this problem.


texts: Areal interpolation with mostly missing data in R

texts:
I am attempting to preform an areal interpolation using R with the sf package's st_interpolate_aw.
Of note, my sampling scheme for my areal units is such that I have many missing values. When using the st_interpolate_aw, I am only returned missing values. Ideally, I am looking for predicted values based on these missing data. See the reprex below:
library(sf)
library(tmap)
url="https://www2.census.gov/geo/tiger/GENZ2018/shp/cb_2018_us_zcta510_500k.zip"
tdir=tempdir()
if(file.exists(paste(tdir,"/cb_2018_us_zcta510_500k.shp",sep=""))==F){
  download.file(url,destfile = file.path(tdir,"ZCTAs.zip"))
  unzip(file.path(tdir,"ZCTAs.zip"),exdir=tdir)}
ZCTAs=read_sf(paste(tdir,"/cb_2018_us_zcta510_500k.shp",sep=""))
df1 = data.frame(val=seq(1:33144))
df1$val[c(501:33144)]=NA
ZCTAs=cbind(ZCTAs,df1)

g = st_make_grid(ZCTAs, n = c(50,50))
test=st_interpolate_aw(ZCTAs["val"], g, extensive = FALSE)
tm_shape(test)+tm_fill(col="val")

Is there a better way to make areal level predictions with many missing values using R?


texts: world city database (with lon,lat) & population per year

texts:
I'm building a visualisation with cities population around the globe for a website, found the MaxMind free database that includes most cities around the globe with coordinates and population, which is awesome.
What I'd like to do is have something like a timeline, to show which cities have increased their population in the last years (even from 2005ish to 2012 would be great).
I looked at the Maxmind site but couldn't find any older versions of their database. Do you know a place where I can find such a database? Thanks!


texts: How to summarize Excel data within the zones of another dataset by using ArcGIS 10.2

texts:
I have the following Excel data, I want to summarize the P-value within the zones of map dataset.

I have sereached a lot and learned I might need to convert Excel data into XY table, and then into a raster before using Zonal Statistics tool.  So my question might eventually be how to convert the Excel data into a raster.

Thanks!

LAT LONG    P-value
42  19      0
-32 -57     0.000199835
-25 -56     0.000148499
-11 -73     4.36478E-06
-11 -72     1.30944E-05
-11 -36     0
-8  -75     0
2   -58     3.41092E-05
.
.
.


texts: Use files with local Danish coordinate system 34 (system 34, S34) in QGIS

texts: The old System 34 is still being used in Denmark even though authorities changed to UTM/ETRS89. We regularly get files, especially when asking for information about local underground cables. Can anyone recommend a free tool/programm which can transform DXF/DWG or shapefiles in S34 to another - in QGIS usable - coordinate system?


texts:
I am trying to pass the geometry of a polyline from a map service layer to the elevation profile widget. Normally, when using a polyline generated from the Draw tool, we can grab the event geometry upon drawing end:

this._toolbar.on("draw-end", (e: any) => this.refreshElevationProfile(e));

refreshElevationProfile(evt): void {
    this._elevationProfileWidget.set("profileGeometry", evt.geometry);
};

Using an identify task, I can get the layers in my map service layer upon map click. But does anyone know how to get specifically the geometry from a feature from one of the layers? (Say the top-most layer for example)


texts: ArcGIS Publisher PMF Version Compatibility

texts:
Does anyone know of version compatibility issues of creating a PMF? Or links to that documentation? I have not seen anything from Esri. I have asked Esri Support, and will update with their response.

Scenario 1:

ArcMap 10.2 with publisher
Microsoft SQL Sever Enterprise Geodatabase
10.1 Geodatabase
Contractors use ArcReader 10.1 64 Bit

Scenario 2:

ArcMap 10.2 with Publisher
9.3 File Geodatabase
Contractors use ArcReader 10.1 64 Bit

Scenario 3:

ArcMap 10.2 with Publisher
10.2 File Geodatabase
Contractors use ArcReader 10.1 64 Bit

I have the above 3 scenarios available. What compatibility issues would I have (or have to consider) when creating a PMF in my ArcMAP 10.2 for users with ArcReader 10.1?

Obviously it would be appropriate for the users to upgrade their ArcReader to the latest version or for whatever version the creators of the PMF are using, but I do not have that control (contract limitations).

Thoughts?


texts: QGIS 1.7.1 and PyQT4 conflict: what am I doing wrong?

texts:
I've just started playing with QGIS with a view to eventually developing a few simple plugins in Python. I downloaded QGIS 1.7.1 for Windows and started working through these workshop notes.
If I've understood correctly, developing plugins often uses QT Designer and something called pyuic4. These don't seem to be included in the QGIS install, so I downloaded PyQT4 separately. This enabled me to finish working through the workshop, but now when I start QGIS I get the following error:
Traceback (most recent call last):
  File "", line 1, in
RuntimeError: the PyQt4.QtCore module is version 1 but the qgis.core
module requires version -1

If I uninstall PyQT4, the problem goes away, but then I don't have QT Designer or pyuic4, so how do I develop plugins? I'm not sure what version -1 of PyQT4 is.
Does anyone have any suggestions, please? I'm obviously doing something daft, but I'm not sure what. Perhaps QT Designer is actually included somewhere with QGIS 1.7.1, but if so I can't find it...
Thanks!


texts: Importing data from a web service as a layer in CartoDB

texts:
I'm trying to figure out if it's possible to bring web services into CartoDB as a layer. In particular, I'm currently trying to import an ArcGIS Feature Service: http://data.nconemap.gov/geoportal/catalog/search/resource/details.page?uuid=%7B20945466-5E79-458C-AB52-5E2899961E45%7D

I'm sure there ought to be some way of doing this, but I can't figure out how...


texts: QGIS Polygon layer table of adjacent polygons

texts: In a polygon layer, how can I generate a table of adjacent polygons?


texts: Convert MapInfo .TAB to Latitude and Longitude Points

texts:
I have been given MapInfo data with .TAB, .MAP, .ID and .DAT files and I am wanting to convert this information into latitude and longitude points.  The goal is to have a set of latitude and longitude points that I can use in an iPad app to create overlays.

I found an application that would load the MapInfo and generate png images suitable for viewing in Google Maps.

Can anyone point me in the right direction?

Update

@Darren Cope helped me find this to work:

ogr2ogr -f "GeoJSON" c:\temp\GeoJsonOutput C:\temp\myMapData.TAB


texts: What file did esriControls.ControlsOpenDocCommand open?

texts:
I'm working on a Windows application, using ESRI's ArcObjects.

In the ArcGIS 10.4 Engine ToolbarControl, I'm placing a esriControls.ControlsOpenDocCommand. With this, the user can select any .mxd file in his local filesystem.

I would like to save the path to this file, so I can restore it, should the user exit the application and the restart.

Does anyone know how I can determine this?


texts: How to create an area cartogram in Shapefile format

texts:
I am working with state level data on health insurance to create a map with each state color-coded based on that data. This is fairly straight forward but some of the largest states geographically are often the least populated. So I want to display my data as an area cartogram map based on population.

Essentially I want to use a map with the states shaped like below as my base shapefile layer.

I found some resources on how to make an area cartogram map but none of them seem to output into a shapefile.


texts: How to calculate a line from a position and bearing in SQL Server 2012 (Spatial)

texts:
Currently I've a table filled with measurements from a certain point. All the rows contain a location in Lat/Lon, a bearing in degrees and some other data containing sensor updates etc. At this moment I'm using the "bearing distance to line tool" in ArcGIS Pro to display the data.
I'm looking for a good method to generate the lines in SQL Server (Spatial). Ultimate goal is to display large amounts of measurements without having to use the tool of ArcGIS Pro.
How can I solve this problem?


texts: Re-referencing raster in rotated coordinate system by "de-rotating" it

texts:
My Trying to define a "rotated regular grid" coordinate system in QGIS question about this got closed because there was already an answer to it, which I understand.  I went to the answer in question and have done my absolute best to implement it, but it was for a slightly different scenario and it is still not working so I seek help in seeing where I went wrong.

I'm not even sure if my pre-processing leading up to implementing the solution is correct.  I am also not completely sure if my attempted solution will even work.  Geodesy was never my strong suit and my last class on it what quite some time ago.  But in order to explain it I need a lot of space as well as diagrams.

The problem is this:

I have downloaded a raster that is defined in a rotated coordinate system with a virtual north pole at 39.25 N, 162.00 W.  I charted this rotation on Google Earth to visualize it:

Now, the solution involves transforming a set of coordinates in the geographic representation (lat, long) to cartesian (x, y, z), performing the rotation, then converting back to geographic representation.  At the end of the answer, a simple "one-liner" set of equations was given as follows:

lat = arcsin(cos(ϑ) sin(lat') - cos(lon') sin(ϑ) cos(lat'))
lon = atan2(sin(lon'), tan(lat') sin(ϑ) + cos(lon') cos(ϑ)) - φ

What I have done is save the image in .asc format, which has the ability of having the xy coordinates of the lower-left corner of the raster defined in the headers of a plain text file.  Since I had to display it in QGIS before converting it, the coordinates in WGS84 are:

xllcorner    -28.342000000000
yllcorner    -23.342000000000

Which is, of course, in the Southern Atlantic Ocean, as can be seen when looking at the raster in WGS84 (remember, the lower-left corner is what we are interested in):

So, I took these coordinates and plugged them into a spreadsheet which implements the "one-liner" equations in the solution that I was directed to:

The equations for calculation lat and lon, respectively are:

=ASIN(COS(B7)*SIN(C13)-COS(C14)*SIN(B7)*COS(C13))
=ATAN2(SIN(C14), TAN(C13)*SIN(B7)+COS(C14)*COS(B7))-B8

These numbers are obviously not correct, but here's an image to drive home how just how not correct it is:

I'm beginning to wonder if the theory behind my implementation of the method is even sound.

Does it make sense to first display the image in WGS84 to grab its "apparent" latlong coordinates in WGS84 and then "de-rotate" them?

I feel like I am in way over my head here but I can't stop thinking about it, I feel like there is some simple piece of the puzzle that is eluding me.

The data that I am working with, as well as a PDF metadata document is here


texts: Variogram model - SAGA on QGIS Processing

texts:
QGIS 2.14 (Win) comes with SAGA 2.12, which Ordinary Kriging supports 6 variogram model options.

However, if I run SAGA Ordinary Kriging from QGIS Processing Toolbox, there is no drop down menu. Only an entry box in which an equation a + bx is shown.

QGIS 2.14 user guide Kriging explains these options, so I can expect them coming. (future...not now).

Do you know any documentation about how to set variogram models in SAGA on QGIS Processing Toolbox? (Ideally, is there any working example?)

At the moment I am playing with Meuse river dataset, so a reference to Spherical model is desirable. (It is not important for now, though.)

So far I tried to manually put equation in the box (like above), which gave me a Kriging-like map. But I have zero confidence in what I am doing.

I installed SAGA 3.0 on QGIS 2.18 Processing, which was the same.


texts: Error importing TIFF file using raster in R: "Failure during raster IO"

texts:
I've been trying to import a TIFF formatted raster image into R using the raster::raster() function with no luck, and I would like to know if the file is corrupted, or if there is something I am missing here (e.g. GDAL or raster setting).
library(raster)
elevation <- raster('raster_image.tif')

A RasterLayer object is created successfully, but when I try to plot it I get this error.
> plot(elevation)
Error in rgdal::getRasterData(con, offset = offs, region.dim = reg, band = object@data@band) :
Failure during raster IO

A link to download the .tif file is here:
https://www.dropbox.com/s/fpxdzpi6ikngr0b/raster_image.tif?dl=0
Also, if I try and access its contents via indexing, I can get output for some values:
> elevation[1:3]
[1] 24.20100 24.25353 24.30607

but other elements I get the same error as for plot():
> elevation[1:300]
Error in rgdal::getRasterData(con, offset = offs, region.dim = reg, band = object@data@band) :
  Failure during raster IO

R Version: 4.0.3
rgdal Version: 1.5-18
GDAL Version: 3.0.4, released 2020/01/28

For context, the TIFF file was downloaded from the publicly available ELVIS elevation data repository in Australia. The actual TIFF is for a box around (151.0642° east, -33.9399° south) and is the 2020 version.
Things I have tried so far
I tried to use GDAL directly via the command line:
gdal_translate raster_image.tif output.tif

But I get the following error related to LZW compression:
Input file size is 2000, 2000
0...10.ERROR 1: LZWDecode:Not enough data at scanline 1024 (short 240 bytes)
ERROR 1: TIFFReadEncodedTile() failed.
ERROR 1: raster_image.tif, band 1: IReadBlock failed at X offset 2, Y offset 1: TIFFReadEncodedTile() failed.

I did manage to render it in QGIS where it doesn't appear to be corrupted:

but when I use the export tool to save to a GeoTIFF the resulting file has a large region missing.

Additionally, there is a warning in the "Rendering" tab of Log Messages (have removed actual file path):
2021-03-09T10:48:34     WARNING    raster_image_39b1418a_1419_4ea7_94ca_2c1dff6d2067 :: file/path/to/raster_image.tif, band 1: IReadBlock failed at X offset 2, Y offset 1: TIFFReadEncodedTile() failed.

I've also tried the gdalwarp solution in this related question:
>gdalwarp raster_image.tif raster_image-fix.tif
Creating output file that is 2000P x 2000L.
Processing raster_image.tif [1/1] : 0Using internal nodata values (e.g. -9999) for image raster_image.tif.
Copying nodata values from source raster_image.tif to destination raster_image-fix.tif.
ERROR 1: LZWDecode:Not enough data at scanline 1024 (short 240 bytes)
ERROR 1: TIFFReadEncodedTile() failed.
ERROR 1: raster_image.tif, band 1: IReadBlock failed at X offset 2, Y offset 1: TIFFReadEncodedTile() failed.

but the resulting file when loaded into R has completely missing values:
> elevation2 <- raster('raster_image-fix.tif')
Warning message:
In showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
  Discarded datum Unknown based on GRS80 ellipsoid in CRS definition
> summary(elevation2)
        raster_image.fix
Min.                  NA
1st Qu.               NA
Median                NA
3rd Qu.               NA
Max.                  NA
NA's             3994240
Warning message:
In .local(object, ...) :
  summary is an estimate based on a sample of 1e+05 cells (2.5% of all cells)


texts: Can feature classes have the same name in two different feature datasets?

texts:
I have an SDE Database called "Vector" which contains many feature datasets eg. Region_1, Region_2 etc. I have to reproject the feature classes inside these datasets.

Instead of creating a new geodatabase, I want to simply create a new feature dataset inside the existing database but call it Region_1_ETRS. All the reprojected feature classes will have their old names, thus making it very easy to reconstruct the map documents.

Unfortunately for me, when I reproject I get an error telling me that the name already exists. Does this mean that I can´t have two feature classes with the same name in one geodatabase even though they are in different feature datasets?


texts: Automatically update records in layer based on changes in another

texts:
I am using QGIS and I have two layers the first one is "Reading" layer with the attribute details

The second layer is "Reading_History" layer with the attribute details

I want to use the expression functions in the default value field to take the last date reading from "Reading_History" layer attribute and update the "Reading" layer with it based on the "ID"` field.

So when a new record for a specific "ID" added to "Reading_History" attribute the "Reading" attribute will be automatically updated to the latest "Reading".
I wonder if there is such a function code to use in the expression calculator that may give the result above.


texts: Failing to create a grid with hexagonal units

texts:
I followed this tutorial step by step.

However, when running the algorithm I get the following error:

  "Horizontal spacing is too small for the covered area"

I'm using EPSG 3857 (WGS 84 Pseudo Mercator) and I'm trying to obtain a 1 HA just as stated in the question I linked. This happens whether I increase or decrease the grid extent or horizontal and vertical spacing. The same error occurs with other projections (4326, 32644, 3763)

Using QGIS 2.18.16


texts: geopandas + shapely. Point within a shape

texts:
I am using geopandas to import a shapefile and then check if a point is within that shape. I am using the shapefile found in http://opendata.atlantaregional.com/datasets/cities-georgia?geometry=-87.766%2C33.036%2C-81.448%2C34.633 to get the shape of Atlanta. However, whenever I "ask" if a point is within that shape (the point is surely within the shape), the function returns False. I have tried creating a Polygon and then asking if a point is within that polygon and it works as expected.

I don't know what else to do to troubleshot or debug this issue.

I am doing all of this using a jupyter notebook. Here is the MWE: https://github.com/vcaballero-salle/mwe-shapefile

The code looks like (full code with imports can be found at the github repo, in the README.md):

plt.clf()
shape = gp.read_file(r'./Cities_Georgia/Cities_Georgia.shp')

atlanta_shape = shape[shape["Name"] == "Atlanta"]

atlanta_shape.plot()

plt.show()
plt.ion()

# In[10]:

def within_shape(df, shapes):
    in_shape = []
    for sh in shapes.geometry:
        within = df.within(sh)
        in_shape.append(within)
    return in_shape

print(within_shape(Point(33.751079, -84.443274), atlanta_shape)) # This returns [False]


texts: testing if a Geodesic polygon contains a point c++

texts:
I am looking for an algorithm, preferrably in C++, that can take a point (lon,lat) and determine whether that point is inside a geodesic polygon.

Is there such an algorithm freely available?


texts: Is there something like SQLselect from MapInfo in QGIS?

texts:
I am a user of MapInfo that will need to adapt increasingly to QGIS.

Is there something like "SQL select" in MapInfo which allows me to easily query between two tables or spatial queries?

MapInfo SQL Window allows you to link easily and efficiently the data in two tables.


texts:
I need to create vector grid 1km x 1km, but I can't figure out what parameters I should set. When I set some value, for example 1000.0, it creates grid with cells with 1000.0 x 1000.0 map degrees instead of meters.

How can I set x - spacing and y - spacing in meters?


texts: How can we merge the two selected features in a vector layer

texts:
I have taken one feature and find which are the neighbouring features to the feature which i have taken and then i have combine all those neighbouring features with my feature to make one feature. i have used the following code for finding neighbouring features.

import itertools
layer = qgis.utils.iface.activeLayer()
for geom1,geom2 in itertools.permutations(layer.getFeatures(),r=2):
          if geom1.geometry().intersects(geom2.geometry()):
              print geom1.attributes(),geom2.attributes()
              layer.select(geom1.id())
              geom = geom1.id.combine( geom2.id)

But when i am using the .combine it is throwing the following errror

***Traceback (most recent call last):
  File "<input>", line 1, in <module>
AttributeError: 'QgsFeature' object has no attribute 'combine'.***

So is there any possible way to combine/merge features programatically by using python programming?


texts: How can I add a custom toolbox to ArcMap, and see its script?

texts: I'd like to add a custom toolbox to ArcMap, and also see the script it's written in. To do the latter, I'm trying to open it with my python editor, but it doesn't open it.


texts: A way to search street intersection on openstreetmap

texts: I need to find intersection by searching name of 2 streets; I found overpass API to query on osm data; the way I found was to find a node near a way; but the problem is all nodes don't have name and in most cases the answer isn't accurate, is there a better way?


texts: How to download products from ESA's scihub using OData protocol?

texts:
How one can download files programatically using scihub OData protocol and cURL (or other library/program)?

You can get a descriptive xml response of the API for a  dataset can be obtained by this example:

curl -u username:password "https://scihub.esa.int/dhus/odata/v1/Products('177db57e-56ca-4588-830a-cfb8712ecf6f')/$value/"

But how you can actually use it to pull actual data ?


texts: QGIS Attribute Table on print layout not sorting

texts:
I've added in an attribute table to a print layout and am trying to sort numerically ascending using Plot Order

For some reason it doesn't change order whatever field I use to sort the items in the table.
It seems to be stuck sorting on Feature_ID which I can't edit.
I'm using QGIS version 3.14.16-Pi
Can someone direct me where I might be going wrong?


texts: Distance in Meters/Kilometers from pgr_drivingdistance

texts:
I am having trouble calculating distance (driving distance) along a route with the pgr_drivingdistance function. I have a straightforward workflow, I download OSM data, and put that into postgresql with osm2pgrouting, and am trying to calculate a route between two nodes. I just want the distance for that route. This is my query:
select sum(cost) from pgr_drivingdistance('select gid as id, source, target, cost from ways', 84216, 84216, true);

I know cost is not the distance, but when I try to pass length_m as cost, I don't get the correct answer.


texts:
Hello I was trying to but together a quick and dirty bing/osm  map using leaflet and leaflet plugin. My problem comes in when adding a layer from cartodb.The points do not show and I am not getting any errors, however when I change opacity on the tile where my points should be they show up. So it looks like the base maps are taking priority over my layer. Is there a way to fix this?

    function main() {
    var map = new L.Map('map', {
            center: new L.LatLng(39.095,-75.54),
            zoom: 2
    });

    var osm = new L.TileLayer('http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png');
    var bing = new L.BingLayer("AulNwbEkrE4P0VQrC8JI3SKoJDC4A9V63I8q4tixTNo6HUkfUyjWJuXSgpgysdZ5");
        map.addLayer(bing);bing.bringToBack();
        map.addControl(new L.Control.Layers({'OSM':osm, "Bing":bing}, {}));
          cartodb.createLayer(map, 'http://centuryeng.cartodb.com/api/v2/viz/a7fd6eb2-f6ee-11e3-b807-0e73339ffa50/viz.json')
            .addTo(map)
            .on('done', function (layer) {

            });
    console.log(map);
    }

    window.onload = main;