texts: Virtual layer functionality does not seem to support memory layers

texts:
I have a memory layer loaded in the QGIS virtual layer in the data source manager. I am trying to run a very simple query:

select * from layer

Here is a snapshot:

While this works for the physical layer, it does not work for the memory layer. They output virtual layer is blank, and I know that my memory layer has attribute data.

Is this simply a bug? I want to keep this a memory layer to reduce the number of steps I need to do to arrive at my final result. Should I resort to another tool or is there a workaround?

EDIT

Also found this quote in this link:

  Creating a virtual layer with the layer_ref key allows to directly access already loaded QGIS layers (including memory layers).


texts:
i'm expected output picture attached to him by sample google map.

given the following my leaflet map script. but cannot enable to achieve i will expected output.

function single_vehicle_route(data){

            var latlng_val =[];

            var defaultCoords = [11.016844,76.955832];
                var StartMarker = {};
                var EndMarker = {};
                //set up our map
                var map = L.map('single_map').setView(defaultCoords,6);
                L.tileLayer('http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png',
                    {
                    maxZoom: 16
                    }).addTo(map);

            for (i = 0; i < data.length; i++) {

                var latlngs = [[data[i].lat,data[i].lng]];

                            var polyline = L.polyline(latlngs, {color: 'red'}).addTo(map);

                            map.fitBounds(polyline.getBounds());

            }
            //console.log(latlng_val);

    }


texts: Modifying GDAL rasterize tool using Python

texts:
I would like to take a look and maybe modify the code of the gdal tool rasterize and I don't find the python file of the tool.

Where can I find the gdal python files so that I can open it and make some changes?


texts: When I query for larger features I get esri.config.defaults.io.proxyUrl is not set error. Can someone please explain why proxyurl is needed, how I should set it up, what should it point to, and any other relevant information about it? I am running Websphere and IBM http server.


texts: Participative maps for collecting geographic information?

texts:
I'm looking for a simple way to collect information for a participative plans (regarding landscape) and show them on Google Maps.

I'm thinking to a small web form for asking information to public (es: point of interest, ), then they have to digtize the refered point, line or polygon. After submit, the information are showed in a map.

Any ideas? Like using Joomla extension or similar.


texts: Exporting attachments gives RuntimeError?

texts:
I found this code on a youtube tutorial (https://www.youtube.com/watch?v=u26Kpe4NBqU) for batching exported attachments from arcgis online.

from arcpy import da
import os

attachTable = arcpy.GetParameterAsText(0)
origTable = arcpy.GetParameterAsText(1)
nameField = arcpy.GetParameterAsText(2)
fileLocation= arcpy.GetParameterAsText(3)

attachCursor = da.SearchCursor(attachTable,['DATA', 'ATT_NAME','REL_GLOBALID'])

for attRow in attachCursor:
    binaryRep = attRow[0]
    fileName = attRow[1]
    relID = attRow[2]

    originCursor = da.SearchCursor(origTable, ['GLOBALID', nameField])

    for origRow in originCursor:
        origID = origRow[0]
        origName = origRow[1]
        if origID == relID:
             break

    open(fileLocation + os.sep + origName + ".jpg",'wb').write(binaryRep.tobytes())

    del originCursor

BUT when I enter the script into python, I get an error here:

attachCursor = da.SearchCursor(attachTable,['DATA','ATT_NAME','REL_GLOBALID'])
Runtime error
Traceback (most recent call last):
  File "<string>", line 1, in <module>
RuntimeError: cannot open ''*

Does anyone know why I receive this error?

I think it has something to do with the Cursor (either attach or search), but I'm stuck and can only scratch my head in confusion.

...
I've also tried to enter it as a script, when I do it will run but usually never to completion. It doesn't crash but just keeps going. The one time it worked took an hour for a handful of attachments and though it says it was successful, I couldn't find the exported files. I don't have enough reputation points to add all the steps I took as a pictures so here is a link to a pdf with screenshots:

link to steps taken to create and run script


texts: Problem with importing shapefile to PostGIS table

texts:
I have problem with importing shapefile to Postgis table.

PostGIS table has:

the_geom geometry(LineString,102067) for geometry

Shapefile has:

INFO: Open of `postup_doby.shp'
using driver `ESRI Shapefile' successful.
1: postup_doby (Line String)

ogr2ogr command is:

PGCLIENTENCODING=WIN1250 ogr2ogr -skipfailures --config PG_USE_COPY YES -f PostgreSQL PG:"host=localhost port=5432 dbname=myname password=secure active_schema=povis user=sa" -lco DIM=2 postup_doby.shp postup_doby -append -nlt GEOMETRYCOLLECTION -lco GEOMETRY_NAME=the_geom -lco FID=envi_id -nln postupove_doby PRECISION=NO

Error message is:

Insertion is likely to fail
ERROR 1: COPY statement failed.
ERROR:  Geometry type (MultiLineString) does not match column type (LineString)
CONTEXT:  COPY postupove_doby, line 1, column the_geom: "0105000020B38E01000200000001020000006C000000E4E2809372B61CC13F49396A56C030C1C24C8D97FAB41CC175C2C6B5..."

Problem is in shapefile, because has 10 features, which one feature has MULTILINESTRING geometry, others has LINESTRING geometry.

ogrinfo of feature with MULTILINESTRING, others has LINESTRING:

MULTILINESTRING ((-470428.644046349683776 -1097814.414936616783962,-470334.64799995359499 -1097731.710 ...

My solution is convert feature with MULTILINESTRING to LINESTRING and import to PostGIS, but I dont know, how I make it in comand line.

How can I resove this problem?


texts: Installing ECW for QGIS in Ubuntu 12.04?

texts:
I have followed various online help guides on how to make ecw for QGIS 1.8 work in Ubuntu 12.04. I had it all working in Ubuntu 10.04 but since migration to 12.04 no such luck.

Is there a guide that works for Ubuntu 12.04?

At the moment it seems easier to install the windows version where ecw support is included in the install, but I would like to stay with Ubuntu.


texts: Export SQL Server to Shapefile with OGR, difference between 2014 and 2012

texts:
I'm trying to export some data from a SQL Server table to a shapefile.

I have tested this on a 2014 SQL Sever and it works but when I try and run it on a 2012 SQL Server it creates a blank output.

The original that works SQL Server, 2014. Result is as expected with a shapefile created:

ogr2ogr -overwrite -f "ESRI Shapefile" "C:\temp\sqlexport_Works.shp" "MSSQL:server=[SERVER];database=[DB];trusted_connection=yes;GeometryFormat=native;" -sql "SELECT TOP 1000 [ID] ,GEOMETRY::STGeomFromWKB([geom].STAsBinary(),3857) As [geom] FROM [DB].[dbo].[RoadNetwork] WHERE [geom].STGeometryType() = 'MULTILINESTRING'" -lco "SHPT=ARC" -mapFieldType Binary=String

The alternative version, that does not work on SQL Server, 2012. Result is an empty shapefile with only a populated dbf file created:

ogr2ogr -overwrite -f "ESRI Shapefile" "C:\temp\sqlexport_DoesNotWork.shp" "MSSQL:server=[SERVER];database=[DB];trusted_connection=yes;GeometryFormat=native" -sql "SELECT TOP 1000 [ID] ,GEOMETRY::STGeomFromWKB([geom].STAsBinary(),3857) As [geom] FROM [DB].[dbo].[RoadNetwork] WHERE [geom].STGeometryType() = 'MULTILINESTRING'" -lco "SHPT=ARC" -mapFieldType Binary=String

Can anybody explain how to get this working on SQL Server 2012?


texts:
I have a relatively large worldwide raster mbTiles file from which I want to convert individual countries into new mbtiles. Unfortunately i fail to use tilelive which is actually predestined for this purpose.
Has anyone ever run tilelive on an Ubuntu server and got the tilelive-copy command implemented? I only get the error message.
i am going to do the following:

installation of node js
installation of npm
npm install tilelive
npm install @mapbox/mbtiles
npm install mbtiles

After executing the command: tilelive-copy input.mbtiles ouput.mbtiles I only get the error message Error: Invalid tilesource protocol: mbtiles.
I have now tried different nodejs versions etc. Nothing leads to success.
Does anybody have an idea what the problem is or an alternative approach e.g. with GDAL?
Thanks!


texts: Sampling from all the images in a collection and combining them into one FeatureCollection

texts:
I'm trying to collect a set of sample points from all the images in a collection. However instead of outputting a Feature Collection that consists of a group of Features with the training data as properties, I get a Feature Collection that contains another nested Feature Collection for each image I process, and cannot tell if these contain the training data I want.

How can I output the data as one featureCollection containing nothing but features? Do I have to use the .iterate()? From the examples I am unsure how to apply this for my purpose.

    // Create Example Collection
var collection =  ee.ImageCollection("LANDSAT/LT05/C01/T1_TOA")
  .filterDate('2001-01-01', '2001-12-31')
  .filterBounds(ee.Geometry.Point(-92.94,37.52))
  .filterMetadata('CLOUD_COVER', 'equals', 0);

print(collection);

// function for sampling training data
function sampleImage(image) {
    var bands = ['B1', 'B2', 'B3'];
    var param = {region: clust, scale: 30, numPixels: 500};
    var SamPix = image.select(bands).sample(param);
    return SamPix;
}
var training = collection.map(sampleImage);
print(training, "sample points");


texts: How to create a Button with custom Icons in Mapinfo 10?

texts:
I am using mapinfo 10 and I have created dll with Pelles C with 16*16 and 24*24 bmp images. I placed the dll file into folder containing .mbx file but I get error :

  Unableable to locate external icon resource.

Am I missing something.

Code:

Create ButtonPad "mybuttonpad"  As
        ToggleButton
        ID 1001
        Icon 8001 File "icon.dll"
          Calling vendor_ButtonAction
          HelpMsg "ICONTEST"
Show


texts: QGIS 2.8 problem accessing GRASS commands

texts:
QGIS 2.8.1 Win7 64 bit installed with osgeo4w-setup.

I’m trying to run GRASS commands on a LIDAR raster layer but I get the following message:

  Missing dependency.... GRASS not installed in your system` message each time.

GRASS appears in the processing tools pane and the “Processing options” for GRASS shows the Msys and GRASS folders locations correctly (which was a problem for a similar questioner). I note that GRASS has been 7.0 stable since mid February but osgeo4w-setup only shows GRASS 6. All manuals refer to a ‘GRASS plugin’ but none seem to be listed now.

Is there 'best-practice' way of installing GRASS functionality in QGIS 2.8 on Windows7 64 bit?

As a supplementary question is it possible to install a GRASS add-on that is not part of the core GRASS setup (I need r.skyview for LIDAR visualisation of archaeological features)?


texts: Looking for a GeodesicDistance() technique for getting accurate distances using NetTopologySuite for points initially in 3857

texts:
I am looking for a technique for getting accurate distances between two points in the 3857 projection.  I believe solutions exist which transform the coordinates to 4326 and compute the distance from there.  I believe this is called the "geodesic distance."  I cannot find how the same is done with NetTopologySuite.
Some thoughts:

I can transform the points to 4326 and compute a distance (maintenance concerns)
Perhaps there is a library which selects a suitable local projection given a 3857 coordinate?  If so I could project to that SRS before calculating the distance (maintenance concerns).

There is code readily available online for computing a geodesic distance:
private double distance(double lat1, double lon1, double lat2, double lon2, char unit) {
  if ((lat1 == lat2) && (lon1 == lon2)) {
    return 0;
  }
  else {
    double theta = lon1 - lon2;
    double dist = Math.Sin(deg2rad(lat1)) * Math.Sin(deg2rad(lat2)) + Math.Cos(deg2rad(lat1)) * Math.Cos(deg2rad(lat2)) * Math.Cos(deg2rad(theta));
    dist = Math.Acos(dist);
    dist = rad2deg(dist);
    dist = dist * 60 * 1.1515;
    if (unit == 'K') {
      dist = dist * 1.609344;
    } else if (unit == 'N') {
      dist = dist * 0.8684;
    }
    return (dist);
  }
}

But I don't want to write code that might already be available and it seems like NetTopologySuite would already have this as a built-in solution.


texts: Select points with maximum attribute value per category on a spatial join with PostGIS

texts:
I'm trying to select the points with the maximum altitude value contained in the polygons of a grid. grid and point with altitude

I've managed to get the max value per category:

SELECT grid.pagename, max(pnt.altitude) AS max_altitude
FROM grid_150m as grid LEFT JOIN pnt_cote_index_proj as pnt
ON st_contains(grid.geom, pnt.geom)
GROUP BY grid.pagename

The hard part is to get the actual geometry, without success I've tried:

SELECT distinct on (pnt.altitude) pnt.altitude, grid.pagename, pnt.geom
FROM grid_150m as grid
LEFT JOIN pnt_cote_index_proj as pnt
ON st_contains(grid.geom, pnt.geom)
group by grid.pagename, pnt.altitude, pnt.geom
order by pnt.altitude;


texts:
Why do I keep receiving this error message? I am using clip raster by extent as a batch process and am receiving this error each time. Any clue?

              File "C:/PROGRA~1/QGIS3~1.6/apps/qgis/./python/plugins\processing\gui\BatchAlgorithmDialog.py", line 88, in runAlgorithm
              parameters[param.name()] = wrapper.parameterValue()
             AttributeError: 'NoneType' object has no attribute 'parameterValue'```


texts:
I have already built a desktop application for the management of archaeological settlements. Currently the data are stored in shapefiles (spatial data) and MYSQL (descriptive information). The application is based on MapWindow and MapWindow plugins. In the database there are some relationships between tables (one-to-many, many-to-many).

I need to convert my application to a web mapping application. I intend to share my spatial data with GeoServer.  The web application should have forms to search the descriptive data of mysql and tools to select and zoom to the corresponding spatial data. Almost all the existed web mapping solutions have query capabilites but only on the spatial data (i.e. fields of a shapefile). Unfortunately my data are stored separately in mysql.

What's the best solution to implement my web mapping application?
Do I have to make custom tools?
I intend to use Geomoose and some special plugins or Mapbender.
Are they the appropriate platforms?

I need a platform to join the spatial data and my mysql data in a easy and usefull way. I try to avoid extra programming for the forms and the connection between spatial and descriptive data.


texts: Plot points using ArcGIS javascript api from co-ordinates

texts:
I would like to map co-ordinates that I got from the exif in photo's using my own command line tool that I made. As opposed to uploading your private photo's for everyone to see on flickr.

How to I get the coordinates plotted on a map like this:
http://help.arcgis.com/en/webapi/javascript/arcgis/jssamples/#sample/map_graphics

Using a file generated from the stdout of my command line tool, here is an example:
http://dpaste.com/880696/plain/

I can't find any info on how esri.request works in the Javascript API for arcGIS. I also tried looking at this loading text sample in arcgis javascript api's sample section, the example doesn't seem to work.

Basically all I'm asking is what is the simplest method to plot co-ordinates on the map as they are doing in the example with flickr, but from a text file. Perhaps I should first turn them into the format that is being used in the example with flickr? I tried using this format: on georss dot org but using sed to replace the values with my own values but it didn't work.


texts:
The documentation for Map.querySourceFeatures says it takes a parameter for sourceID - the ID of the vector tile or GeoJSON source. In their example, they add a source manually:

map.addSource('counties', {
    "type": "vector",
    "url": "mapbox://mapbox.82pkq93d"
});

And so the sourceID is counties.

If I create a Mapbox style in Mapbox Studio and publish it with a url, say:

mapbox://styles/mycoolcompany/mycoolmap123

and pass that directly to the Map constructor like this:

this.map = new mapboxgl.Map({
    container: 'map',
    style: 'mapbox://styles/mycoolcompany/mycoolmap123',
    center: [-123.12, 49.28],
    zoom: 5,
});

How can I get (or set?) the styleID programatically? I'm trying different combinations of things from Mapbox Studio but querySourceFeatures returns an empty array every time with no exceptions or indication of what it's missing.

More succinctly: how do I use querySourceFeatures on a home-made Mapbox Studio Style?

(mapbox-gl 0.36.0)


texts: OSGeo4W QGIS 3.8: don't understand setting of internal Path variable

texts:
As I understand it, when QGIS is started from qgis.bat some OSGeo4w specific folder paths are prepended to the OS PATH variable. When I do this and go look in Settings>Options>System>Current environment variables>Path I see

C:\OSGEO4~1\apps\qgis\bin;C:\OSGEO4~1\apps\Python37;C:\OSGEO4~1\apps\Python37\Scripts;C:\OSGEO4~1\apps\qt5\bin;{app};C:\OSGEO4~1\apps\Python27\Scripts;C:\OSGEO4~1\bin;C:\WINDOWS\system32;C:\WINDOWS;C:\WINDOWS\system32\WBem;C:\Program Files\R\R-3.6.1\bin\x64;C:\OSGEO4~1\apps\Python37\lib\site-packages\pywin32_system32;C:\OSGEO4~1\apps\Python37\lib\site-packages\numpy\.libs

If I re-start QGIS using qgis-bin.exe instead and go look at the Path variable, it's the same value. I would have thought that the box would be empty and that the Python console would be interrogating the OS Path-listed instance of Python instead, which on my computer is, in fact, under an ArcGIS Pro folder (this isn't my goal, just part of why I don't understand what qgis.bat is supposed to be doing).

In fact, when I edit the batch file by commenting out the call to py_env.bat, eg.

rem call py3_env.bat

and start again, the Path value is still the same. Why is that?


texts: Why does QGIS crash when extracting pixel values by GDAL?

texts:
In below code I want to print pixel values, but when run in QGIS python console QGIS software crashes and closed! what is the problem?

from osgeo import gdal
import numpy
import struct

def pixel2coord(x, y):
    xp = (pixelWidth * x) + originX + (pixelWidth/2)
    yp = (pixelHeight * y) + originY + (pixelHeight /2)
    return QgsPoint(xp, yp)

# Open tif file
ds = QgsRasterLayer("/Data/Sample_Polygon.tif","Study")

pixelWidth = ds.rasterUnitsPerPixelX()
pixelHeight = ds.rasterUnitsPerPixelY()

# extent of the Raster layer
ext = ds.extent()

originX ,originY = (ext.xMinimum(),ext.yMinimum())
src_cols = ds.width()
src_rows = ds.height()
drive = gdal.GetDriverByName( 'GTiff' )
src_ds = gdal.Open("/Data/Sample_Polygon.tif")
outBand = src_ds.GetRasterBand(1)

outData = numpy.zeros((src_cols, src_rows), numpy.float32)

pntRstList = []

for i in range(0, src_cols):
    for j in range(0, src_rows):
        rspnt = pixel2coord(i,j)
        pntRstList.append(rspnt)

for x, y in pntRstList:
    structval=outBand.ReadRaster(x,y,1,1,buf_type=gdal.GDT_Float32)
    intval = struct.unpack('f' , structval)
    val=intval[0]
    print val

Also I test the @xunilk code in this link, and closed the QGIS software again. I use QGIS 2.18.9 in Linux Mint 18.


texts: Including variable in where clause of arcpy.Select_analysis()?

texts:
I am trying to loop through a shapefile, selecting each feature in turn and copying it to a temporary shapefile to by included in a union analysis. I'm using a cursor to find the ID name for each feature which I'm setting to a varible 'Name'. Whenever I try to use this variable as part of the where clause in arcpy.Select_analysis I get an error:

ExecuteError: ERROR 999999: Error executing function.
An invalid SQL statement was used.
An invalid SQL statement was used.
Failed to execute (Select).

The code I'm using is:

Name = 101
where = "\'\"StudyID\" = \\'"+str(Name)+"\\'\'"
arcpy.Select_analysis("C:\\input.shp", "C:\\output.shp", where)

If I type it out without using the variables:

arcpy.Select_analysis("C:\\input.shp", "C:\\output.shp", '"StudyID" = \'101\'')

it works fine

What to I need to do to fit the variable into the sql statement?


texts:
I'm stuggling to add a KML layer to a simple test OpenLayers map - no errors are shown, but the KML layer does not display on the map.

I've tried multiple samples including the Sundials on a Spherical Mercator sample and the Layers help page. I can open this example in Firefox and see the KML points on the map without error.

I copied this sample's source code and adjusted the paths to the CSS, JS and KML files to use full paths.

However, when opening this map the KML points are not displayed. Firebug shows GET sundials.kml on the working example from OpenLayers, but OPTIONS sundials.kml on my version.

Any advice on what I'm doing wrong? Thanks

<!DOCTYPE html>

    <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
    <meta name="viewport" content="width=device-width, initial-scale=1.0, maximum-scale=1.0, user-scalable=0">
    <meta name="apple-mobile-web-app-capable" content="yes">
    <link rel="stylesheet" href="http://openlayers.org/dev/theme/default/style.css" type="text/css">
    <link rel="stylesheet" href="http://openlayers.org/dev/examples/style.css" type="text/css">

KML Layer Example

          kml, popup, feature

          Demonstrates loading and displaying a KML file on top of a basemap.


texts: Join layers in PyQGIS, on multiple fields

texts:
Is it possible to join two layers in PyQGIS, on multiple fields?

ie, JOIN ON layer1field1 = layer2field1 AND layer1field2 = layer2field2

If so, how?


texts: How can I add features to a GRID raster dataset?

texts:
I am using a high resolution GRID, but unfortunatelly have realized I have portions of no data within it. As a solution I have decided to create and add the elevations myself for some of the areas. In order to do this I have created a separate feature class containing the additional features with coresponding elevations.

Is there a way I can incorporate the feature class into the GRID or do I have to convert it into a TIN?

( i know I can then use the add feature to tin TOOL). dUEto the high resolution, TIN will not work out, process keeps on giving me errors. (I am using ArcMap 10.1)


texts:
i'm working with a hudge amount of data, that's why i decided to use clustermarker on my leaflet map :

function CickOnPoint (e) {
var pm = e.target.feature.properties.distCoupe;//recupère dist coupe depuis geojson
var ind = findPosSlideIndex(pm);
currentSlide(ind);
}

function onEachPoint (feature,layer){
layer.on({click: CickOnPoint});
//marker pour cluster
var coordinates = feature.geometry.coordinates;
var latlng = new L.LatLng(coordinates[1], coordinates[0]);
var marker = L.marker(latlng);
marker.bindPopup("SURPRISE MOTHERFUCKER");
return marker;
}

var geo_point = L.geoJSON(point_layer_WGS84_dist, {
onEachFeature: onEachPoint
});

var clusters = new L.MarkerClusterGroup({
    chunkedLoading: true,
    disableClusteringAtZoom: 19,
    spiderfyOnMaxZoom: false
});
clusters.addLayer(geo_point);
map.addLayer(clusters);

I can't understand why no popup is showing when i click on markers, but the fonction "clickOnPoint", which control a picture viewer does work.

I'm maybe missing a comprehension point in here.


texts:
Here is my code and (console.log(feature)) gives me URL that has data in console. How can I fetch data and use it as var?
map.on('singleclick', function(evt) {
    var view = map.getView();
    var viewResolution = view.getResolution();
    var feature = UnTiled.getSource().getFeatureInfoUrl(
        evt['coordinate'],
        viewResolution,
        'EPSG:3857', 'EPSG:4326',
        {'INFO_FORMAT': 'text/html'}
      );
      console.log(feature);

     });

Here is one of the URLs:
http://localhost:8080/geoserver/shomal/wms?0=E&1=P&2=S&3=G&4=%3A&5=4&6=3&7=2&8=6&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetFeatureInfo&FORMAT=image%2Fpng&TRANSPARENT=true&QUERY_LAYERS=shomal%3Ageotiff_coverage&STYLES=&LAYERS=shomal%3Ageotiff_coverage&exceptions=application%2Fvnd.ogc.se_inimage&X=50&Y=50&SRS=EPSG%3A3857&WIDTH=101&HEIGHT=101&BBOX=5518113.346630408%2C4109935.178502069%2C6012202.297465787%2C4604024.129337448


texts: Can you edit an SDE table via ArcObjects when not added to ArcMap

texts:
I have a table that is in SDE, what I want to do is in my userform click on the "Save" button, and the table will be edited even though it is not within the map document. Is this possible?
I have tried the following code, which runs, but nothing happens to the table (I have also added the table to the map document and still nothing happens or no error messages occur).
Here is that I have code.

Private Sub cmdSaveDev_Click()

Dim pTable As ITable
Dim pCursor As ICursor
Dim pRow As IRow
Dim pRowBuffer As IRowBuffer

Dim pPropset As IPropertySet
Set pPropset = New PropertySet
    With pPropset
        .SetProperty "Server", "***"
        .SetProperty "Instance", "***"
        .SetProperty "User", "***"
        .SetProperty "Password", "***"
        .SetProperty "Version", "SDE.DEFAULT"
    End With

Dim pSdeWS As IWorkspace

Dim pSdeWF As IWorkspaceFactory
Set pSdeWF = New SdeWorkspaceFactory

Dim pSdeFeatWS As IFeatureWorkspace
Set pSdeFeatWS = pSdeWF.Open(pPropset, 0)

Set pSdeWS = pSdeFeatWS

Dim pWkspEdit As IWorkspaceEdit
Set pWkspEdit = pSdeFeatWS

Set pTable = pSdeFeatWS.OpenTable("SY_CSY_TERRIER.DEVELOPMENT_DEVELOPER")

Dim i As Integer

If lstSelectedDevelopers.ListCount = 0 Then
    MsgBox "Add a developer to the Selected Developer box and try again"
Else

    pWkspEdit.StartEditing (True)
    pWkspEdit.StartEditOperation

    Set pRowBuffer = pTable.CreateRowBuffer
    Set pCursor = pTable.Insert(True)
    For i = 0 To lstSelectedDevelopers.ListCount - 1
        Set pRow = pRowBuffer
        pRow.value(1) = frmCityDev.txtRecordNo
        pRow.value(2) = lstSelectedDevelopers.List(i)

        pCursor.Flush
    Next i

    pWkspEdit.StopEditing (True)
    pWkspEdit.StopEditOperation

End Sub

Thanks!


texts: How do you embed iFrame content into a CartoDB infoWindow?

texts:
All,

I'm trying to embed Youtube videos and SoundCloud audio clips into a CartoDB infoWindow. I'm trying to do this by simply using the {{embed}} column and entering it into the custom infoWindow html...

  x

name

{{name}}

description

{{description}}

                {{embed}}

To see the result, click one of the two point features in the map:

Result Map:

http://cdb.io/K8yZ9V

SoundCloud embed code:

Youtube embed code:


texts: Problem with PostGIS in Opengeo Suite in Ubuntu 13?

texts: Everything is working fine after installing the Opengeo Suite in Ubuntu except PostGIS issues. In the dashboard the PostGIS manage is disabled(grey) but PostgreSQL is working via pgadmin3. And also when we create a new database the tables "Geometry Columns" & "Spatial Ref_System" isn't getting created. Also we are not able to import any shape files.


texts:
I have a simple OpenStreetMap for which I want to display a series of circles. Each circle represents readership for that area, so its radius will change based on the numbers.

However, when I create the circle and zoom in the map, the circle zooms in as well - it seems like OpenLayers assumes the circle is "glued" to the map and represents a piece of it, when really it's more like a programmatically created marker.

How can I get the circle to stay the same pixel radius no matter what the zoom level of the map? Or is there some way to get a marker to be created as a Vector instead of an image?

var map = new OpenLayers.Map({
  div: 'geomap',
  projection: "EPSG:900913"
});
map.addLayer(new OpenLayers.Layer.OSM());

var epsg4326 = new OpenLayers.Projection("EPSG:4326"); // WGS 1984 projection
var projectTo = map.getProjectionObject(); // The map projection (Spherical Mercator)
// test point
var lonLat = new OpenLayers.LonLat(-74.265982, 40.665443).
    transform(epsg4326, projectTo);

var vectorLayer = new OpenLayers.Layer.Vector("Overlay");
var point = new OpenLayers.Geometry.Point(lonLat.lon, lonLat.lat);

var mycircle = OpenLayers.Geometry.Polygon.createRegularPolygon
    (
        point,
        50000,
        40,
        0
    );

var featurecircle = new OpenLayers.Feature.Vector(mycircle);
vectorLayer.addFeatures([featurecircle]);

map.addLayer(vectorLayer);
map.zoomToExtent(vectorLayer.getDataExtent());

UPDATE: Here is the final answer (thanks Rene)!

var map = new OpenLayers.Map({
    div: 'geomap',
    projection: "EPSG:900913"
  });
  map.addLayer(new OpenLayers.Layer.OSM());

  var epsg4326 = new OpenLayers.Projection("EPSG:4326"); // WGS 1984 projection
  var projectTo = map.getProjectionObject(); // The map projection (Spherical Mercator)
  var lonLat = new OpenLayers.LonLat(-74.265982, 40.665443).
      transform(epsg4326, projectTo);

  var vectorLayer = new OpenLayers.Layer.Vector('FSA', {
    styleMap: new OpenLayers.StyleMap({
      "default": new OpenLayers.Style({
        pointRadius: "${radius}",
        fillColor: '#ff0000',
        fillOpacity: 1.0,
        strokeColor: '#000000',
        strokeWidth: 2,
        strokeOpacity: 1.0
      })
    })
  });
  map.addLayer(vectorLayer);

  var feature = new OpenLayers.Feature.Vector(
      new OpenLayers.Geometry.Point(lonLat.lon, lonLat.lat),
      {
        radius: 10
      }
  );
  vectorLayer.addFeatures([feature]);

  map.zoomToExtent(vectorLayer.getDataExtent());


texts: Can ArcGIS Online navigate and identify large GIS datasets?

texts:
For my upcoming project, I want to see if ArcGIS Online can provide a viable solution.

Users will need to have access to webmap on their desktops with basic navigation and identify. This will be an internal application.

As far as I know, hosting data on the web is not an issue from security point of view.

I want to know if ArcGIS Online can handle large GIS datasets, it will range from 50,000-250,000 records for some of the basemap data layers.


texts: Limit construction tools listed for an edit template to custom edit tool

texts:
Is it possible to limit the construction tools available for a specific feature class to a custom edit tool?  I'm porting a 9.3 extension to 10, and would like to move tools off an old toolbar and into the Construction Tools panel.  I've ported one tool over.  However, I want the tool to only show when the feature class it is inteded to work with is set as the current edit template.  In addition, I want to hide the standard construction tools that are available.

The feature class has a class extension associated with it.  I've been able to use the ClassID within the IEditEvents5.OnCurrentTemplateChanged event to determine if the current edit template is the correct feature class.  But all I can do with the tool is enable it or disable it.  I want to be able to hide the custom tool all together for all other feature classes.  In addition, I want to hide the standard tools that are displayed as well, limiting the users options to just the custom edit tool.

If I can't make this work I'll just go back to placing the tools on a toolbar again, but I would like to be able to intergrate the tools into the native editing environment if possible.


texts: SLD function/transformation based on table values

texts:
I'm trying to accomplish an apparently simple task, to display some lines from a postgis table, but I don't know how to do it.
The table has, say, 3 fields, one for the geometry and the other two have some value.
I'd like to set the color of the line based on one of the fields, I'll call it field1, and to set the dash pattern of the line based on the second field (field2).

field1
value  color
1      black
2      red
3      green
...

field2
value  dash pattern
1      solid  (e.g. values for stroke-dasharray)
2      dot line
...

I  know how to create a style to accomplish that but mixing the filters/rules for managing 2 fields is a mess.
Using another language I'd use a function with 2 parameters (the values of field1 and field2) returning the values to use for stroke and stroke-dasharray.
I searched the docs and the only thing I found out is the function Recode but, as I understand, it uses only one input parameter. I read the docs about Transformations but I assume that they must be created outside SLD, I think it'd be overkill for my needs.
Is there a way to accomplish my task?


texts: Resolution for different layers using Python GDAL

texts:
I'm extracting data from a GRIB file using GDAL in Python 3 on CentOS.

The model is 40 km resolution, and is usually 151x93.

When I try to extract data from band 195, I get the following error:

  Error 4: 'vsimem/work_grib_0x2a194da.jpc' not recognized as a
  supported file format.

  dec_jpeg2000: Unable to open JPEG2000 image within GRIB file.

Is the JPEG2000 driver available?

  Warning 1: Band 295 of GRIB dataset is 108x94, while the first band
  and dataset is 151x113. Georeferencing of band 295 may be incorrect,
  and data access may be incomplete.

Could the resolution warning be coming from using the wrong driver?

When I call band.XSize, it reports 151, rather than 108.

Can this be resolved in Python?

Edit: the source for the warning message is found here:

http://even.rouault.free.fr/coverage/frmts/grib/gribdataset.cpp.gcov.html


texts: Editing symbology of contour lines with ArcPy

texts:
I just have a problem with my Python script.  I had run the Contour-tool to generate polylines.

Contour("rectExtract", "C:/fakepath/Class10.shp", 20,0)

Now I want to mark the 100 meter lines in another colour.


texts: Cannot get 2 raster files of same cell size to align ArcMap 10.2

texts:
I have two raster files that are the same cell size but I cannot get them to align properly. (I want the black to line up with the pink)

When creating the black raster file (from using polyline to raster), I imported the pink raster in the Snap Raster section in the Environments setting. I thought that Snap Raster was supposed to make it line up?

Am I doing something wrong or is it a bug in Arc (I remember reading that this was a problem about two years ago), and if it is a bug is there a work around?

I have also read there is a Raster Shift, but I'm not sure what to enter in the X Coordinate and Y Coordinate. Since I only want the black raster to shift down, I know that the X Coor = 0 and the Y Coor should be negative something, but what are the units?

Thanks.


texts: How to generate legend in leaflet webmap?

texts:
Although I have checked the "Create Legend" option in the qgis2leaf plugin, it did not generate any legend in the output webmap.

What am I missing here ?


texts: Publishing layer (shapefile or postgis) with TileMill moves it 200 meters?

texts:
I am working with the TileMill application to publish layers.

Well, first I have posted directly loading the shp file., But when the charge goes I moved about 200 meters.

Then as I have the option, I imported the layer from postgis and my surprise was that the movement was more than 1000 km.

After choosing SRS: Custom and put the appropriate reference system, the movement returned to be 200 meters.

Does anyone have any idea why this might happen?


texts: Fix Black Border in Python Batch Raster Conversion

texts:
I have a Raster image like this

as you can see the road are not connected between 2 raster (up and down) its because the scanned raster use different projection (dgn utm 50s) than my data frame (utm 50s) when I'm doing Georeferencing. Somebody in here (Michael Mils) help me with some python code to create batch raster convertion without change raster name (I modified it a little to match the original raster projection "dgn utm 50 s").
Code :
Code are Fixed by Michael Mills, Problem Solved, Here's the Complete Code

import arcpy, sys, os

InFolder = r"d:\00_Start"
OutFolder = r"d:\01_End"
OutSR = arcpy.SpatialReference(23880) # "DGN_1995_UTM_Zone_50S"

arcpy.env.workspace = InFolder

for Ras in arcpy.ListRasters():
    Name,Ext = os.path.splitext(Ras)
    OutRas = Name + ".tif"
    arcpy.AddMessage("Sementara Proses " + Ras + " Coy")
    arcpy.ProjectRaster_management (InFolder + "\\" + Ras , OutFolder + "\\" + OutRas,OutSR)
    arcpy.SetRasterProperties_management(OutFolder + "\\" + OutRas,"GENERIC","", "", "3 0;2 0;1 0")
    arcpy.AddMessage("So Klar")

and the result is like this :

as you can see the road are connected now and everything's are perfect in place but there are still anoying black border in each raster (all the raster I converted result black border, there are 170 raster I need to process). How to convert using that python script without getting that black border?


texts: Having trouble getting specific items from a list. What am I doing wrong?

texts:
I am a total beginner when it comes to python.  I received some help already getting this code started but I hit a wall.  I have a directory with 200 folders, each with a geodatabase, feature dataset, and point line and polygon feature classes.  Each feature class has a field "Label" in the attribute table.  I need to be able to run a search through each feature classes' attribute table and return a specific string value.  For example, I need to find all feature classes where diabase is in the label field.  The output would give me only the feature classes that contain diabase in the label field.  I created a temp feature layer from each feature and did a select by attribute to see if get count returns anything.  Here's my code:

import arcpy
arcpy.env.workspace = "X:"
ws = arcpy.ListWorkspaces("*","Folder")
arcpy.env.overwriteOutput = True
tempLayer = r"in_memory\templayer"
field = "Label"
value = "diabase"
where = "%s = '%s'" % (field,value)
dhList = []
for w in ws:
     arcpy.env.workspace = w
     gdb = arcpy.ListWorkspaces("*", "Access")
     for fc in gdb:
         arcpy.env.workspace = fc
         fcl = arcpy.ListDatasets("*", "Feature")
         for fcc in fcl:
             arcpy.env.workspace = fcc
             fccl = arcpy.ListFeatureClasses("*", "All")
             for stil in fccl:
                 arcpy.MakeFeatureLayer_management(stil, tempLayer, where)
                 arcpy.SelectLayerByAttribute_management(tempLayer,"SUBSET_SELECTION")
                 count = str(arcpy.GetCount_management(tempLayer).getOutput)
                 print count
                 if count > 0:
                    dhList.append(stil)
for item in dhList:
     print item

The output is giving me a list of all of the feature classes in the directory instead of only the feature classes with diabase.  What am I doing wrong?

Thanks for any and all help in advance.


texts: What are good out of the box alternative geocoding solutions to Esri's services?

texts:
Esri's legacy geocoding services will retire at the end of the year, what are good GIS services that can replace the North American and US Streets geocoders?

I have totally excluded Esri's replacement, the ArcGIS Online Geocoding Service as it is a subscription service.


texts: Are Layer extensions supported in ArcGIS Pro?

texts:
Can I store some custom information about layer, project, etc?

For instance, featurelayer has no unique ID(I'm not 100% sure, but it seems like that) I want to add unique ID to each layer and store this information in aprx project file.

Can I do it using ArcGIS Pro SDK or any other way?


texts: Selecting by expression n features with lowest attribute value

texts:
I have a point layer in QGIS. The points have a distance attribute to a fixed point. I would like to select the three points with the shortest distance by an expression.

I want to do it by expression because I want to use this expression in some code I will use in PyQGIS.

The minimum function doesn't work over different features, so I am not sure how to do this.


texts: MapServer merging features from layers in Template-Driven Output

texts:
I want to get WMS GetFeatureInfo result in GeoJSON format from MapServer. Result should be valid GeoJSON FeatureCollection containing array of features from one or more layers.

I use Template-Driven Output (https://mapserver.org/output/template_output.html) with external template single file.

Template file (my.tmpl) is referenced in map file:

OUTPUTFORMAT
  NAME "geojson"
  DRIVER "TEMPLATE"
  FORMATOPTION "FILE=my.tmpl"
END

Template file (my.tmpl) looks like (two layers with Point geometry, having simple properties set):

// MapServer Template
{
 "type": "FeatureCollection",
 "features": [

[resultset layer=layer1]
[feature trimlast=","]
{
  "type": "Feature",
  "geometry": {
      "type": "Point",
      "coordinates": [ [shpxy precision=6] ]
  },
  "properties": {
    "id": [num]
  }
},
[/feature]
[/resultset]

<------ How to put ',' optional separator here?

[resultset layer=layer2]
[feature trimlast=","]
{
    "type": "Feature",
    "geometry": {
        "type": "Point",
        "coordinates": [ [shpxy precision=6] ]
    },
    "properties": {
      "id": [number],
      "name": "[name]"
  }
},
[/feature]
[/resultset]

]}

It works, but there is problem to put ',' separator between resultsets from several layers.
Using trimlast option in [feature] tag is essential to remove ',' separator after last feature in layer results.

Is there any support for conditional template output in MapServer? I.e. output string if some condition (for example, 'featureNumber == featureCount') is satisfied ?


texts: Setting value to defaults Layer properties using PyQGIS 3

texts:
I have a vector layer, how to set defaults value in fields from the Python console.

Layer Properties / Attributes Form / Defaults |

In option:
Default value |______________|


texts: Constrained Delaunay Triangulations in PostGIS?

texts:
It is possible to make Constrained Delaunay triangulations in PostGIS (PostgreSQL)? I have a Polygon and I need to make internal triangulation like this:

I tried to use ST_DelaunayTriangles but is not possible to add constraints like "triangulations only internal polygon"

Can someone help me?


texts: How can I clip this shapefile in Arcmap?

texts: I'm trying to clip one shapefile with another in ArcMap 10.1, but every time I try I get an error saying "overlay not polygon". I'm pretty new to GIS - is clip the right tool to be using here? I've also tried intersection (as suggested in other threads) but that doesn't work either.


texts: Creating Feature Layer from Parcel Fabric gives Error 000732?

texts:
Currently I am writing a python script that creates a feature layer from our parcel fabric and then appends that data into a different SDE feature class.  The problem is that when executing arcpy.MakeFeatureLayer the 000732 error is returned; dataset does not exist or is not supported.  I know the path to the fabric data is valid.  The fabric data can be accessed manually using that connection without trouble.  I also have a model for this type of process that uses the Make Feature Layer tool and uses the same SDE connection file which works without trouble.  I have found an abundance of information about this error online, but none of the suggested fixes seem to resolve my issue.

The code I am using is as follows:

    import arcpy
    parcels = r"\\networkMachineName\folder\projects\lgim\fabric.sde\fabric.DBO.ParcelEditing\fabric.DBO.ParcelFabric\fabric.DBO.ParcelFabric_Parcels"
    layer = "parcelLyr"
    whereClause = "SystemEndDate is null and Type = 7"
    arcpy.MakeFeatureLayer_management(parcels, layer, whereClause, "", "")

The parcel data from the fabric that is being input to the tool does not have any subtypes.

In the same network drive location I have other SDE connection files to different SDE geodatabases on the same server.  These connection files have been used as inputs to arcpy.MakeFeatureLayer using the same raw string convention with no trouble.  I only receive the 000732 error from the fabric SDE connection.  The only difference that I can see between the connections is that the other SDE geodatabases are SDE owner, but the fabric SDE geodatabase is DBO owner.  All of the connection files reference a user with full select, insert, update, delete privileges.

The ParcelEditing dataset is registered as versioned.

The error occurs in stand alone python script and in the ArcMap interactive python window.

When the model is first opened from ArcToolbox, the path to the fabric data shows the same error, but if the user browses to the data again the error is resolved and the model executes successfully.

Because the parcel fabric is a somewhat complex data type, I also attempted to execute arcpy.MakeParcelFabricLayer using the same connection file, but it failed with the same 000732 error.

We are using ArcGIS 10.2.2 for desktop, geodatabases, and fabric.  The RDBMS is SQL Server 2012.

I feel like this should be a straightforward task in python, but cannot figure out what is going wrong for some reason.  I am hoping that someone here can point me in the right direction.


texts:
I have a classification of vegetation communities across a landscape in a .tif format. I wanted to use the GRASS Zero-Crossing Edge Detection tool in QGIS to help define boundaries between my classes-at the moment the edges of each class are rather speckled and I would like a clear defined boundary. However when I run the tool, I get the following message: "The above files could not be opened, which probably indicates that they were not correctly produced by the executed algorithm". It then goes on to say that "GRASS seems to be correctly installed and configured".

I think the problem is possibly with the input file as in the help page it says the type should be "string". However if I convert my file to an .asc the tool will not let me select it as an input.

Has anyone used this tool before? If so, do you have any idea of what the problem is?


texts: Why is IQueryGeometry crashing when ITextSymbol:Text  contains more then one AMPERSAND (&) symbol?

texts:
I have a bit of an issue with my addin which converts annotations (text) to polygons.  To do this I use the ITextSymbol obtained from the Annotation Feature which is then implemented by the IQueryGeometry interface.

The actual conversion then occurs on this line:

pPolygon = pTextQuery.GetGeometry(pScreenDisplay.WindowDC, pDisplayTransform, pTextPointGeo)

All works good but I noticed it crashes on text strings that contain more then one instance of the ampersand symbol. (chr(38) I understand that this is a special character used for concatenation and such but I would like to retain it as a part of the string. I know how to test whether the string contains the ampersand but what should I do to retain the ampersands without crashing the IQueryGeometry interface?

In label expressions when I often encounter a similar issue I use FindLabel = Replace(FindLabel ,"&","&amp;") to remedy this. However when I try this in Visual Studio it still crashes the iQueryGeometry.  Strings that contain only one ampersand pass but strings with more then one ampersands crash.  By default REPLACE replaces all searched instances so I don't know what's going on here.


texts: I'm trying to use the modelbuilder to convert some feature classes to rasters using iterators. The names on each feature class contains year_xxxx with xxxx being an extension number that I don't want my output raster names to have. For example: 1845_6002pt.sh as input and I want 1845 only in the output.


texts: Download OpenStreetMap style sheet or mss file

texts:
I have download a shapefile from OpenStreetMap (OSM), then I open it using TileMill.  The next step is create stylesheet ".mss" file.

How can I download the exact OSM's stylesheet?


texts: I am new to mapping. Is there an online resource that lists both the geocentric and geodetic latitudes of various world landmarks?


texts: Importing/Plotting Lat/Long points from single column of *.csv using Python?

texts:
I have data in a .csv that looks something like this:

Entry A   Xa0,Ya0   Xa1,Ya1   Xa2,Ya2
 Entry B   Xb0,Yb0   Xb1,Yb1   Xb2,Yb2
 Entry C   Xc0,Yc0   Xc1,Yc1   Xc2,Yc2

The ultimate goal is to plot the points in each entry and connect them with lines, resulting in shape A, shape B, shape C, etc.

For now I am hung up on how to import these columns using arcpy. In the ArcGIS documentation there is code like this (and a tool that does the same thing) but it assumes lat/long are split up into separate columns:

in_Table = "table.csv"
x_coords = "POINT_X"
y_coords = "POINT_Y"
z_coords = "POINT_Z"
out_Layer = "points_layer"
saved_Layer = r"c:\output\points.lyr"

# Make the XY event layer...
arcpy.MakeXYEventLayer_management(in_Table, x_coords, y_coords, out_Layer, spRef, z_coords)

How do I modify this code to work with my data in the format it is in, that is, where a single column represents an X,Y pair?


texts: NDVI calculation of Color-Infrared Image in Arcpy

texts:
I have a color infrared image of an agricultural field. The color infrared image contains 3 different bands; NIR, Red, and Green (order maintained here).
To calculate the NDVI, I am currently using the "Image Analysis" tool of ArcMap (shown below).

However, I am writing a larger python script where calculating the NDVI from the color infrared image is a part. If I had two different single band rasters, I could have used the raster calculator. But, in this case, I only got one Color Infrared raster which contains 3 different bands. Is there any way to calculate the NDVI using the Color Infrared raster with ArcPy?


texts: Why graphical models from GRASS 7.6 don't work in GRASS 7.8?

texts:
I was working in some graphical model with grass and developed in GRASS GIS 7.6. Now I did an update with osgeo instaler from GRASS 7.6 to GRASS 7.8. When I try to run a graphical modeler nothing happen, It does not work.

Do I have to make my models again?

Does he graphical modeler from GRASS 7.6 work in GRASS 7.8?

I work in Windows 10.

This is the error I had.

Traceback (most recent call last):
File "C:\OSGEO4~1\apps\Python37\lib\sre_parse.py", line
1021, in parse_template

this = chr(ESCAPES[this][1])
KeyError
:
'\\E'
During handling of the above exception, another exception
occurred:
Traceback (most recent call last):
File "C:\OSGEO4~1\apps\grass\grass78\gui\wxpython\gmodeler
\frame.py", line 621, in OnRunModel

self.model.Run(self._gconsole, self.OnModelDone,
parent=self)
File "C:\OSGEO4~1\apps\grass\grass78\gui\wxpython\gmodeler
\model.py", line 689, in Run

self.RunAction(item, params, log)
File "C:\OSGEO4~1\apps\grass\grass78\gui\wxpython\gmodeler
\model.py", line 596, in RunAction

log.RunCmd(command=item.GetLog(string=False,
substitute=params),
File "C:\OSGEO4~1\apps\grass\grass78\gui\wxpython\gmodeler
\model.py", line 1211, in GetLog

cmd[idx] = pattern.sub(value, cmd[idx])
File "C:\OSGEO4~1\apps\Python37\lib\re.py", line 309, in
_subx

template = _compile_repl(template, pattern)
File "C:\OSGEO4~1\apps\Python37\lib\re.py", line 300, in
_compile_repl

return sre_parse.parse_template(repl, pattern)
File "C:\OSGEO4~1\apps\Python37\lib\sre_parse.py", line
1024, in parse_template

raise s.error('bad escape %s' % this, len(this))
re
.
error
:
bad escape \E at position 2


texts: Rasterio: error when resampling a raster

texts:
I'm trying to resample a Sentinel 2 band using an upscale factor. I'm using the code exactly as in Rasterio documentation but I keep getting an error.

Here is my code:

import os
import rasterio
from rasterio.enums import Resampling

Workspace = '.../SnowCoverProject'
S2_path = os.path.join (Workspace, 'raw_data/S2B_MSIL2A_20191208T184749_N0213_R070_T11UNQ_20191208T205518.SAFE')

#Browse through the S2 product folder and retrieve Green and SWIR bands
for root, dirs, files in os.walk(S2_path):
    for name in files:
        if name.endswith("B03_10m.jp2"):
            Gpath = os.path.join (root, name)
        elif name.endswith("B11_20m.jp2"):
            SWIRpath = os.path.join (root, name)
        else:
            pass

#Read the bands as arrays. Resample (upscale) SWIR to 10m.
upscale_factor = 2
with rasterio.open(Gpath) as green:
    GREEN = green.read()

with rasterio.open(SWIRpath) as swir:
    # resample data to target shape
    SWIR = swir.read(out_shape=(swir.count, int(swir.width * upscale_factor), int(swir.height * upscale_factor)), resampling=Resampling.bilinear)
    # scale image transform
    transform = swir.transform * swir.transform.scale((swir.width / SWIR.shape\[-2\]),(swir.height / SWIR.shape\[-1\]))

ndsi= (GREEN.astype(float) - SWIR.astype(float)) / (GREEN+SWIR)

profile = GREEN.meta
profile.update(driver='GTiff')
profile.update(dtype=rasterio.float32)

NDSI = os.path.join (Workspace, 'NDSI.tif')
with rasterio.open(NDSI, 'w', **profile) as dst:
    dst.write(NDSI.astype(rasterio.float32))][1]

I'm getting this error:

TypeError: read() got an unexpected keyword argument 'resampling'

Any ideas?


texts:
I have two set of polygons

select place_id,geom from table1;
select place_id,geom from table2;

I want to see which all polygons meet between these two tables. How do I go about this problem?


texts: Deselecting original feature from selection of both original and duplicate features?

texts:
I have a field which is supposed to hold a unique identifier, however due to issues with Attribute Assistant I have multiple features which have the same unique ID. I want to write a script which will select duplicates and remove the originals from the selection. The duplicated field is FID, the originals will have a lower OID (also, there can be more than one duplicate).

I have used SelectLayersByAttribute_management to select features where there are more than one feature with the same FID:

arcpy.SelectLayerByAttribute_management("MyLayer", "NEW_SELECTION", "FID In (SELECT FID FROM MyFeatureClass GROUP BY FID HAVING Count(*)>1))

This selects both originals and duplicates, but how do I now remove the features which have the original FID value?

The originals will have a lower OID value:

OID   FID
 1    x101    ORIGINAL
 2    x101    DUPLICATE
 3    x102    ORIGINAL
 4    x102    DUPLICATE
 5    x102    DUPLICATE

I'm using ArcMap 10.2.1


texts:
I'm making a plugin were i want to get the min and max value of a raster. When i tried the following code in the qgis python console it worked fine but in my plugin code i get the following error

 ver = provider.hasStatistics(1, QgsRasterBandStats.All)
 NameError: global name 'QgsRasterBandStats' is not defined

The code i m using is:

    layer = self.iface.activeLayer()
    renderer = layer.renderer()
    provider = layer.dataProvider()
    extent = layer.extent()

    ver = provider.hasStatistics(1, QgsRasterBandStats.All)

    stats = provider.bandStatistics(1, QgsRasterBandStats.All,extent, 0)

    minum= stats.minimumValue
    maxim = stats.maximumValue

    self.dlg.lineEdit.addItems(minum)
    self.dlg.lineEdit_2.addItems(maxim)

It look like QgsRasterBandStats isn't included in my python plugin. How can i include this?


texts:
I have downloaded numerous hdf files from an ftp site. Half of them were level 1 and the rest were level 2. Lvl1 format is Afile.A1_AC.sample.hdf and the lvl 2 format is: Afile.A1_AC.sample.hdf.

How will I do a batch process in projecting all the level2 files only?

How will I filter these files and create a new directory?

Here's what I have so far:

def inputData(self, event):

    def findSDS (path, filter):
        for root, dirs, files in os.walk(path):
            for file in fnmatch.filter(files, filter):
                yield file

     for sds in findSDS('/path/to/file', '*.hdf'):
        inSDS = 'HDF4_SDS:' + '/path/to/file' + '/' + sds + ':01'
        path = 'D:/new/dir'
        outTIFF = path + sds + '.tif'
        cmd = 'gdalwarp -geoloc -t_srs "+proj=longlat +datum=WGS84" -te 113.205 1.120 157.105 2.005 %s %s' % (inSDS, outTIFF)

        os.system(cmd)
        os.mkdir(path)

(tried calling .bat from the python script and run the script. And it just created an empty/new directory(out_path). Tried calling the .bat from the cmd line and here's the output:

C:\sample>project.bat
C:\sample>set in_path=D:\path\to\file
C:\sample>set out_path=D:\newfolder
C:\sample>md D:\newfolder
A subdirectory or file D:\newfolder already exists.
C:\sample>cd D:\path\to\file
C:\sample>FORFILES /m *A1_AC*.hdf /C "cmd /c gdalwarp -geoloc -t_srs "+proj=longlat +datum=WGS84" -te 113.205 1.120 157.105 2.005 @path D:\newfolder@fname.tif"
ERROR: Invalid argument/option - '+datum=WGS84 -te 113.205 1.120 157.105 2.005 @path D:\newfolder@fname.tif'.
Type "FORFILES /?" for usage.

C:\sample>project.bat
C:\sample>set in_path=D:\path\to\file
C:\sample>set out_path=D:\newfolder
C:\sample>md D:\newfolder
C:\sample>cd D:\path\to\file
C:\sample>FORFILES /m *A1_AC*.hdf /C "cmd /c gdalwarp -geoloc -t_srs "+proj=longlat +datum=WGS84" -te 113.205 1.120 157.105 2.005 @path D:\newfolder@fname.tif"
ERROR: Files of type "*A1_AC*.hdf" not found.

C:\sample>projected.bat
C:\sample>set in_path=D:\path\to\file
C:\sample>set out_path=D:\newfolder
C:\sample>md D:\newfolder
C:\sample>cd /d D:\path\to\file
C:\sample>FORFILES /m *A1_AC*.hdf /C "cmd /c gdalwarp -geoloc -t_srs "+proj=longlat +datum=WGS84" -te 113.205 1.120 157.105 2.005 @path D:\newfolder/@fname.tif"
ERROR 1: Translating source or target SRS failed:
'+proj=longlat

Tried running the gdalwarp in the cmd line and it was successful:

C:\sample>gdalwarp -geoloc -t_srs "+proj=longlat +datum=WGS84" -te 113.205 1.120 157.105 2.005 HDF4_SDS:hdf:"Afile.A1_AC.sample.hdf":01 sample.tif
Creating output file that is 1941P x 1670L.
Processing input file HDF4_SDS:hdf:"Afile.A1_AC.sample.hdf":01.
0...10...20...30...40...50...60...70...80...90...100 - done.


texts: Creating a label from two Attribute table column titles with a space between them

texts: would anyone happen to know the expression I need to type so that i'd have a space between two attribute columns titles in QGIS' Label expression dialogue?


texts: Get Substring of Linestring by nearest point

texts:
I need get custom segment of any Linestring. For better understanding I've added image.

Blue line is my Linestring which contains 3 points, and on the base of red point I need cut or fill line in same level as red point. How can I do that?

I use dijsktra function like that:

SELECT e.seq, e.path_seq, e.node, e.cost, e.agg_cost, w.name, w.the_geom, w.x1, w.y1, w.x2, w.y2, st_astext(w.the_geom)
FROM pgr_dijkstra('
    SELECT gid AS id,
         source,
         target,
         cost_s AS cost,
         reverse_cost_s AS reverse_cost
        FROM pgr.ways',
    35959,
    19303,
    directed := true)e
    left join pgr.ways w on w.gid = e.edge
    order by e.seq

and on image are my first and last edges of result which are not corresponding with ending vertexes (red points), so I need edit this edges.


texts: Determining distances between multiple points without crossing roads?

texts:
I'm working with civic addresses (points) in the county that I work for (there are nearly 10,0000 points). I am trying to find out how many of the points are within certain distances or less from each other. (Ex: 400 feet or less from each other).

The other challenge is that I do not want the calculation to include going across roads. They need to be on the same side of the road to meet the search criteria.

I've never done this sort of query so would really appreciate any help. A small snapshot of what I'm referring to is below.

I have access to full ArcGIS suite if necessary so have full range of tools at my disposal. I've done some coding before and can work with it, but I'm not good at creating it myself necessarily.

I'm really just starting to tackle this problem and haven't put too much into it yet. It seems that the NEAR and POINT DISTANCE might be the way to go?

I need to get access to the license to do these so am trying to do my research in advance to best utilize my time.


texts:
Just before submitting bug: in Diagrams, there is possibility to use Position AroundPoint or OverPoint. Whatever I choose Line Options selection stays gray: I can't choose any line.

So, I suggest that this Line Option -option will draw line between actual point and re-located diagram. However, it doesn't work in my QuantumGIS 1.8.0 (win or Ubuntu). I haven't test it on 1.9

Rgs,

Pekka


texts: Polygons (world countries) over 180 degrees gets broken in rendering

texts:
I have a table with world country polygons, and some polygons (as in real life) cross 180 degree longitude line, so they have both 179... and -179... coordinates.

As long as it is in PostGIS or exported to Shapefile, then it is displayed ok in QGIS.
Problem 1: When I converted it to GeoJSON, then it got broken, and geojson.io and QGIS both shows it broken. See image below
Problem 2: When I imported GeoJSON or Shapefile to CartoDB, then it was broken also. I suspect CartoDB does internal conversion to GeoJSON which distorts it.

My ultimate goal is to get it working in CartoDB, so it is rendered correctly. How can I do it?


texts: Merge multiple feature datasets while maintaining shapefiles for respective point, line, and polygon files nested ArcGIS

texts:
I've been doing some digging and a few days of playing around with my data in ArcInfo but I haven't been able to find a solution to merging datasets, which consist of many different point, line, and polygon shapefiles.

The files are electronic navigation charts in a vector format (some examples of these are freely available: http://www.charts.noaa.gov/ENCs/ENCs.shtml). Each of these charts covers a different spatial extent or resolution and contains information relevant to mariners, fishermen etc. (e.g., depth, substrate type, buoys, lights, shipwrecks, obstructions, and so on). There is a standard for the display of these data called S-57 in which each feature is universally coded (e.g., obstructions are OBSTRN in each ENC file: http://www.s-57.com/). The coding of the attributes of each feature is also standardized. Therefore, each file you open in Arc will have identically named point, line and polygon layers for the features on the chart (although some charts may have features not found on other charts [e.g., some areas might not contain shipwrecks]). I am wanting to merge many of these ENC files, which cover portions of my region of interest, into a large regional layer. My hope is to use the respective shapefiles (and attributes) for this whole region rather than for the smaller areas the charts cover individually. The merge has to maintain all the respective point, line, and polygon shapefiles (e.g., depth, substrate type, buoys, lights, shipwrecks, obstructions) (i.e., it is basically a merge of all identically named shapefiles nested within each chart).

I know the merge can be done between two respective shapefiles from 2 different ENC chart datasets (after a name change of one of them), however, I have potentially 20-30 charts to merge and there are 20+ shapefiles in each chart file.

Any direction would be greatly appreciated. Thanks in advance.

(I attached a screenshot of the ENC file in ArcInfo. Each chart file has these identical layers covering a different spatial extent I wish to merge).


texts: Info (data from fields) from multiple overlapping features in popup element in ArcGIS JavaScript API

texts:
I'm just trying to do a simple task. I have a point file with numerous overlapping points. When I load the service onto ArcGIS.com (their online mapping site) and enable popups and click on a point that is lying on top of other points, it gives me a popup, and gives me a little back and forward arrow to run through each of the points selected through my click.

You can sort of see it here in the image: https://developers.arcgis.com/en/javascript/jshelp/intro_popuptemplate.html

The image on the right hand side of that page shows (1 of 2) in the popup window title element.

The way my selection is done is by simply adding to the layers as a featureclass to my dojo map element.

So basically, what my question is, how do I enable and do this? Where when I click on a point with underlying points, the popup allows me to view those as well (instead of the point that's just on top).

Thanks


texts: Get field names of shapefiles using GDAL

texts:
I use GDAL in Python for importing shapefile. I want to know the field names for the file, my current way is:

fields = []
for i in range(1, layer.GetFeature(0).GetFieldCount()):
    field = layer.GetFeature(0).GetDefnRef().GetFieldDefn(i).GetName()
    fields.append(field)

But this way, I am getting the feature for the first layer. Does it mean it is possible that different layers can have different features?

If not, is it possible to get the field names at once, instead of getting into this deep? If yes, is there any easier way of getting the field names?


texts: How to change symbology on layer only available on Esri REST service?

texts:
I'm trying to work with a layer I can only find from an Esri REST service and I need to change the color of the layer, as the default blue it comes in doesn't show up well against an already-established map we're adding it to. When I open layer properties, I don't see the usual symbology tab or any other symbology options on any tabs.

The symbology/legend seems to be built into the layer API from what I see at this page on the file details: https://gis.wim.usgs.gov/arcgis/rest/services/AIR_NDGA/AIR_NDGA/MapServer

Is there an option to change symbology on layers like this? Here's the layer of interest that seems to have fixed symbology: https://catalog.data.gov/dataset/bia-indian-lands-dataset-indian-lands-of-the-united-states I have not found any other sources of this dataset other than through Esri REST API links.

One option I've explored to no avail is converting / saving the API-provided layer into a regular shapefile saved locally, figuring I could change the symbology on that.


texts: GEE - Add date to band names in ImageCollection

texts:
I'm new to the GEE platform and I'm struggling with something very basic. I want to extract pixel values by points from an ImageCollection and I want my ouput to be a table in a wide format. To do this and apply my reducer, I exploded the ImageCollection and, instead of having a number X of images with Y bands for each image, treated all rasters as bands of only one image.

Here is the code:

var band_Names = ee.List(['Blue', 'Green', 'Red', 'NIR', 'SWIR 1', 'SWIR 2', 'pixel_qa']);

var LC08_collection = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
  .filterDate('2013-01-01','2019-12-31')
  .filterBounds(points)
  .map(function(image){
  return image.rename(['Aerosol', 'Blue', 'Green', 'Red', 'NIR', 'SWIR 1', 'SWIR 2', 'TIRS 1', 'TIRS 2', 'sr_aerosol', 'pixel_qa', 'radsat_qa']);
  })
  .select(band_Names);

var mergeBands = function(image, previous) {
  return ee.Image(previous).addBands(image);
};

var collection_asImage = LC08_collection.iterate(mergeBands, ee.Image([]));
var collection_asImage = ee.Image(collection_asImage).float();

var extraction = collection_asImage.reduceRegions({
  collection: points,
  reducer: ee.Reducer.first(),
  scale: 30,
});

Export.table.toDrive({
  collection: ee.FeatureCollection(extraction),
  description: 'Pixel values',
  fileFormat: 'CSV',
});

All is well and good, but, in this way, the field names of my dataset will look like this:

'Id' | 'Class' | 'Blue' | 'Green' | 'Red' | 'NIR' | 'SWIR 1' | 'SWIR 2' | 'pixel_qa' | 'Blue_1' | 'Green_2'...
1
2
3
4
5
6

Which is great, except for the fact that after the first set of bands there is the number that underlines the iteration of an element ('Blue_1'), which I don't want. Instead, I want to have the date of the image from which the value of that band comes from ('Blue_2013-05-20' for example). I got the date as a list in this way:

var dates = LC08_collection
  .map(function(image) {
    return ee.Feature(null, {'date': image.date().format('YYYY-MM-dd')});
  })
  .aggregate_array('date');
print(dates, 'Date list');

But I don't know how to rename the bands. I thought of renaming them before exploding the collection and merging all the bands, writing a function that will be mapped over the ImageCollection which returns the band names already there and adds the element of dates array, but I don't know how to achieve this result. Does anybody know how to help me?

Thanks in advance.


texts: Calculating forest canopy density and advanced vegetation index in ArcGIS Raster Calculator?

texts:
While performing raster calculator operation in ArcGIS for Advanced Vegetation Index ( in calculating forest canopy density), the results show only 2 values i.e high (1) and low ( 0), instead of a range of values.The low values considered are of that of the background no-data values.

What could be the possible problem in the operation ?

I am using the formula :

AVI = (("ref5"+1) * (256-"ref4") * ("ref5"-"ref4")) ^ (1/3)

where ref5, ref4 are bands 5, 4 corresponding to Landsat 8 image, converted into TOA reflectance values.


texts: ArcPy PrintMap is cutting off part of layout?

texts:
I am running a script that pulls certain .mxd files out of a folder and prints them for me.  For some reason, when I use the script, part of the layout is cutoff; on the other hand, when I open each .mxd individually and print, everything comes out fine.  Is there a way to control the print setup within my script to correct this?

    import arcpy, os
    from arcpy import env
    env.workspace = r"folder directory"

    ws = env.workspace

    locator = arcpy.ListFiles("*LOCATOR.mxd")
    for i in locator:
        filepath = os.path.join(ws, i)
        mxd = arcpy.mapping.MapDocument(filepath)
        arcpy.mapping.PrintMap(mxd)

    vicinity = arcpy.ListFiles("*VICINITY.mxd")
        for k in vicinity:
        filepath = os.path.join(ws, k)
        mxd = arcpy.mapping.MapDocument(filepath)
        arcpy.mapping.PrintMap(mxd)


texts: Convert polyline to polygon Using ArcObjects

texts:
I have a polyline and I want to convert it to a single polygon feature, I used this code but I get error while I can convert the same shapefile using feature to polygon on ArcGIS.
Where is my mistake?

ESRI.ArcGIS.DataManagementTools.FeatureToPolygon ftp = new ESRI.ArcGIS.DataManagementTools.FeatureToPolygon(path + "lines.shp", path + "polygons.shp");
        ftp.attributes = "NO_ATTRIBUTES";
        GP.Execute(ftp, null);

The path is set to

string path= System.IO.Path.GetDirectoryName(Application.ExecutablePath) + @"\temp\Conversations\";

and the error that raises is:

Error HRESULT E_FAIL has been returned from a call to a COM component.

I tried to change output path to

"E:/MASTER/Project/Semirom/file gis semirom/New Folder/gisiran.gdb/tester"

But I get the same error... Do you think I must set Environment workspace or something like that? I also use:

ESRI.ArcGIS.ConversionTools.GPXtoFeatures and ESRI.ArcGIS.DataManagementTools.PointsToLine(path + "points.shp", path + "lines.shp"); with the same settings and they works fine..:(

Is there any better approach to create a polygon?


texts:
I am taking the Learning Spatial Analyst course, which is very good but written for Version 9. I and my students are using Version 10, so some of the syntax for Raster Calculator have changed. Can anyone assist with translating these commands to Raster Calculator or point me to a document that updates them for Learning Spatial Analyst 9? The concepts are being emphasized in the course but its hard to get them without seeing the exercise results.

CON([DamBasin] AND ([DamElev] <= 790), 790)

CON(ISNULL([Reservoir]),[DamElev],790)

SetNull([Ownership] <> 200 OR [Soil] == 901, 1)


texts: Can you use ArcObjects in an MVC service at v10.1

texts:
I'm slighlty confused with the changes at 10.1 and understand this version no longer supports DCOM or remote arcobjects.

Can you use ArcObjects in a webservice (in my case a .NET MVC 4 service) that is deployed on the ArcGIS Server host (so not remote) and uses ArcObjects to perform editing on a gdb featureclass?

Or is the only way to go through the SOE framework?

If you can use ArcObjects in the .NET MVC service are there any ramifications of doing so over using an SOE?


texts: QGIS temporary layer to PostGIS table database

texts:
I have an empty polygon table in PostgreSQL/PostGIS database with fields (id,field_1).

after some processing in QGIS using pyqgis and processing toolbar I create some temporary layer like this case :

outputs_GDALOGRPOLYGONIZE_1=processing.runalg('gdalogr:polygonize', raster,'DN',None)
outputs_GDALOGRCLIPVECTORSBYPOLYGON_1=processing.runalg('gdalogr:clipvectorsbypolygon', outputs_GDALOGRPOLYGONIZE_1['OUTPUT'],cut,None,clip)

how to add the temporary layer outputs_GDALOGRPOLYGONIZE_1['OUTPUT'] into the empty polygon table from PostGIS table using pyqgis?


texts: How to clip a spatial polygon to a large matrix (that is plotted as a large trellis, originated as a netcdf) in R?

texts:
I would first like to issue a disclaimer that I am a novice at GIS, and I have only just begun learning R. So if I use the wrong terminology for something or my question is confusing, please help clear that up for me.

So for my work, we are planning on mapping tons and tons of weather data all across the globe. For practice, the guy providing our data provided us with some "example" data to play with. My boss would like for me to create my maps in R. These maps are intended to be return period maps (https://en.wikipedia.org/wiki/Return_period). The data is supposed to be in netcdf format.

I've managed to find several tutorials that have helped me to actually map the netcdf data. Here's the code:

#set working directory
setwd("D:\\stuff")

#add libraries
library(maptools)
library(ncdf4)
library(RColorBrewer)
library(lattice)
library(rasterVis)

#add madagascar shapefile
proj <- CRS('+proj=longlat +ellps=WGS84')
madagascar <- readShapeSpatial("MDG_adm0.shp", proj4string = proj)

#read netcdf data for 10 year return period
ncdf.data <- nc_open("swio_rpmaps_200_83.nc")
returns <- ncvar_get(ncdf.data, "y010")

#define latitude and longitude
lon <- ncdf.data$dim$longitude$vals
lat <- ncdf.data$dim$latitude$vals

grid <- expand.grid(lon = lon, lat = lat)
cutpts <- c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90)
test <- levelplot(returns ~ lon * lat,
          data = grid,
          at = cutpts,
          cuts = 10,
          pretty = T,
          col.regions = rev(brewer.pal(10, "RdBu")),
          ylab = "Latitude",
          xlab = "Longitude",
          main = "10 Year Return Period")

#return periods with madagascar plotted on top
test + layer(sp.polygons(madagascar))

Here's what the map looks like:

Problem: I want to be able to clip out JUST the Madagascar portion of the windspeed data presented in this map. I cannot for the life of me figure out how to do this. I attempted to use gIntersection from the rgeos package, but the netcdf data (which shows up as a large matrix in my workspace; when I define it as "test" and use levelplot, it becomes a large trellis in the workspace) is apparently not spatial, so it won't let me do anything.

Is there any way to perform this clip with the code I've written? If not, how would I be best suited to overlay a spatial polygon with this netcdf data? I'm in WAY over my head here!


texts: How to generate a map projection in which equidistance is maintained along routes?

texts:
In this SE.History question, an answer had shared a map from this PDF,     ,showing two shipping routes between London and the Arabian Gulf.

I assume that these shipping routes aren't equidistant.  This is, I assume that, if we trace a route at equal lengths, it wouldn't correspond to equal travel distances across the Earth's surface.

I'm curious what such a projected map would look like, but I don't seem to see an appropriate projection in Wikipedia's list of map projections.

Specifically, the projection should:

Maintain equidistance along arbitrarily drawn routes.
Seek to minimally distort the other parts of the map, especially those closest to the routes.

Since such a projection would depend on the routes selected, it'd need to be more of a dynamic projection where the map of Earth looks different depending on the routes being displayed.

Question:  Is there an algorithm that creates projected maps in which a selection of routes maintains equidistance?


texts: Displaying a map in EPSG:26915 projection in OpenLayers 3

texts:
My GIS experience is limited to having worked on several Geospatial Humanities projects. That is to say that I still struggle significantly with the technology (and terminology) in getting to my goal. I am working on moving a project from OpenLayers2 to OpenLayers3. I have a historical map served as a WMS from an Iowa State ArcServer that is in EPSG:26915 projection. The following code is what I use to successfully display the map in OpenLayers2:

var map;
function init() {
    map = new OpenLayers.Map('map',{
        projection:'EPSG:26915',
        maxExtent: new OpenLayers.Bounds(202159,4470000,737457.523439407,4822691.7712555),
        fractionalZoom: true
    });
    var wms = new OpenLayers.Layer.WMS(
        'OpenLayers WMS',
        'http://ags.gis.iastate.edu/arcgisserver/services/Ortho/andreas/ImageServer/WMSServer',
        {layers: 'basic'},
        {isBaseLayer: true },
        {}
    );
    map.addLayer(wms);
    if (!map.getCenter()){
        map.zoomToMaxExtent();
    }
}

I have a number of unsuccessful attempts to display this using OpenLayers3, trying to piece together what I need from the ol3 examples and documentation. Here is my most recent unsuccessful attempt:

var map;
function init() {
    var andreas_url = 'http://ags.gis.iastate.edu/arcgisserver/services/Ortho/andreas/ImageServer/WMSServer';
    var andreas_extent = [202073.808723, 4470598.397186, 736849.210393, 4822673.786539]
    var andreas_projection = new ol.proj.Projection({
        code: 'EPSG:26915',
          // The extent is used to determine zoom level 0. Recommended values for a
          // projection's validity extent can be found at http://epsg.io/.
        extent: andreas_extent,
        units: 'm'
    });
    ol.proj.addProjection(andreas_projection);
    var andreas_view = new ol.View ({
        projection: andreas_projection,
        center: [469163,4653689],
        zoom:1
    });
    map = new ol.Map({
        target:'map_container',
        renderer:'canvas',
        view: andreas_view
    })
    var andreas_layer = new ol.layer.Tile({
        source: new ol.source.ImageWMS({
            url: andreas_url,
            params:{
                'LAYERS':'basic',
                'TILED':true
            }
        })
    });
    map.addLayer(andreas_layer);
}

The WMS Capabilities URL for this map is here:

http://ags.gis.iastate.edu/arcgisserver/services/Ortho/andreas/ImageServer/WMSServer?Request=GetCapabilities&service=WMS&

Can anybody see what I might be doing wrong, or does anybody have a good resource for learning about projections in OL3?


texts:
How can I merge shapefiles in a special order in QGIS?
What I have and what I need to be merged in that order:

What happens when using merge by mmqgis:

As u can see not usable at all. Attributes are gone.

What happens when i use "merge shapefiles" tool by qgis:

As u can see some time zones are covering smaller time zones. So i tried to order them by using mmqgis again:

Attributes are ordered now but geometry is still overlapping.

Is it possible to merge shapefiles in a special order using qgis respectively to sort their geometry afterwards?


texts:
I have a large point dataset that I am publishing as a service. When zoomed out I would like to show the data as a heat map, only showing the actual points when zoomed right in. I know there are a couple of API's that can do heat map rendering client side (java, leaflet, Silverlight, etc) and the esri maps for office has a nice example as well. But is there any way to utilise a dynamic heat map renderer within ArcMap, maybe using python?

Currently I am using the point density tool to generate multiple raster layers with varying radius's and using scale dependencies to flick through them as you zoom in, you get a similar effect but this seems a little clunky and I am hoping there is a better way?


texts:
I have read Overpass API by example, but still I am unable to adapt convert statement examples to my need. What I would like to do is to get all marked routes in the area using Overpass Turbo, drop some tags and synthesize some new tags.

My query looks like this:

relation[route][route!="power"][route!="road"]({{bbox}});

out body;
>;
out skel qt;

I would like to:

add tag color=green when osmc:symbol contains green or when tag kct_green exists
add tag color=blue when osmc:symbol contains blue or when tag kct_blue exists
add tag kind as a value of any tags of kct_green, kct_blue, kct_red, kct_green
remove some tags I no longer need to reduce output file size (like kct_green)

My initial naive attempt is:

relation[route][route!="power"][route!="road"]({{bbox}});

foreach (
  convert rel ::=::,
    ::id = id(),
    color=t["osmc:symbol"]
  ;

  out;
);

out body;
>;
out skel qt;

This kind of works, but lists relations without including their nodes. Moreover, I would like to perform a bit more processing on the color value:

I want the color to be created only osmc:symbol is present
I want to extract some information from osmc:symbol to color
I have no idea how to create the combined value from those kct_green... tags.

I guess I need something more complex, however I have found almost no documentation or examples about convert. Can someone point me in the right direction?


texts: What is the maximum number of features a SDE featureclass can hold?

texts:
I have a SDE 10.0 Geodatabase based on MS SQL 2008.

I'm running a continuous process where about 6400 point features are written to a feature-class every hour.  I need it to be a featureclass, since I'm doing some krigging on the feature-class, after applying an attribute query based on a time field in the data.

I was wondering what is the limitation on the number of features a SDE Feature-class can hold?


texts: Calculation based on features in shapefile

texts:
I hope I get the terminilogy right! Correct me if I mess something up.
I need to perform some calculations for a scoring system for different habitat types.
The Problem is that I have two different tables in two different shapefiles. They need to be divided through each others somehow (yet again, all needs to be done automatically in the modelbuilder).

For the calculation I need to Sum up the Size of all the Types with a common/matching Name in Table A (I tried to clarify this through the field Sum).

Shapefile A with Table A
Name    Type    Size   Sum
A       BO      10     ?
A       BE      20     ?
C       BA      30     ?
D       BO      10     ?
D       FW      50     ?

The Sum for Name A would be 30 (10+20), for D it would be `60 (50+10).

In my second table, I do need not to do any calculations. It looks likes this:

Shapefile B with Table B
Name    Type    Size
A       BO      20
A       BE      50
C       BA      60
D       BO      90
D       FW      10

I know the Size of each Type in in each Name.

The Size of Table B needs to be divided through the Sum of all the Types for each Name in Table A.
I hope you can see the problem I am running in. How do I tell the Modelbuilder to summarize all the Size of all Types that have a matching Name? And how would I populate Fields with this outcome for further calculations?


texts: Comparing raster - pixel colors using PyQGIS?

texts:
I would like to compare satellite pictures (raster) infra - red.

I have 2 rasters from the same area. I would like the compare pixels with the same xy coordinate using PyQGIS.

10 different values between the same= O and totaly different = 10.
10 different red colors.

So I can compare leaf green the state of the vegetation.

How do I start and which why to go?


texts: PyShp error: Object does not implement __geo_interface__

texts:
I'd like to convert geometry to WKT with pyshp and pygeoif using this suggestion, but pyshp throws a following error:

Traceback (most recent call last):
  File "E:/!Projects/imScrape/metadata_extractor.py", line 82, in <module>
    wkt = extract_wkt_from_shp(shp_fileobj, shx_fileobj, dbf_fileobj, prj_fileobj)
  File "E:/!Projects/imScrape/metadata_extractor.py", line 36, in extract_wkt_from_shp
    print(pygeoif.geometry.as_shape(geometries[0]))
  File "C:\Users\lobanov\AppData\Local\Programs\Python\Python36\lib\site-packages\pygeoif\geometry.py", line 1163, in as_shape
    raise TypeError('Object does not implement __geo_interface__')
TypeError: Object does not implement __geo_interface__.

What might be the case?

Here's the code:

shape_file = shapefile.Reader(shp=shp_fileobj, shx=shx_fileobj, dbf=dbf_fileobj)
geometries = shape_file.shapes()
geom = [pygeoif.geometry.as_shape(geometries[0])]
wkt = pygeoif.MultiPolygon(geom)


texts: Importing large tiled datasets (e.g. Ordnance Survey) into MapInfo?

texts:
I am currently using MapInfo and have several large Ordnance Survey datasets, which require updating. They are tiled .tiffs and need to first, be registered as .tab files. Second the .tab files need to be imported to mapinfo and thirdly turned into a seamless table.

The process I used is quite longwinded:

I copy the paths of the .tiff files and paste into excel and build the columns around it to create this line for each file e.g.

**Register Table "X:\Ordnance Survey\25k_raster\nx\nx15.tif" type "raster" into    "X:\Ordnance Survey\25k_raster\nx\nx15.tab**

I add this list to the workspace file using a text editor and then I build a line in excel which then opens each of the tab files e.g.

**Open Table "X:\Ordnance Survey\25k_raster\nx\nx15.tab" As nx15 Interactive**

I add this second list to the bottom of the workspace file, save then open the workspace. It takes a while but eventually it has registered and opened all the files. I then proceed to create a new seamless table from the open tables in mapinfo.

I guess my question is. How do you do this kind of import with mapinfo? do you know of or use any tools which would do this process for you? Maybe you have a quicker workaround?


texts: Using ArcGIS Marketplace for focussed apps with limited functionality?

texts:
According to the ArcGIS Marketplace FAQ:

  All apps in the marketplace must require a named user login.

and

  Can people who aren't members of an ArcGIS Online organization purchase an app from the marketplace?

  No. The marketplace is for ArcGIS users, so access to apps or data is
  available exclusively to ArcGIS Online-named users.

According to the ArcGIS Marketplace Help:

  There are two main requirements for building ArcGIS Marketplace apps:

  Leverage ArcGIS Online subscriptions ...
  User login with OAuth 2.0 ...

I have an idea for an app that leverages ArcGIS Online subscriptions and requires user login with OAuth 2.0.

I am hoping to bring it to market at a price of about $20 per annum because it is a focussed app targeted at a potentially large but non-traditional market.

However, when I read about named users on the Identity page it says that:

  With your identity [(named user credential)], you can access a growing
  list of standard and premium capabilities:

  Ready-to-Use Apps, App Builders, and Connectors
  Collector for ArcGIS
  Operations Dashboard for ArcGIS
  ArcGIS Open Data (for publishing)
  Web AppBuilder for ArcGIS (for building apps)
  ArcGIS Maps for Office
  ArcGIS Pro

My understanding is that named users have been priced according to the most sophisticated functionality in this list (ArcGIS Pro) when my app would be using only the most basic functionality in this list (Ready-to-Use Apps, App Builders).

Is there a pathway to participate in the ArcGIS Marketplace using named users in the sense of an identity that makes them known to Esri (OAuth 2.0), which leverages ArcGIS Online subscriptions by enabling its credits to be drawn down, but only allows those named users access to very limited functionality at a commensurate price?

I ask because I have not been able to find one.


texts: What do version numbers in CRS IRIs mean?

texts:
By following the links on http://www.opengis.net/def/crs/ one can download identifiers/IRIs for various CRSes minted by the OGC.  For many CRSes, there seem to be multiple identifiers, differing in what looks like a version number.  E.g., for EPSG:2000:

http://www.opengis.net/def/crs/EPSG/0/2000
http://www.opengis.net/def/crs/EPSG/8.5/2000
http://www.opengis.net/def/crs/EPSG/8.9.2/2000
http://www.opengis.net/def/crs/EPSG/9.4.2/2000
http://www.opengis.net/def/crs/EPSG/9.5/2000
http://www.opengis.net/def/crs/EPSG/9.5.3/2000

Does anybody know what these version-like numbers refer to?  The OGC GeoSPARQL standard uses IRIs with 0 (i.e., the first IRI in the above list), but does not explain what the 0 stands for.


texts: Highlighting Feature From Selection Set using ArcObjects?

texts:
Is there a way in ArcObjects 10.4 that I can highlight a specific feature from a selection set.

I'd like to give visual feedback to my users when they are highlighting a feature in my form. I want the currently selected feature in my selection set to be highlighted (Yellow).


texts: How to convert data from GML to GeoJSON using ogr2ogr?

texts:
How should I convert GML to geojson data using ogr2ogr.
I used :

ogr2ogr -f "GeoJSON" bd.geojson bd.gml

But in this case I have only part of my data , not all.
I tried add -append , but it didn't work properly.


texts: Plotting on MapInfo given edge geometry in row?


texts: Landsat 5 RGB image looks very dark in QGIS

texts:
When I combine B-3,2,1 by "Raster>>Miscellaneous>>Merge" in QGIS 2.14.3 I get very dark looking images but when I use the same bands in Sentinel-2 toolbox it looks much better.

Landsat-8 data works fine with "Raster>>Miscellaneous>>Merge" too.

Any idea/suggestion why is this happening?


texts: Using Aggregate Points to create heat map?

texts:
I'm trying to create a heat map using Point Density but the process keeps failing by cycling through the following errors:

Warning 000597: Skipping Feature because of NULL or EMPTY Geometry

ERROR 010429: ERROR in GRID IO:S_XXfer: Read Error.

ERROR 101244: Could not set the analysis window.

ERROR 010067: Error in executing grid expression.

Failed to execute Point Density.

I tried to aggregate the points then apply the Morans I tool but this has been running for over 2 hours with no sign of ending.

Does anyone have suggestions to get smooth heat map outputs that are not blocky/ pixelated?


texts:
The scale of the canvas plot shows in the status bar as everyone knows. However If I want to change the scale it refuses to work properly.  This occurs whether I select a scale from the drop-down menu of standard scales or type in my own desired scale.  QGIS seems to have its own set of irregular values such as 1:577,790 and 1:144,447 and refuses to accept that I type in or select.
Help, please.


texts: convert lat long coordinates in a table to MGRS

texts: I have an Excel table with lat and long coordinates in separate fields. Now I'd like to convert these fields into a third field with its MGRS coordinate. Any suggestions of how?


texts: How to assign numerical value to polygons?

texts: I have a shapefile of polygons, and a list of mean values of the area in each polygon. I need to assign said mean values to each polygon. How would I go about this using GDAL/OGR in Python?


texts: Scriptable alternative to ArcScan (particularly vectorising lines within tolerances)?

texts:
Sadly the answer to Can you access ArcScan from a python script? suggests that it is not possible to script ArcScan from within a python script.

Therefore, does anyone know of any alternatives to ArcScan (either within ArcGIS, or within another GIS, or anything really...) that are scriptable? The functionality that I am in need of is the ability to convert from raster to vector by 'intelligently' vectorising lines. For example, ArcScan allows you to set a distance tolerance so that any lines that are vectorised and are less than that distance apart will be joined together. It also does various other sensible things with raster lines (for example, if a line is 3 pixels wide it vectorises it to a single line placed in the middle of the raster line).

Any ideas?


texts: Why is my .prj file not being read by QGIS?

texts:
First off, I am new to GIS.
I am using qgis (1.8.0), and have been given a shapefile containing .shp, .prj, .shx, .dbf, however qgis will not load the .prj file stating an invalid data source.
Altering the project properties>CRS to WGS84 Epsg:4326 gives correct latitude for the project but incorrect longitude.

The .prj file reads:

  PROJCS["British_National_Grid",GEOGCS["GCS_OSGB_1936",DATUM["D_OSGB_1936",SPHEROID["Airy_1830",6377563.396,299.3249646]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",400000.0],PARAMETER["False_Northing",-100000.0],PARAMETER["Central_Meridian",-2.0],PARAMETER["Scale_Factor",0.999601272],PARAMETER["Latitude_Of_Origin",49.0],UNIT["Meter",1.0]]


texts:
I have the following script:

import arcpy
from arcpy import env
import datetime
import os

env.workspace = r"Q:\\SDE\Direct Connection to Contours.sde"
dataFC = arcpy.ListFeatureClasses("","","Contours")
outpath = r'C:\Users\TEST\Desktop\SDE Conversion Project\TEST.gdb'
query = "objectid > 0 AND objectid <= 200000 "

for fc in dataFC:
        if fc == 'Contours':
            output = outpath + os.sep + fc
            arcpy.Select_analysis(fc, output+'200k', query)

Currently it selects records that have objectid between 0 and 200,000 and creates a feature class called Contours 200k in the TEST gdb.

What I want to be able to do is make a multiprocessing script that utilises one core to calculated between object id 0 and 200k, another that calculates betweek 200,001 and 400k and so on until the whole feature class is copied over as separate feature sets to the TEST gdb.

How would I go about doing that?
I'm trying to do this since the bottleneck of copying from the sde is the network connection and testing with separate files where I had the ranges set manually improved the overall speed a lot.


texts: Convert all rasters in a rasterstack to polygons

texts:
I need to convert all rasters in a RasterStack to polygons and give the polygons names according to its respective raster.

library(raster)

## r is a single raster
rpolygon <- rasterToPolygons(r, dissolve = TRUE)

## rstack is a raster stack
lapply(rstack, rasterToPolygons)

But this lapply function returned the following error.

Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘nlayers’ for signature ‘"numeric"’


texts: How do you create a button in ArcGIS and have it run code?

texts:
I'm brand new to ArcGIS and have only a year or two of programming experience. I know Java and a little C#. I have a Java app that processes geographical data. I need to have a button in ArcMap that, when you click it, sends map data to this Java app.

This task was already implemented in VBA; it's my job to translate it to something else since VBA is no longer supported in ArcGIS. Any thoughts on how to accomplish this?

I'm not sure where to begin. I'm guessing either using the embedded Python console or the ArcObjects SDK for Java or .NET. Just want to see what someone with experience would recommend.


texts: Schema Lock Test in Model Builder for ArcGIS 10.0 Desktop

texts: I have made a model (Model Builder) that exports a Polygon Layer from a Parcel Fabric and Exports it as a set of Features Classes into a set of Feature Datasets. However, there are often other people working within the datasets. I would like to be able to check the dataset for a schema lock before trying to save the data, and then if possible keep checking until the schema lock is removed and the model can complete the process. Is this possible within model builder, and how?


texts: Conducting Land Cover Signature Classification on Google Earth Engine?

texts:
I was wondering if there is a way to conduct Land Cover Signature Classification, as specified on page 122 of the following user manual: https://buildmedia.readthedocs.org/media/pdf/semiautomaticclassificationmanual/latest/semiautomaticclassificationmanual.pdf, on Google Earth Engine.

Essentially, my problem consists of the following: using a single image, I am trying to find pixels whose spectral curves are contained within a certain threshold from my training ROI's spectral curve. It is almost equivalent to performing a one-class supervised classification, but a bit more rigorous due to the threshold requirement. This is similar to Parallelepiped Classification, with the exception that spectral regions are defined by user, and can be assigned independently for the upper and lower bounds.


texts: Equivalent of GetFeature POST request using GET? GeoServer

texts:
Right now I'm making GetFeature requests using WFS with http POST method. Due to project requirements I need to make these requests using GET method.

I'm facing serious issues traslating the POST body to a GET url query string.

This is the body:

<?xml version="1.0" encoding="UTF-8"?>
<GetFeature xmlns="http://www.opengis.net/wfs" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" service="WFS" version="1.1.0" outputFormat="application/json" xsi:schemaLocation="http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.1.0/wfs.xsd">
   <Query typeName="TRAC:v_expedientesmapa" srsName="EPSG:25830">
      <Filter xmlns="http://www.opengis.net/ogc">
         <Contains>
            <PropertyName>geom</PropertyName>
            <Point xmlns="http://www.opengis.net/gml" srsName="EPSG:25830">
               <pos srsDimension="2">574026.178090231 4753824.833116741</pos>
            </Point>
         </Contains>
      </Filter>
   </Query>
</GetFeature>

And this is what I have achieved so far:

  http://foo.foo?service=wfs&version=1.3.0&request=GetFeature&typeNames=TRAC:v_expedientesmapa

Don't know how to set the contains filter on the URL. I don't know how to add the srs either.