texts:
I am using "Point Sampling Tool" in QGIS 3.2 to obtain values from a raster file to a point layer.

After selecting the layer with the points, the values to get values from, and defining the Output point vector layer, I run the process and I get this error:

I already have all layers with the same projection, and the point layer is a shapefile. I think the problem could be during the process of saving the new shapefile because it seems that the values were obtained (see bellow in the image, "Status: Processing point 180 out of 180").

Update: I tried the same procedure but this time with QGIS 2.18, and I get a similar error: "The new layer seems to be created, but is invalid. It won't be loaded"

(which ends up being the same problem as this one: QGIS Point Sampling Tool: Invalid layer)

Does anyone know how to solve this error?


texts: Classification algorithm for landsat classification with training sites biased towards only few tiles

texts:
I would like to classify several LANDSAT tiles (18 of them), but only have training sites for 3 of the tiles.

Should I merge/mosaïc all the tiles before running a supervised
classification in this case?
Or should I run an unsupervised
classification on the tiles since I do not have training sites for
all of them?


texts: Mercator projection that shrinks toward pole

texts:
I've noticed that lots of map viewer (especially those nautical chart viewer) will shrink a mercator projected chart while panning towards poles.

Is this some compensation for mercator distortion?

Is this a dedicated projection?

How to implement that?


texts: Iterating over Twitter data using TwitterSearch creating array in Arcpy for shapefile

texts:
Below is the section of code I'm having trouble with. How do I designate the x and y coordinates from the (lat, lng) variable in the twittersearch for loop in the top part? These coordinates will be inserted into a feature class/shapefile to then be loaded into ArcGIS as point data.

for tweet in ts.search_tweets_iterable(tso):
        print( '@%s tweeted: %s' % ( tweet['user']['screen_name'], tweet['text'] ) ).encode('ascii', 'ignore')
        #ignore unicode here, show error on map??
        if tweet['place'] is not None:
            (lat, lng) = geo(tweet['place']['full_name'])
            print 'And they said it from (' + str(lat) +', ' + str(lng) + ')'
        else:
            print "And their place wasn't specified..."

    with arcpy.da.InsertCursor("D:/UW/Winter_2016/501/w7_twitter/twitterSearchFC.shp", ["SHAPE@XY"]) as cursor:
        x = float(lat).xcor()
        y = float(lng).ycor()
        varPoint = arcpy.Point(x,y)
        array.add(varPoint)
        twitterPoint = arcpy.Point(array)
        cursor.insertRow([twitterPoint])

After running the script this is what I get from Twitter, if that helps.

@12voltarts tweeted: #chihuly bowl @ Chihuly Garden and Glass https://t.co/Hf0RqE5nU3
And they said it from (47.6038321, -122.3300623)

@twanderingeater tweeted: [new post] @ChihulyGG the home of many gorgeous #Chihuly pieces in #Seattle https://t.co/OPPwvd7hQ2 https://t.co/C4HSdHZs5c
And their place wasn't specified...

My thinking is that I will need to identify each lat and lng as a new variable and call that in the cursor piece. I'm almost positive there is an easier way to do this...

also the error:

x = float(lat).xcor()
AttributeError: 'float' object has no attribute 'xcor'

from geopy import geocoders

def geo(location):
    g = geocoders.Nominatim()
    loc = g.geocode(location)
    return float(loc.latitude), float(loc.longitude)

output is a shapefile with twitter location data (tweets) spatially referenced to be able to see it in ArcGIS


texts: Using AddJoin_management gives ERROR 000840?

texts:
While i'am trying to execute this code IN IDLE, i got error :

import arcpy
import math
import arcview
from arcpy import env

arcpy.env.workspace = "D:/Users/saadiya/Documents/ArcGIS/Default.gdb"
arcpy.env.overwriteOutput = True

#Create a table of neighbors
arcpy.PolygonNeighbors_analysis("Parcelles_class_copy","Voisins_table") #Parcelles_class_copy est le résultat de la conversion de dwg en polygon

#Select each table neighbors
for i in range(1,8) :
    expression = '"src_OBJECTID"=' + str(i)
    output_table = "output_table" + str(i)
    arcpy.TableSelect_analysis("Voisins_table", output_table, expression)
    #print expression

#***************** supprimer la table Voisins_table **********************************************************************
arcpy.Delete_management("Voisins_table")

#Joint table of neighbors to feautre class
arcpy.AddJoin_management("Parcelles_class_copy","OBJECTID", output_table,"nbr_OBJECTID","KEEP_ALL")

del expression

#Select a target parcel

target_parcel = "target_parcel" + str(i)
whereclause1 = '"Parcelles_class.OBJECTID"=' + str(i)
arcpy.Select_analysis("Parcelles_class_copy",target_parcel,whereclause1)
arcpy.RemoveJoin_management("Parcelles_class_copy",output_table)

i got this error :

ERROR 000840: The value is not a Raster Catalogue Layer
ERROR 000840: The value is not a Raster Layer.
ERROR 000840: The value is not a Mosaic Layer
Failed to execute (Addjoin)


texts: How do I use a DHM in tiff format for distance and slope measurements?

texts:
I am a student and I have to work out a project about soil loss on a plot. I need basic measurements of this plot such as length, width, area and slope for my erosion calculations.

I have got my hands on a 4 GB large DHM measured by LIDAR in a TIFF file format.

My question is now, how do I preform these calculations on this file format? I guess I need to convert it so I can use it in IDRISI or ArcMap?


texts: Clipping raster with polygon for band-by-band output using ArcGIS Desktop?

texts:
How do you clip a raster image with a polygon shapefile so that the output is multiple clipped rasters, each one a band (i.e., clipped raster band 1, band 2, and so on)?

I've tried the clip tool for raster, but that doesn't output each individual band.


texts: WFS Spatial Operator Query - Crosses AND Within

texts:
I'm trying to query my wfs service with a polygon. I want the service to return both features that are inside the polygon and partially inside the polygon. However, when I use the "Within" spatial operator, I only get features that are completely inside. When I use the "Crosses" operator, I only get the features that are partially inside the polygon. I have tried using And/Or to query for features that are either one or the other, but WFS does not support multiple spatial operators in a query. It gives me this error:

  "Unsupported filter - filter has more than one spatial operator."

Does anyone have an idea of how to get around this without sending two separate ajax requests to the service?

Here's an example that returns data with the "within" filter.

var xml = '<ogc:Filter xmlns="http://www.opengis.net/ogc">\n'
             +      '<ogc:Within>\n'
             +        '<ogc:PropertyName>Shape</ogc:PropertyName>\n'
             +        '<gml:Polygon>\n'
             +          '<gml:outerBoundaryIs>\n'
             +            '<gml:LinearRing>\n'
             +              '<gml:coordinates>' + queryObject.geometryString + '</gml:coordinates>\n'
             +            '</gml:LinearRing>\n'
             +          '</gml:outerBoundaryIs>\n'
             +        '</gml:Polygon>\n'
             +      '</ogc:Within>\n'

             +    '</ogc:Filter>\n'

Here's the capabilites of the wfs service:

<ogc:SpatialOperators>
<ogc:SpatialOperator name="BBOX"/>
<ogc:SpatialOperator name="Equals"/>
<ogc:SpatialOperator name="Disjoint"/>
<ogc:SpatialOperator name="Intersects"/>
<ogc:SpatialOperator name="Crosses"/>
<ogc:SpatialOperator name="Touches"/>
<ogc:SpatialOperator name="Within"/>
<ogc:SpatialOperator name="Contains"/>
<ogc:SpatialOperator name="Overlaps"/>
</ogc:SpatialOperators>


texts: CalculateField by dividing with a variable that holds a number

texts:
I'm trying to calculate a field (using idle python to write my scipt and working with ArcMap 10.3) by dividing two fields then dividing by a variable that holds a number, which the number was calculated by dividing 2 sums.

So essentially I'm trying to divide two ratios but I'm having trouble going about that. Basically the equation is ( a / b ) / ( sum(x) / sum(y)).

I have the variable set to hold the number of the second ratio of (sum(x) / sum(y)) but can't figure out how to incorporate the variable into the calculate field expression or how to set up a double division in the expression. As you can see I used TableToNumPyArray to figure out the sum ratio.

field = arcpy.da.TableToNumPyArray (fc, ["E_Total", "E__Black"], skip_nulls=True)
LQC = 1.0 *(field["E__Black"].sum()) / (field["E_Total"].sum())
#Fields are in long format, 1.0* forces it to be in float format and this section of the code works!

field = arcpy.da.TableToNumPyArray (fc, ["E_Total", "E__Black"], skip_nulls=True)
Lratio = 1.0 *(field["E__Black"]) / (field["E_Total"])
expression = '(!'+Lratio+'!/!'+LQC+'!)'
arcpy.CalculateField_management(lyr, "LQ_Black", expression, "PYTHON")


texts: Deselecting previous selection in OpenLayers?

texts:
Previously we did when you click on a geometry on the map to color in another color and zoom on geometry. But now you should deselect the previous selection in the original color in which the map was loaded and selected recently selected object.

 map.on('click', function(e){
         console.log(e);
        //coordinate = evt.coordinate;

        var feature = map.forEachFeatureAtPixel(e.pixel, function(feature) { return feature; });

        var coordinate = e.coordinate;
        var hdms = ol.coordinate.toStringHDMS(ol.proj.transform(coordinate, 'EPSG:31700', 'EPSG:3857'));

        _myStroke = new ol.style.Stroke({
           color : 'yellow',
           width : 4
        });

        _myFill = new ol.style.Fill({
           color: 'rgba(255,0,0,1.0)'
        });

        selected_polygon_style = new ol.style.Style({
           stroke : _myStroke,
           fill : _myFill
         });

       feature.setStyle (selected_polygon_style);
       map.getView().fit(feature.getGeometry(), { duration: 2000 }, {"maxZoom":11} );

    })


texts: How to apply a Cloud Mask in Google Earth Engine? (Landsat 5 TM 8-Day NDVI Composite)

texts:
I need to download ndvi values for some pixels for 2000-2001 using the Landsat 5 TM 8-Day NDVI Composite, but since there is no quality band as in mod13q1 there's no way to filter by quality. Therefore, I want to make a mask using an image collection with information about clouds or shadows so I can mask those pixels with high o low values.
I have the script to download the ndvi values but I don't know how to apply the mask neither how to create the mask..

any experience with that?

Code:

var puntos = ee.FeatureCollection('ft:11qUO03jhF4MiECWgd5Rev536WFERNiqWf5qpqXwb');

var landsat = ee.ImageCollection('LANDSAT/LT5_L1T_8DAY_NDVI').select('NDVI')
                                     .filterDate('2000-01-01', '2000-02-28');

var f_extraer_puntos = function(i){return i.reduceRegions({collection: puntos, reducer: ee.Reducer.mean(), scale: 30, })};

// #Extract the values by running reduceRegions over each image in the image collection.
var valores_ = landsat.map(f_extraer_puntos);

var ext_puntos = ee.FeatureCollection(valores_.flatten());
// #Turn the result into a feature collection

Export.table.toDrive({ collection: ext_puntos, description:'LT5_NDVI_puntos', fileFormat: 'CSV' });


texts: Elevation profile and chemical concentrations

texts:
I have a IDW layer showing the pattern of a chemical concentration. I have also created an elevation profile across the same area.

Is there any way to include the IDW measurements with my elevation profile?

Using Arc 10.4.1, I can create an elevation profile like this one below, can I 'fill in' the area below the line with the IDW layer?


texts: QGIS Raster layer query elevation at point

texts:
In QGIS 3.14 I try to get elevation data from denmarks elevation model using "kortforsyrningen" plugin. The elevation model (layer name = 'DHM/Terræn') is added as a raster layer (QgsRasterLayer with qgis._core.QgsSingleBandColorDataRenderer)
I use this code (copied from here) to extract elevation at a point in Denmark:
from qgis.core import QgsProject
dhmLayer = QgsProject.instance().mapLayersByName('DHM/Terræn')[0]
x = 724501
y = 6175727
val, res = dhmLayer.dataProvider().sample(QgsPointXY(x, y), 0)

but all I get is val == nan and res == False.
I also tested the other way:
ident = dhmLayer.dataProvider().identify(QgsPointXY(x, y), QgsRaster.IdentifyFormatValue)

With similar results (ident.result() == {}).
ident.error() yields

<QgsError: WMS provider Format not supported>

Does anyone have any hint to how to solve this?


texts: Restricting edits for specific fields in a feature class for particular users

texts:
How to restrict edits for specific fields in a feature class for particular users in ArcSDE.

For example, there is a feature class called “Tower”. It has 20 Fields. Out of 20 fields, 10 are updated by interface.
So these “10” attributes should not be editable by “Editor” group of desktop GIS users. They may edit the other attributes in this feature class.


texts: Clipping TIFF raster image using bounding box with GDAL in Java

texts:
I would like to know how can I clip a raster image using a bounding box with GDAL in Java.

I can do it with the binary file gdal_translate:

gdal_translate -projwin -75.3 5.5 -73.5 3.7 -of GTiff original.tif new.tif


texts: mean coordinates and measurement of distances

texts: I am trying to find a location that minimises the total distance to a series of points, and compare it with the distance to an existing location. I used the QGIS tool mean coordinate(s) to calculate the centroid and then calculated a distance matrix to the centroid and the existing location. I expected the sum of distances to the centroid to be shorter than the sum of distance to the existing location but this is not the case. I've tried with a WGS and a projected coordinate system (I guess results are only valid for the second) and in both cases the centroid did not minimise distance. Is my procedure flawed ?


texts: QGIS display map on a certain lat/lon range

texts: I really want to make a screen capture of New York city (have min_lat, min_lon, max_lat, max_lon). Later I will use that picture as background and plot other information in python. Is there a way to let QGIS display the certain range?


texts: Python GeoJSON - TypeError: string indices must be integers

texts:
I am trying to iterate JSON objects and append selected features to new GeoJSON using geojson Python library.  When hardcoded the script works fine. But when append method is used in a for loop I get this error.

  TypeError: string indices must be integers

Why does this work (polygon variable hardcoded and taken from print statement in the loop)...

polygon = Polygon([[[-112.30701948248064, 33.65181795979867], [-112.30701798409073, 33.651669879089], ...[-112.30701948248064, 33.65181795979867]]])

features = []
features.append(Feature(geometry=polygon, properties={"country": "Spain"}))
feature_collection = FeatureCollection(features)

But this does not...

with open(subs) as f:
   data = json.load(f)

for feature in data['features']:
   if feature['properties']['PUC'] == '0261':
      polygon = 'Polygon(' + str(feature['geometry']['coordinates']) + ')'
      features.append(Feature(geometry=polygon, properties={"country": "Spain"}))


texts:
Can i add two sources to Vector layer, i want regular feature source, and geojson source?

Maybe can I do something like:

var layer = new ol.layer.Vector({source: VectorSource, GeoJsonSource});

My reason why i need to sources to one Vector layer because i have two different sources to load the features geojson and regular lat long.


texts:
I got this 42 tiles that I'm trying to show each tiles layer name on top of it like a label.

So I can know each tiles name, any Ideas on how to do that?
I have searched the web and found nothing


texts: Converting raster projection info in .TAB file to proper .prj4 file?

texts:
I have rasters in ecw format and they are accompanied by tab files which contain their crs information. If I try to translate these rasters to a different format with raster --> conversion --> translate tool (gdal_translate), I get following error message and the process won't start:

  ERROR 1: No PROJ.4 translation for destination SRS, coordinate transformation initialization has failed.

  ERROR 1: No PROJ.4 translation for source SRS, coordinate transformation initialization has failed.

So assumably the process can't read TAB files for CRS info.

How to proceed?

gdalinfo for input raster (ecw).

Contents of TAB file (I think you can't run gdalinfo for a TAB file).

  !table
  !version 300
  !charset WindowsLatin1

  Definition Table   File "raster_bw.ecw"   Type "RASTER"
  (23425461.980317321,6759132.2255573412) (0,0) Label "Pt 1",
  (23471612.780317321,6712981.4255573414) (307672,307672) Label "Pt 2",
  (23471612.780317321,6759132.2255573412) (307672,0) Label "Pt 3",
  (23425461.980317321,6712981.4255573414) (0,307672) Label "Pt 4"
  CoordSys Earth Projection 24, 115, "m", 23, 0, 1, 23500000, 0   Units
  "m" begin_metadata "\IsReadOnly" = "FALSE" "\MapInfo" = ""
  "\MapInfo\TableID" = "3ff0e4af-5e08-48c9-86ac-b8e63ca4cf00"
  end_metadata

My gdal_data file location is in: C:\OSGeo4W64\share\gdal? Though, running echo %GDAL_DATA% in cmd gives me a file location in my PostgreSQL folder which seems to have an identical Gdal-data folder.


texts:
I have created a database from QGIS Cloud free account. I have added a Post GIS layer into QGIS desktop. i have installed versioning plugin and added a group. After Selecting a Group when Hostorize the following error message pop out. how to fix this?
An error has occurred while executing Python code:
AttributeError: 'NoneType' object has no attribute 'write'
Traceback (most recent call last):
  File "C:/Users/parul Jain/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\qgis_versioning\versioningDB\utils.py", line 94, in execute
    self.cur.execute(sql)
psycopg2.ProgrammingError: relation "versioning_constraints" already exists

During handling of the above exception, another exception occurred:
Traceback (most recent call last):
  File "C:/Users/parul Jain/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\qgis_versioning\plugin.py", line 940, in historize
    versioning.historize(self.pg_conn_info(), schema)
  File "C:/Users/parul Jain/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\qgis_versioning\versioningDB\versioning.py", line 64, in historize
    pcur.execute(sql.format(schema=schema))
  File "C:/Users/parul Jain/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\qgis_versioning\versioningDB\utils.py", line 97, in execute
    sys.stderr.write(traceback.format_exc())
AttributeError: 'NoneType' object has no attribute 'write'

Python version: 3.7.0 (v3.7.0:1bf9cc5093, Jun 27 2018, 04:59:51) [MSC v.1914 64 bit (AMD64)]
QGIS version: 3.6.0-Noosa Noosa, 58734527ab

Python Path:
C:/PROGRA~1/QGIS 3.6/apps/qgis/./python
C:/Users/parul Jain/AppData/Roaming/QGIS/QGIS3\profiles\default/python
C:/Users/parul Jain/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins
C:/PROGRA~1/QGIS 3.6/apps/qgis/./python/plugins
C:\Program Files\QGIS 3.6\bin\python37.zip
C:\PROGRA~1\QGIS 3.6\apps\Python37\DLLs
C:\PROGRA~1\QGIS 3.6\apps\Python37\lib
C:\Program Files\QGIS 3.6\bin
C:\Users\parul Jain\AppData\Roaming\Python\Python37\site-packages
C:\PROGRA~1\QGIS 3.6\apps\Python37
C:\PROGRA~1\QGIS 3.6\apps\Python37\lib\site-packages
C:\PROGRA~1\QGIS 3.6\apps\Python37\lib\site-packages\win32
C:\PROGRA~1\QGIS 3.6\apps\Python37\lib\site-packages\win32\lib
C:\PROGRA~1\QGIS 3.6\apps\Python37\lib\site-packages\Pythonwin
C:/Users/parul Jain/AppData/Roaming/QGIS/QGIS3\profiles\default/python
C:\Users\parul Jain\AppData\Roaming\QGIS\QGIS3\profiles\default\python\plugins\mmqgis/forms
C:/Users/parul Jain/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\qgis_versioning
.
C:\Users\parul Jain\AppData\Roaming\QGIS\QGIS3\profiles\default\python\plugins\QRealTime


texts: UK regional shape map with Scotland, Wales and Northern Ireland

texts:
We have obtained boundary line data for England from the Office for National Statistics portal, which displays English regions like:

We want boundary line data for regions within the United Kingdom. So the English regions and also the following countries:

Scotland,
Wales,
Northern Ireland.

So to display something like this:

I would like help with where to source appropriate boundary line data/shape maps?
Or the Best method of combining existing files into one map?


texts: NameError: global name 'QObject' is not defined - When I add a second Toolbar Button and Dialog to a QGIS plugin

texts:
I have a problem  when I add a second Toolbar Button and Dialog to a QGIS plugin, I load my plugin in "Manage an installing plugins..." but this error appears

This is a part of my .py

What is the problem?

I use this link: How can I add a second Toolbar Button and Dialog to a QGIS plugin built with Plugin Builder?


texts: Using arcpy to find geometry within another feature

texts:
In Arcpy, I would like to start with a point and then using its geometry method to find out if it is inside my polygon feature (Urban Areas) deeper than the buffer distance.

This is a follow up to this post:  Zonal Analysis using geometries

distances=[100,200,300,400]
for i in range(4):
    buf=point.buffer(distances[i])
    if buf.within(UrbanAreaFC):  # Error here
        arcpy.AddMessage("Inside")
    else:
        arcpy.AddMessage("Outside")

The error is:  type 'exceptions.ValueError'  and the path to the geodatabase where Urban Areas is found is given.

What is that error and how do I fix it?


texts:
Is it possible to open a large image (20 to 30 GB file size) in ArcMap?

I wanted to know this in advance before doing a project in which I would use UAV to capture very high resolution images and then produce a composite/mosaic final single image.


texts: Looking for a Table of Contents sample for the ArcGIS Server JavaScript API

texts:
Where can I find an example of a Table of Contents control to plug into the ArcGIS Server JavaScript API 2.3? Ideally it should support feature/dynamic/tiled layers; show symbology; and allow grouped layers. Transparency and zoom-to functionality would be ideal.

The Esri samples page is strangely lacking in any TOC examples - the nearest example is the extremely rudimentary Toggle Layer Visibility.

Answers to Kirk's question show a few templates but they don't contain a decent TOC. The ROK sample is close but doesn't show symbology or support feature layers (as far as I can see).

The TOC on the ArcGIS.com My Map viewer is almost exactly what I need. Has anyone managed to build something similar, or get this TOC working locally?


texts: How to register table with geodatabase with a command line?

texts:
Is there a way to register with geodatabase a table stored in an Oracle 11g database, without using ArcMap or Python script?
I'd like to do it with a command line like when registering with ArcSDE (sdelayer -o register ...)


texts: How to compare the results of the random forest algorithm with the other SDMs?

texts:
I'm using the random forests algorithm in species distribution context. I have presence-absence as a dependent variable, so I decided to work with classification. At this point, I recovered the probability of presence this way :
predict(RF_MODEL, test_, type='prob')[,2]

If I want to compare the results of this model with other models "glm", "gam" and "brt", can we do this or we have to work with regression to have the estimation of occurrence probabilities as the other models?


texts: When I try to open a project in QGIS unexpectedly quits immediately without any error messages. Now, I open a new project with the same shapefile, when I try to edit or make filters QGIS closes promptly. What should I do to recover the project? Should I reinstall the program?


texts: I want create a crystal report in gis with writing code

texts:
I am working in ArcGIS 9.3.1 and I try to write code for creating a crystal report and work with it in VBA, but I recieved an error when the program arrived at this line of code:

Dim pProj As CRAXDDRT.Application
Dim pReport As CRAXDDRT.Report

The error is "User defined type not defined". Why does the program give this error?
Please help me with this problem.


texts: How to import ECW into my postgis database with raster2pgsql (windows)

texts:
I'm working on windows 7 64 bits, postgis 2.01. I have to import ECW files into my postgis database. Unfortunatly, my postgis seems to has been compiled without gdal ecw support.When i do a raster2pgsql -G, ECW does not appears.

I have tried to follow instructions erected in this following post :
How to add ECW support to GDAL bundled with PostGIS Windows

So I paste the raster2pgsql.exe into my C:\Program Files\PostgreSQL\9.1\gdal-data folder. And i paste gdal19.dll from C:\Program Files (x86)\Quantum GIS Lisboa\bin toward C:\Program Files\PostgreSQL\9.1\bin. I rename the dll like libgdal-1.dll.

But it does not work and i can't import ecw with raster2pgsql.I have an error message saying : application error.For information/. Moreover, when i do a SELECT postgis_full_version();

"POSTGIS="2.0.1 r9979" GEOS="3.3.5-CAPI-1.7.5" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.9.1, released 2012/05/15"

And my qgis 's gdal has been also compiled with gdal 1.9.1.

Moreover, for information, my GDAL_DATA (environment system variable) pointed to C:\Program Files\PostgreSQL\9.1\gdal-data

Have you got any idea? How to use ECW with postgis? I have to use this ECW format. My client give to me data with this format and i did not convert raw data in ECW with a proprietary software..

Thank you!


texts: Use less temporary space in DWithin join

texts:
I only have 12GB of free space. Is there a way to create this table without creating a bunch of temporary files? It runs out of space and never completes. I want to know how many cities are within 80km of each city.

SELECT a.k, count(a) x
FROM     cities AS a
JOIN     cities AS b
  ON     ST_DWithin(a.geog, b.geog, 80000)
 AND     a.k != b.k
WHERE    a.l = true
group by a.k
;

QUERY PLAN

-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
--
 Finalize GroupAggregate  (cost=1282860.42..1285716.94 rows=24280 width=15)
   Group Key: a.k
   ->  Gather Merge  (cost=1282860.42..1285372.97 rows=20234 width=15)
         Workers Planned: 2
         ->  Partial GroupAggregate  (cost=1281860.39..1282037.44 rows=10117 width=15)
               Group Key: a.k
               ->  Sort  (cost=1281860.39..1281885.69 rows=10117 width=348)
                     Sort Key: a.k
                     ->  Nested Loop  (cost=0.41..1281187.39 rows=10117 width=348)
                           ->  Parallel Seq Scan on cities a  (cost=0.00..11417.42 rows=111968 width=380)
                                 Filter: l
                           ->  Index Scan using citiesn_geog_idx on cities b  (cost=0.41..11.33 rows=1 width=39)
                                 Index Cond: (geog && _st_expand(a.geog, '80000'::double precision))
                                 Filter: ((a.k <> k) AND (a.geog && _st_expand(geog, '80000'::double precision)) AND _st_dwithin(a.geog, geog, '80000'::double precision, true)
)
(14 rows)

when I only limit it to 1000 it is really fast

Limit  (cost=0.83..127479.34 rows=1000 width=15) (actual time=0.858..440.869 rows=1000 loops=1)
   ->  GroupAggregate  (cost=0.83..3095179.08 rows=24280 width=15) (actual time=0.856..440.550 rows=1000 loops=1)
         Group Key: a.k
         ->  Nested Loop  (cost=0.83..3094814.88 rows=24280 width=348) (actual time=0.711..423.748 rows=42526 loops=1)
               ->  Index Scan using citiesk on cities a  (cost=0.42..47380.58 rows=268722 width=380) (actual time=0.027..4.349 rows=1043 loops=1)
                     Filter: l
                     Rows Removed by Filter: 79
               ->  Index Scan using citiesn_geog_idx on cities b  (cost=0.41..11.33 rows=1 width=39) (actual time=0.083..0.389 rows=41 loops=1043)
                     Index Cond: (geog && _st_expand(a.geog, '80000'::double precision))
                     Filter: ((a.k <> k) AND (a.geog && _st_expand(geog, '80000'::double precision)) AND _st_dwithin(a.geog, geog, '80000'::double precision, true))
                     Rows Removed by Filter: 11
 Planning Time: 0.480 ms
 Execution Time: 448.505 ms
(13 rows)

after vacuum analyze and with count(a.k)

GroupAggregate  (cost=3058755.01..3059179.53 rows=24258 width=15)
   Group Key: a.k
   ->  Sort  (cost=3058755.01..3058815.66 rows=24258 width=7)
         Sort Key: a.k
         ->  Nested Loop  (cost=0.41..3056988.28 rows=24258 width=7)
               ->  Seq Scan on cities a  (cost=0.00..12999.33 rows=268413 width=39)
                     Filter: l
               ->  Index Scan using citiesn_geog_idx on cities b  (cost=0.41..11.33 rows=1 width=39)
                     Index Cond: (geog && _st_expand(a.geog, '80000'::double precision))
                     Filter: ((a.k <> k) AND (a.geog && _st_expand(geog, '80000'::double precision)) AND
 _st_dwithin(a.geog, geog, '80000'::double precision, true))
(10 rows)


texts:
Im looking for a way to add a marker on a position that is not in a database using cartodb.js.

Something similar to this in Mapbox would be great:

L.mapbox.featureLayer({
    // this feature is in the GeoJSON format: see geojson.org
    // for the full specification
    type: 'Feature',
    geometry: {
        type: 'Point',
        // coordinates here are in longitude, latitude order because
        // x, y is the standard for GeoJSON and many formats
        coordinates: [
          -77.03221142292,
          38.913371603574
        ]
    },


texts: HCS resolution merge (artifacts and dead pixels)

texts:
I try to use ATCOR and HCS resolution merge to preprocess urban Worldview-2 data. Applying the HCS algorithm without initial atmospheric/radiometric correction works just fine, but if I apply ATCOR beforehand, results get severely damaged. This happens:

Also, after using ATCOR the data also looks alright, and I can't seem to find the source for these distortions. If needed, I can deliver further information on the input data, histograms or whatever might help.


texts:
Given a peninsula stretching north south and the two coast lines connecting in the south, I want to find out which is the place, or rather line, which is furthest away from both coastlines on the peninsula? What I am looking for is the red line in the below image. The black and pink lines are the eastern and western coastlines.

One way could be to generate a polygon and do a negative buffer. An other one could be to add the longitudinal values, and divide by two. The latter one seems more rational - still it demands some programming and coordinate handling.

But perhaps there is a third way? Any suggestions? ArcGIS and QGIS are both relevant tools of the trade here.


texts: Determine XY offset in meters between two lat/long points (EPSG:4326), and accurately recreating point locations using XY offset

texts:
I have two tables in EPSG:4326, one for home locations and one for site locations. I need to be able to generate the XY offset in meters between the home location, and the site locations. I then need to be able to accurately regenerate the site locations, based on the XY offset from the home location. Currently when I use ST_Project to recreate a site location by using the specified XY offset from the home location. The results are slightly off in the Y direction. Generating lines using the same logic used to calculate the offset, is accurate. My environment is PostgreSQL 9.4.1 PostGIS 2.1.

The SQL script is my test data for this issue.  The home and sites test tables are the two tables stored in EPSG:4326. The sites table contains a delta_x, and delta_y field representing the X & Y offset in meters. The offset_lines view contains lines generated using similar logic as was used to calculate the offset. These lines appear accurate and terminate at the exact location as the site. The points_using_offsets view is a view of the site locations generated using the offsets from the home location and the ST_Project command. These locations are all slightly shifted in the Y direction.

CREATE TABLE home
(
  gid integer,
  geog geography(Point,4326)
);

CREATE TABLE sites
(
  gid integer,
  geog geography(Point,4326),
  delta_x float,
  delta_y float
);

insert into home (gid, geog) values (1, ST_MakePoint(-78.9570, 41.6241)::geography);

insert into sites (gid, geog) values (1,ST_MakePoint(-79.1133, 41.7690)::geography);
insert into sites (gid, geog) values (2,ST_MakePoint(-78.7726, 41.7801)::geography);
insert into sites (gid, geog) values (3,ST_MakePoint(-78.7753, 41.4041)::geography);
insert into sites (gid, geog) values (4,ST_MakePoint(-79.4226, 41.5032)::geography);

UPDATE sites AS st set delta_x = foo.delta_x, delta_y = foo.delta_y
FROM
(
 SELECT
  CASE
    WHEN ST_X(s.geog::geometry) < ST_X(h.geog::geometry) THEN
      ST_Distance_Spheroid(ST_MakePoint(ST_X(h.geog::geometry), ST_Y(h.geog::geometry) ),ST_MakePoint(ST_X(s.geog::geometry), ST_Y(h.geog::geometry)),'SPHEROID["WGS 84",6378137,298.257223563]')*-1
    ELSE
      ST_Distance_Spheroid(ST_MakePoint(ST_X(h.geog::geometry), ST_Y(h.geog::geometry)),ST_MakePoint(ST_X(s.geog::geometry), ST_Y(h.geog::geometry)),'SPHEROID["WGS 84",6378137,298.257223563]')
    END AS delta_x,
  CASE
    WHEN ST_Y(s.geog::geometry) < ST_Y(h.geog::geometry) THEN
      ST_Distance_Spheroid(ST_MakePoint(ST_X(s.geog::geometry), ST_Y(s.geog::geometry)),ST_MakePoint(ST_X(s.geog::geometry), ST_Y(h.geog::geometry)),'SPHEROID["WGS 84",6378137,298.257223563]')*-1
    ELSE
      ST_Distance_Spheroid(ST_MakePoint(ST_X(s.geog::geometry), ST_Y(s.geog::geometry)),ST_MakePoint(ST_X(s.geog::geometry), ST_Y(h.geog::geometry)),'SPHEROID["WGS 84",6378137,298.257223563]')
    END AS delta_y,
    s.gid
FROM
  home as h,
  sites as s) as foo
where st.gid = foo.gid;

CREATE OR REPLACE VIEW offset_lines AS
SELECT
   s.gid as gid,
   ST_SetSRID(ST_MakeLine(ARRAY[ST_MakePoint(ST_X(h.geog::geometry), ST_Y(h.geog::geometry)),ST_MakePoint(ST_X(s.geog::geometry), ST_Y(h.geog::geometry)),ST_MakePoint(ST_X(s.geog::geometry), ST_Y(s.geog::geometry))]),4326) as geom
FROM
  home as h,
  sites as s;

CREATE OR REPLACE VIEW points_using_offsets AS
SELECT
CASE
  WHEN delta_x > 0 and delta_y > 0 THEN
      ST_Project(ST_Project(h.geog, abs(s.delta_x), radians(90)),abs(s.delta_y),radians(0))
  WHEN delta_x > 0 and delta_y < 0 THEN
      ST_Project(ST_Project(h.geog, abs(s.delta_x), radians(90)),abs(s.delta_y),radians(180))
  WHEN delta_x < 0 and delta_y < 0 THEN
      ST_Project(ST_Project(h.geog, abs(s.delta_x), radians(270)),abs(s.delta_y),radians(180))
  WHEN delta_x < 0 and delta_y > 0 THEN
      ST_Project(ST_Project(h.geog, abs(s.delta_x), radians(270)),abs(s.delta_y),radians(0))
  END as geom,
  s.gid as gid
FROM
home as h,
sites as s;


texts: Error while rendering a layer from postgresql in mapserver

texts:
I got a point layer into QGIS from postgresql database. After this I used a plugin to export it to a map file.After exporting the following is the starting of the map file.

 NAME "QGIS-MAP"
 # Map image size
 EXTENT 59.396660 7.316000 105.983338 35.204001
 Imagecolor 255 255 255
 SIZE 100 100
 UNITS meters
 Imagetype jpeg

 FONTSET 'C:\ms4w\Apache\htdocs\fontset.txt'

PROJECTION
'proj=longlat'
'datum=WGS84'
'no_defs'

Now when I change the size parameter to 600 400 and the units to dd. It gives me the following error while rendering the layer.

msDrawMap(): Image handling error. Failed to draw layer named 'growth_1'. msPostGISLayerWhichShapes(): Query error. Error executing query: ERROR: column "%growth_1_id%" does not exist LINE 1: ...ary(ST_Force_2D("the_geom"),'NDR'),'hex') as geom,"%growth_1... ^

First of all I am not sure why is it taking units as meter.But if I dont change it to dd , the legends come up but the layer is not rendered. Also the extent which i get is in the range of 100's and  not the one which I have defined.

Could someone tell me where is this thing going wrong.


texts:
I have this service: http://sdi.gdos.gov.pl/wfs?SERVICE=WFS&REQUEST=GetCapabilities&VERSION=2.0.0 in version 2.0.0.

When I do a describe feature request for a certain feature type, the key to select the relevant feature type in the KVP should be TYPENAMES. TYPENAME is used for earlier versions.

For example: http://sdi.gdos.gov.pl/wfs?SERVICE=WFS&VERSION=2.0.0&REQUEST=DescribeFeatureType&TYPENAMES=GDOS:SpecjalneObszaryOchrony should work, but it does not. It is not selecting only the feature type GDOS:SpecjalneObszaryOchrony.

But if I used http://sdi.gdos.gov.pl/wfs?SERVICE=WFS&VERSION=2.0.0&REQUEST=DescribeFeatureType&TYPENAME=GDOS:SpecjalneObszaryOchrony with TYPENAME for earlier versions, then it works.

Why TYPENAMES does not work? I added the version in the KVP, so it should be fine!


texts:
I am moving some applications written in VBA that have been working in ArcGIS 10.1 to ArcGIS 10.2.2.  I have the ArcGIS VBA Compatibility patch installed, in addition I have gotten several applications to work.  Eventually this VBA code will be rewritten in the C# environment, however at the moment I am just trying to get the code to work.

In the following example the code gets down to the

Set pPoint = pActiveView.ScreenDisplay.DisplayTransformation.ToMapPoint(X, Y)

line and errors out with
runtime error 13 Type Mismatch

I believe it is some sort of reference error, however I cannot figure out where it is hanging.

Private Sub pixel3x3picker_MouseDown(ByVal Button As Long, ByVal Shift As Long, ByVal X As Long, ByVal Y As Long)
'Code to get a 3 x 3 sample average of pixels similar to the sample tool in photoshop
Dim pMxDoc As IMxDocument
Set pMxDoc = ThisDocument

Dim pActiveView As IActiveView
Set pActiveView = pMxDoc.FocusMap
Dim pPoint As IPoint
Set pPoint = pActiveView.ScreenDisplay.DisplayTransformation.ToMapPoint(X, Y)

I got it figured out, I had to mess with the references in VBA, I turned on the correct one, however I do not know which one was the "correct one".

VBA References


texts:
My colleagues and I are using different versions of QGIS (2.8.2, 2.86 & 2.14) on windows 8.1 and Windows 10.

When we open existing projects we get a "Handle Bad Layers" box coming up. We have each tried to reset the datasource link and when we do this we get a message saying "The specified network name is no longer available."

When we copy & paste the shapefiles onto our C drive we can open them in the project. After doing this, the original file path then seems to work but is lost again when the project is closed down & reopened.

We have referred this to our company IT consultant who has advise it is a QGIS problem and not a server problem.

I also have projects & shapefiles on Dropbox and these work fine. Can anyone clarify if this is a QGIS problem or a server problem?


texts: Curvature computation for different scales

texts:
I would like to assess terrain curvature, i.e. the 2nd derivative of a DEM, on different scales.

The most common approach to curvature seems to be the one suggested by Zevenberg (1987) where a fourth-order polynom is fitted to a 3x3 window around each grid cell, as it is implemented e.g. in ArcGis. The nice thing about it is that non-directional curvature can be computed as a convolution of the DEM with the kernel

    0.0  0.5  0.0
k = 0.5 -2.0  0.5  *  1/s^2
    0.0  0.5  0.0

with s the cell size. However, I'd like to derive curvature on larger scales, too, so that the result is less affected by small terrain features.

I can see two immediate ways to do so:

Smooth or scale the underling DEM and re-use the same kernel, or
to use a larger kernel.

To me, (1) feels pretty ragged and for (2) I'm not sure what kind of kernel to use.

Any suggestions to get me in the right way? Or is it rather a problem of defining "larger scale curvature"?


texts:
My code has to check if a value in a shapefile attribute table belongs to a certain group of landuses. Below I paste my code:

cursor = arcpy.SearchCursor(buf_intersect)
fields = ["GTYS_Branc", "landuse", "landuse_1", "landuse_12", "landuse_13", "Shape_Area"]

for row in cursor:
    book = xlrd.open_workbook("MyDirectory.xls")
    first_sheet = book.sheet_by_index(0)
    index = row.getValue(fields[0])

    # Rural

    if (row.getValue(1) == 'field') or (row.getValue(1) == 'farmland') or (row.getValue(1) == 'grass') or (row.getValue(1) == 'quarry') or (row.getValue(1) == 'paddock') or (row.getValue(1) == 'vineyard') or (row.getValue(1) == 'reservoir') or (row.getValue(1) == 'orchard') or (row.getValue(1) == 'greenfield') or (row.getValue(1) == 'allotment') or (row.getValue(1) == 'conservation') or (row.getValue(1) == 'moorland') or (row.getValue(1) == 'f'):
        col = 2

    # Suburban

    elif (condition similar to the one above)
        col = 3

And so on. When I run this code, I receive the following error:

Traceback (most recent call last):
  File "Mydirectory.py", line 32, in <module> the first if line
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\arcobjects\arcobjects.py", line 1100, in getValue
    return convertArcObjectToPythonObject(self._arc_object.GetValue(*gp_fixargs(args)))
RuntimeError: ERROR 999999: Error executing function.

What can the reason be? May it be because I have too many conditions on my if statement?


texts: How to get the bounds of a GeoJSON file in GeoPandas as a seperate file?

texts:
I am trying to generate a polygon which circumscribes all shapes in the given GeoJSON (could be either point, line, or polygon shapes).
I figured out how to do this with bounding boxes (.total_bounds). This is not an optimal solution and I am looking for a minimum area that covers all the shapes within the GeoJSON.
import shapely
bound = shapely.geometry.box(*input_json.total_bounds).buffer(buffer_dist)
bound = gpd.GeoDataFrame(gpd.GeoSeries(bound), columns=['geometry'])
bound = bound.set_crs(input_json.crs)

How can I generate a polygon which contains all shapes within the file?


texts: Render land based on relation of ways with geosever and data from OSM

texts:
I have styles for ways and polygons in a geoserver setup. However the land is not being properly rendered.
I have notice that the feature for 'Arabia' (https://www.openstreetmap.org/relation/5631846#map=12/28.6282/48.4909) is a relation of ways.
What I have done is download OSM data and import it into postgis with osm2pgsql, so this relation ends up in the lines table and my current styles for lines don't fit.
So, which SDL will render the relation of lines as an area?


texts: Raster extents incorrect for vrt when changing projection with -a_srs flag

texts:
I want to generate a .vrt in epsg:3857 from a series of raster tiles in epsg:4326, and am trying to do so in one step using gdalbuildvrt with the -a_srs flag:

gdalbuildvrt -a_srs epsg:3857 ortho_3857.vrt ortho-*.tif

The projection is changed to 3857 as desired, however the corner coordinates, center and pixel sizes are incorrect (they are still in 4326, with the geographic extents of the corner coords messed up):

Driver: VRT/Virtual Raster
Files: ortho_3857.vrt
       /Users/jeremyeastwood/Desktop/debug/ortho-0-0.tif
       /Users/jeremyeastwood/Desktop/debug/ortho-0-1.tif
       /Users/jeremyeastwood/Desktop/debug/ortho-1-0.tif
       /Users/jeremyeastwood/Desktop/debug/ortho-1-1.tif
Size is 9191, 7284
Coordinate System is:
PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],
    AUTHORITY["EPSG","3857"]]
Origin = (-122.318962274868440,37.873466153295041)
Pixel Size = (0.000000446142119,-0.000000446142119)
Corner Coordinates:
Upper Left  (-122.3189623,  37.8734662) (  0d 0' 3.96"W,  0d 0' 1.22"N)
Lower Left  (-122.3189623,  37.8702165) (  0d 0' 3.96"W,  0d 0' 1.22"N)
Upper Right (-122.3148618,  37.8734662) (  0d 0' 3.96"W,  0d 0' 1.22"N)
Lower Right (-122.3148618,  37.8702165) (  0d 0' 3.96"W,  0d 0' 1.22"N)
Center      (-122.3169120,  37.8718413) (  0d 0' 3.96"W,  0d 0' 1.22"N)
Band 1 Block=128x128 Type=Byte, ColorInterp=Red
  NoData Value=0
Band 2 Block=128x128 Type=Byte, ColorInterp=Green
  NoData Value=0
Band 3 Block=128x128 Type=Byte, ColorInterp=Blue
  NoData Value=0
Band 4 Block=128x128 Type=Byte, ColorInterp=Alpha
  NoData Value=0

I can achieve the desired result by generating the vrt in the native projection then warping the result to 3857 with gdalwarp afterwards, but then I need to keep both .vrts about for future operations as the warped vrt references the native one, which references the original geotiffs, which is pretty messy.

Any suggestions on correct usage of the -a_srs flag in gdalbuildvrt (or a better understanding of its limitations), or how to generate a single warped vrt in a different projection from the raw geotiffs?


texts: arcpy.Rename_management() returning ERROR 999999 and [Errno 13]

texts:
I have a function where I am trying to rename a geodatabase, both os.rename and arcpy.Rename_management() work in the Python window but not in the script when run as a tool through ArcCatalog.
Traceback (most recent call last):
  File "C:\PecosFieldTools\ProTools_v1.2\PhaseDeliverables\phase.py", line 562, in <module>
    arcpy.Rename_management(in_data=newGdb, out_data=phaseGdb, data_type="Workspace")
  File "c:\program files (x86)\arcgis\desktop10.6\arcpy\arcpy\management.py", line 4528, in Rename
    raise e
ExecuteError: ERROR 999999: Error executing function.
Failed to execute (Rename).

Failed to execute (Phase)

I think this may be a permission issue but I am admin running this on my own machine. When I wrap in try and except I get these error
[Errno 13] Permission denied: u'C:\\\\Projects\\\\PhaseMasters\\\\HasData\\\\1810227C\\\\Master\\\\PHASE_DELIVERABLES\\\\DB13635_AF2136641_AP_00_X04.gdb\\\\EQ_AIR_COOLED_EXCHANGER_POLY.LAPTOP-3B9ORV48.16500.16420.sr.lock'")

But when I open the file geodatabase I do not see any locks and I am deleting all my cursors. Is there a way to await the locks to be removed with Python?
Here is where I am running the tool from:


texts: Calculate a field based on values repeated of an another field

texts:
I'd like to know how can i calculate in QGIS a field (field3) based on how many times a value is repeated (field2) considering an another field (field1).
For example, the field to fill with the field calculator is the "field3":

Any suggestion using the field calculator of QGIS?


texts:
I wish to use a Grass module (v.to.db) in QGis, without the grass toolbox.
I tried with the Python console, but to no avail :

Using Sextante plugin but obviously it didn't know all the grass modules.

from sextante.core.Sextante import Sextante
Sextante.alghelp("grass:v.to.db")
--> Algorithm not found

An other test :

grass.run_command("v.to.db", map='bl@PERMANENT', layer='2', option='start', units='meters', columns='X,Y,Z')
-->Traceback (most recent call last):
File "", line 1, in module
NameError: name 'grass' is not defined

Do you know if there is an other way to use the grass module ?
I don't want to use it in the grass toolbox because I can select only one colums in "attribute field", and I would like to choose several columns.


texts: How can I create a "setback" buffer for coastal lines in QGIS?

texts: I've been going a bit crazy in trying to create a "setback" buffer on QGIS for coastal lines in a map projecting this buffer outwards towards the sea. I've managed to create a normal Bidirectional buffer, but that's not the outcome that I'm looking for.


texts: Importing CSV file (Georeferenced with Google Earth) into ArcMap?

texts:
I had a CSV file that I just georeferenced using Google Earth's import feature. The CSV file contains information about businesses in the United States, such as the business' name, address, and description.

After georeferencing the data in Google Earth, Google Earth displays several points on the globe. When I click on each point, a pop up window appears that displays all of the information for the business in that location, such as the business' name, description, etc.

I would like to add that data into ArcMap, so that I have a shapefile that includes all of the attribute data in the CSV file. I right clicked on the folder that contains all of the point data in Google Earth, clicked "Save Place As". After doing so, Google earth created a .kmz file.

How do I import that .kmz file into ArcMap? When I tried to import the data using the conversion tools menu in ArcToolbox, the points all appeared in the correct locations, but the points weren't associated with any of the data in the CSV file (e.g. the business names and descriptions).

What am I doing wrong? I used Google Earth to georeference the data because it's free.


texts: How to bring workspace layer to one combo box and once layer selected in this attribute of this layer in other combo box

texts:
I want a layers in workspace in combo box.Once i select combo box i want attribute in selected layer to in other combo box. how can i do it. I am using geoserver, openlayers, Geoext and ext js.

Can any one just help me

my combo box code is:

            var simpleCombo = Ext.create('Ext.form.field.ComboBox', {
                fieldLabel: 'Select layer Name',
                displayField: 'name',
                width: 320,
                labelWidth: 100,
                store: store3,
                queryMode: 'local',
                emptyText:'Layer Selection',
                typeAhead: true,
                listeners : {
                    'select' : function(){
                        var selValW = this.getValue();

                            txtquery(selValW);
                }}
            });


texts:
I've added a custom toolbar with some buttons to my QGIS 2 plugin. I'm wondering how to get the button toggled (pushed in-looking effect) using Python. I've gotten this to work for existing buttons like so:

iface.actionAddFeature().trigger()

I can't seem to figure out how to apply it to my own actions and buttons.


texts:
I have an issue with createLayer. The date has duplicate polygons ( countries) so I ha have included a SQL to group and sum up the the values of those duplicated countries :

cartodb.createLayer(map, {
user_name: 'migrantreport',
type: 'cartodb',
sublayers: [{
sql: "SELECT country, SUM(dec_16)dec_16 FROM origins_by_country_of_arrival_sp_copy GROUP BY country",

cartocss: "#origins_by_country_of_arrival_sp_copy{polygon-fill: #FFFFB2; polygon-opacity: 0.8;line-color: #FFF;line-width: 0.5;line-opacity: 1;}",
}]
}).addTo(map)

The SQL query  works ( I already test it externally ) but when using the code above it  doesn't work nor  gives any error...

I could group all the polygons beforehand but  I need to keep it like this to perform different queries later where the data needs to be non grouped ..

any suggestion??


texts: Getting layer number from layer name in MapInfo (vis MapBasic)

texts:
Is there a function to get the layer number from the layer name in MapBasic? I am writing a script to loop through each layer I have and apply styles and labelling. Each layer is named _#_LINE or _#_BUFFER

So I was thinking of assigning the layer number as an integer, applying the styles then incrementing the # by one so it moves onto the next layer.


texts: Altitude extraction from a DEM(.tiff) file using Python

texts:
I have a DEM(.tiff) file and I need to know the altitude of certain points(coordinates).

I opened it like this:

import tifffile as tiff
import matplotlib.pyplot as plt
import rasterio

tfile = tiff.imread('ms.tif')
tfile.shape
tiff.imshow(tfile)

Then I've got this image:

I just need some function to get altitude to my points, not the full code.


texts: QGIS Python script for zonal statistics not running

texts:
I got this Python script from this previous post. I've edited it for my project but nothing happens when I execute it. I downloaded the zonal statistics plugin, and usually when I run it, it adds a new field to the shapefile with the statistics I need, for example minimum. Since I need to do it multiple times, I was hoping to use the script. I tried the batch execute, but for GDAL it doesn't allow you to select one statistic, and Saga always gives an error somewhere.

Here's the script I use:

from PyQt4.QtGui import *
from PyQt4.QtCore import *
from qgis.core import *
from qgis.utils import iface
from qgis.analysis import QgsZonalStatistics

#specify polygon shapefile vector
polygonLayer = QgsVectorLayer('D:/all_structures_tez.shp', 'zonepolygons', "ogr")

# specify raster filename
rasterFilePath = 'D:/slow.tif'

# usage - QgsZonalStatistics (QgsVectorLayer *polygonLayer, const QString &rasterFile, const QString &attributePrefix="", int rasterBand=1)
zoneStat = QgsZonalStatistics (polygonLayer, rasterFilePath, 'pre-', 1, QgsZonalStatistics.Min)
zoneStat.calculateStatistics(None)

QGIS ver 2.18.3


texts: Gdal_Translate -Projwin giving Error (Error while crease Subset from Image)

texts:
"C:\Program Files\QGIS Las\OSGeo4W.bat" gdal_translate -of EHDR -projwin 39
10865.0 330905.0 3894865.0 335885.0 E:\F\CLP_BLD.img C:\Users\kailash\Desktop\gg
g\C_Output_LOWER_CENTRAL_VALLEY_2m_Z11N_1_DTM.bil
Input file size is 18480, 15670

ERROR 1: Error: Computed -srcwin 1.80597e+006 1.79765e+006 -8000 -2490 has negat
ive width and/or height.

Can any body help to run Gdal_tranlate -projwin command


texts: Deploying ArcGIS Web App via windows IIS, domain instead of IP address

texts:
So first off, I'm a newbie at web apps. I have put together a working ArcGIS web app that is hosted on an Amazon EC2 Windows Server 2012 via IIS. The web app is working great, but I can only access it via the IP address of the windows server "http://<ip address>/ParkLocator". I want to use a regular domain name to look more professional http://www.<domain name>/ParkLocator.

My first thought was that I needed to register a domain name so I acquired a domain name via GoDaddy. I thought I could just 'forward and mask' the domain name to forward traffic to the windows server. That worked, but now the web app is no longer mobile responsive, it only works right though the desktop.

Any advice?


texts: Python AddIn button init not working as expected

texts:
I have a toolbar with two buttons. It's a crude toggle that looks
at a text file containing a single character whose values are either 0 or 1.  When the file contents are 0, then button1 ResumeExports() should show up as depressed and disabled while button2 PauseExports() is enabled and normal.  When the contents of the file are 1, then the button2 should be active and button1 inactive--it is used as a switch to allow or disallow other processing to occur.  The toggle works, and the file gets updated.  How can I get the button status (checked/unchecked, enabled/disabled) on the interface to reflect that?  Both buttons always show enabled and unchecked.  I'm at 10.2.1.

import arcpy
import pythonaddins

class ResumeExports(object):
    """Implementation for PythonAddIns_addin.res_btn (Button)"""
    ctl_file = r"C:\Users\myusername\Documents\ArcGIS\switch.txt"
    def __init__(self):
        with open(ctl_file, 'r') as i:
            content = i.read(1)
        if content == '1':          # exports are allowed
            self.enabled = False    # button should be inactive
            self.checked = True     # button should be depressed
        if content == '0': # exports are disallowed, resume button is active
            self.enabled = True
            self.checked = False
    def onClick(self):
        with open(ctl_file, "r+") as f:
            data = f.read()
            f.seek(0)
            f.write("1")
            f.truncate()

class PauseExports(object):
    """Implementation for PythonAddIns_addin.stop_btn (Button)"""
    ctl_file = r"C:\Users\myusername\Documents\ArcGIS\switch.txt"
    def __init__(self):
        with open(ctl_file, 'r') as i:
            content = i.read(1)
        if content == '0':          # exports are paused
            self.enabled = False    # button should be inactive
            self.checked = True     # button should be depressed
        if content == '1': # exports are allowed, pause button is active
            self.enabled = True
            self.checked = False
    def onClick(self):
        ctl_file = r"C:\Users\myusername\Documents\ArcGIS\switch.txt"
        with open(ctl_file, "r+") as f:
            data = f.read()
            f.seek(0)
            f.write("0")
            f.truncate()


texts: How to add an image through URL to arcgis story map?

texts:
I want to add an image to my arcgis story map by Url.
I am having an error while i am uploading my image.

How to add an image to my argis story map by using URL method?


texts: What is the field_map_name in findFieldMapIndex ()?

texts:
I am using ArcPy with ArcGIS Pro and FieldMappings.
FieldMappings are fiddly and I'm struggling to understand a specific function.
I need to replace an 'old' fieldMap with a 'new' fieldMap. FieldMaps have one or more input fields, and one output field (as I understand it). Because there are more than 1 field in a map, one has to refer to the map by Index. So, we have a handy function to find it: findFieldMapIndex (field_map_name)
# https://pro.arcgis.com/en/pro-app/latest/arcpy/classes/fieldmappings.htm
# Find a field map within the field mappings by name
findFieldMapIndex (field_map_name)

Which is great, but what exactly is the field_map_name ?
for fldMap in fieldMaps:
    print (fldMap.name)

AttributeError: 'FieldMap' object has no attribute 'name'


texts: Input a KML file dynamically into MapServer without losing attributes?

texts:
Is there any possibility for MapServer (I'm running v6.2.1) to process the features / styles from a KML file dynamically?

I tried STYLEITEM "AUTO" but it turned out to work only for *.dgn and *.tab files. Actually this'd be exactly the kind of feature I'm looking for tho.

Alternatively i'm searching for a way to get my KML file converted to *.tab (MapInfo) in order to later auto-feed it to MapServer. I tried using ogr2ogr but this always results in losing all the attributes from the KML file, leaving only the coordinates with the geometry.


texts: Change Projection of a Web Mercator PCS to rotate it around center

texts:
I am using the GIS tool SharpMap together with ProjNet(4GeoAPI) as a tool to switch between CRS and transform coordinates (This information are just for context; my question is more abstract).
The app I am building uses multiple layers and I need to rotate one layer which has a Web Mercator PCS.
Currently it looks like this; a projection with the parameters
var parameters = new List<ProjectionParameter>
     {
       new ProjectionParameter("semi_major", 6378137.0),
       new ProjectionParameter("semi_minor", 6378137.0),
       new ProjectionParameter("latitude_of_origin", 0.0),
       new ProjectionParameter("central_meridian", 0.0),
       new ProjectionParameter("scale_factor", 0.9996),
       new ProjectionParameter("false_easting", 0.0),
       new ProjectionParameter("false_northing", 0.0),
    };

var projection = csFac.CreateProjection("WebMercator", "mercator_1sp", parameters);

and from there the Webmercator-PCS gets constructed:
_epsg3857 = csFac.CreateProjectedCoordinateSystem("Google Mercator",
                _epsg4326, projection, LinearUnit.Metre, new AxisInfo("East", AxisOrientationEnum.East),
                new AxisInfo("North", AxisOrientationEnum.North));

Where _epsg4326 is a Geographic Coordinate System using the WGS84 datum and the Greenwhich Prime Meridian.
To have the layer not aligned to north, I  tried to setup the projection with the azimuth parameter and an angle, adding to the parameter list:
new ProjectionParameter("azimuth",45.0)

which does not work. I read about angular parameters in the ArcGIS help. But still the coordinates were the same as in the other layers.
What is the correct way of changing the rotation? Is this the correct way? (Somehow I am doubtful) or is it more easy to apply an Affine transformation to the needed points and "fake" a rotation?


texts: MapInfo Pro 64-bit - 16.0.3 - Snap Frame to Guideline (unable to?)

texts:
I am using MapInfo Pro 64-bit v16.0.3

In the layout, I want to snap my map frame to the guideline. I've noticed that when I create a rectangular box, this snaps to the guidelines once the edge is placed close to it.

When I try to do the same thing, but by pulling the edge of the 'frame box', it won't snap to the guideline.

The reason I am doing this is so that I can align the frame box to the exact specifications of my 'side panel' box (Holds Title, Legend, Scale, North Arrow, Disclaimer etc). I want the map frame and the side panel to be the same height (to line up exactly).

I've noticed that from 32bit to 64 bit they've removed the 'frame properties' that allow you to specify the X1, X2 and Y1, Y2. They instead only allow the top-left corner and then the width and height. Only problem is, when I change both my side-panel and map frame to be the same height exactly, it doesn't look that way in the layout.

Hence, If I can snap my frame box and rectangular (side panel) box to the guidelines, then I won't have to manually adjust the line.

If anyone knows another way of doing this, then I am all ears.


texts:
I'm a bit of a QGIS noob when it comes to anything other than the basics so hoping someone has an easy solution to this. I am trying to use Open Railway Maps ( https://www.openrailwaymap.org/ ) as a basemap for a project and have had no success following multiple guides. I tried the method recommended in the wiki ( https://wiki.openstreetmap.org/wiki/OpenRailwayMap/API#Usage_in_QGIS ) and this guide https://www.xyht.com/spatial-itgis/using-openstreetmap-basemaps-qgis-3-0/ . I have attached a screenshot of the result, it seems to only be displaying a segment of the map, I feel like this is probably some projection issue.
What can I do to sort this out?


texts: QGIS not showing CSV file when trying to join layers

texts:
I am using QGIS 2.8.9 Wien, and trying to add data from a CSV file to an existing vector file.  At the "add vector join" stage the join layer drop down only shows the vector shapefile, so I find I can only able to add data from the shapefile to the CSV, not the other way round.  I would expect both files to appear on the drop down list so I could select the CSV file and go from there.

What am I doing wrong?


texts: Extracting ways from node to node using Overpass QL

texts:
Is there a way to extract ways (basically I need path from node to node in overpass)
For example:
node(1270323047)->.c;
(.c;<;>;)->.d;
way.d;

would recursively find all relations that node 1270323047 belongs to and then extract all ways.
Now my questions is, is there a way to iterate trough all ways and have start and stop conditions for example begin at node(id1) and stop at node(id2) (assuming ways are linked to previous, next)?
I wonder if this is possible at all. Because there is no link between nodes <-> ways.
UPDATE: For example if I have relation: https://www.openstreetmap.org/relation/7916536#map=11/47.2771/8.5530&layers=T,
how can I extract only from Wiedikon - Zurich HB path (from node (train_stop) - to node (train_stop)
Then I want to do this generally from any to train stops, I want a path that connects the 2.
So far I used OSRM(Open Source Routing Machine), but there are some issues if I don't pinpoint exact location of train stations so I have some crazy 180deg turns. (even if I set penalty on lua profile.)


texts: Extract ground elevation from a path of lat/lon

texts:
I have a set of Lat/Lon points and I'd like to know the actual ground elevation of these points.  I have access to Global Mapper, Google Earth, AutoCAD Civil 3D.  Could somebody suggest the bets way to go about doing this.  I'm trying to compare the recorded elevations from a GPS model (at the ground) to actual "ground-truth" (obviouslly based on some sort of terrain set)

Thanks


texts: Downscaling raster images in Google Earth Engine?

texts:
How can I downscale raster images using basic interpolation methods in Google Earth Engine?

For example, I want to resample Landsat image from 30m to 10m, but I cannot find corresponding functions in Google Earth Engine. Right now there are only 'resample' (cannot change the spatial resolution) and 'reduce_resolution' functions there.

Does anyone know something about it?


texts: Extract timestamp not defined using Leaflet Time Slider?

texts:
I am new to coding, so there may be many things wrong with my code.
I am trying to build a Leaflet time slider for some point data. The geojson data is of sites where water was sampled for mercury levels and represents samples taken at different locations every five years.
I want the data for each sample year to show at the same time.

So far I get points on the map, a slider that doesn't work, and an error:
SliderControl.js:183 Uncaught ReferenceError: _extractTimeStamp is not defined
how do I define the extract timestamp?

Here is a sample of the data:

     {"type": "FeatureCollection",{
"type": "FeatureCollection",
"crs": { "type": "name",
    "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },

"features": [
    { "type": "Feature",
    "properties":
        { "Site ID": "JB002",
        "Year": 1993,
        "Hg": 0.050000,
        "Place": "Jamaica Bay",
        "Start": 1993,
        "End": 1993
        },
"geometry": {
    "type": "Point",
    "coordinates": [ -73.797836, 40.63616 ] } },

and the code:

$(document).ready(function () {
var map = L.map('map').setView([40.748297,-73.941422], 11);
L.tileLayer('https://api.mapbox.com/styles/v1/eichners/cjqvwjoam3d7p2smt18p3nssn/tiles/256/{z}/{x}/{y}@2x?access_token=pk.eyJ1IjoiZWljaG5lcnMiLCJhIjoiY2lrZzVneDI4MDAyZ3VkbTZmYWlyejUzayJ9.vEGckM-D3AjV4jXmdibXyw',
{
    tileSize: 512,
    zoomOffset: -1,
    attribution: '© Mapbox © OpenStreetMap'
}).addTo(map);

var hgGeoJSON;
var sliderControl = null
var hg

// 1. ADD Hg DATA
addHg();
function addHg() {
    $.getJSON( "geojson/Hg.geojson", function( data ) {
    var hg = data;
            sliderControl = L.control.sliderControl({
                position: "topright",
                timeAttribute: "Year",
                time: '1000'
            });

    // POINT MARKER STYLE
    var hgStyle=function (feature, latlng) {
            return L.circleMarker(latlng, {
                radius: (feature.properties.Hg)*2,
                 fillColor: '#FF0000',
                 color: '#FF0000',
                 weight: 1,
                 fillOpacity: 0.6
                })
            }
    hgGeoJSON = L.geoJson(hg, {
    pointToLayer:  hgStyle,
    //onEachFeature: hgClick
    });

    //LOAD DATA LAYERS AT END TO CONTROL ORDER
    hgGeoJSON.addTo(map);
    map.fitBounds(hgGeoJSON.getBounds());

    map.addControl(sliderControl);
    console.log("slider added");

//SliderControl.js:183 Uncaught ReferenceError: _extractTimeStamp is not defined
        sliderControl.startSlider();
      })
    };
    });

and, my html with scripts, etc.

        <title>Leaflet time series _HEP Mercury</title>
        <meta charset="utf-8" />

            <link rel="stylesheet" href=jquery/jquery-ui/jquery-ui.css>
            <link rel="stylesheet" href=css/leaflet.css>
            <link rel="stylesheet" href=css/style.css>

Mercury 1993-2013


texts: How to Set Label Expression Using a Long Field Type in Toolbox

texts:
What I am hoping to be able to achieve with this Tool within my Toolbox is to allow the user to input a number of layers and to set all of the labels with the same font settings (font style, color, size, etc.). I started with the font size. Within my getParameterInfo() function, I have defined the following parameters:

        updlayers = arcpy.Parameter(
        displayName = "Layers with features to be updated",
        name = "updlayers",
        datatype = "GPValueTable",
        parameterType = "Required",
        direction = "Input")

        labelFontSize = arcpy.Parameter(
        displayName = "Label Font Size",
        name = "labelFontSize",
        datatype = "GPDouble",
        parameterType = "Required",
        direction = "Input")

Then, within my execute function, I have the following pieces of code:

 updlayers = parameters[1].values
 labelFontSize = parameters[2].value

 for layers in updlayers:
        for labels in layers:
            labels.showLabels = True
            for a in labels.labelClasses:
                    a.expression = '"%s" + [Name] + "%s"' %  ("<BOL><FNT name='Arial' size='labelFontSize'>","</FNT></BOL>")

I know that my code works fine if I replace 'labelFontSize' to an actual numerical value. I also attempted to change my expression to instead use labelFontSize as a string:

a.expression = '"%s" + [Name] + "%s"' %  ("<BOL><FNT name='Arial' size=str(labelFontSize)>","</FNT></BOL>")

However, for some reason, the label either doesn't appear at all or it shows the actual expression text as my label. The issue I am having with is obviously within the expression part of my code.

However, how do I apply my labelFontSize variable into this expression?

Am I going the correct route with this or is there an easier way to do this?


texts: I am trying to calculate the road density in my study landscape. I am using the spatial analyst, and the line density function. By default, the tool is setting a cell size of 3.6km*3.6km, and the window to 30819. The results look correct visually (high density of roads show up as such), but the values dont make sense. The range is from 0-.23, which I am interpreting as .23 km per 12.96km2. These values are very very low, and are definitely not correct. Am I making a very big oversight? The projection is in WGS_1984_World_Mercator, and the study site is in India.


texts: Fixing mercator coordinates that are outside the world bounding box?

texts:
I'm using mercator-projected overlays from a crowdsourced website where volunteers georectify images of maps by hand. Unfortunately, the slippy map that serves as the site's interface wraps and allows users to place points outside [-180, 180]/[90,-90].

I have a bunch of features with bad coordinates, like the following:

var feature = {"geometry": {"type": "Polygon", "coordinates": [[[-231.304241727798, 1682.74422719708], [825.474710540096, 1682.74422719708], [825.474710540096, -272.6456886681], [-231.304241727798, -272.6456886681], [-231.304241727798, 1682.74422719708]]]}, "type": "Feature", "properties": {"id": "7014"}}

Here's some relevant data from the WMS file:

<SRS>EPSG:4269</SRS>
<SRS>EPSG:4326</SRS>
<SRS>EPSG:900913</SRS>
<LatLonBoundingBox minx="-180" miny="-90" maxx="180" maxy="90" />
<BoundingBox SRS="EPSG:4326"
            minx="-180" miny="-90" maxx="180" maxy="90" />
...
<Layer queryable="0" opaque="0" cascaded="0">
    <Name>image</Name>
    <SRS>EPSG:4326</SRS>
    <SRS>EPSG:4269</SRS>
    <SRS>EPSG:900913</SRS>
    <LatLonBoundingBox minx="-231.304" miny="-272.646" maxx="825.475" maxy="1682.74" />
    <BoundingBox SRS="EPSG:4326"
                minx="-231.304" miny="-272.646" maxx="825.475" maxy="1682.74" />

I've tried to wrap the coordinates back into something that a map library can read, but to no avail, first using this:

lng = ((x + 180) % 360) - 180), lat = ((y + 90) % 180) - 90)

then after figuring out problems with mod in Javascript:

lng = ((((x + 180) % 360) + 360) % 360) - 180
lat = ((((y + 90) % 180) + 180) % 180) - 90

but I still get incorrect coordinates. How do I convert these out-of-bounds mercator coordinates back to something that I can put on a map?


texts: Aggregate postal code based aggregated data to district scale (MAUP)

texts: I have data aggregated to postal codes for cities in Germany. Is there anyway to aggregate this data to city districts (which look very different from postal codes spatially) in ArcGIS? It is really important for my project because I have managed to obtain some crucial data on the postal scale whereas all my research is based on the district scale.


texts:
I have installed Spectrum Server/SSA from Pitney Bowes on an AWS server and opened the ports to all traffic. The default site opens when I am logged in on the AWS VM, but not from any other computer.
I have the following in SSA:

The AWS security groups are:

I've tried using the elastic IP and the public DNS but it doesn't work. However, the default IIS install on port 80 works fine.
Any ideas?


texts:
I am using arcgis server 10.0. I published a mxd file using arcgis manager resources.
I see the web url http://cosmic-dellt410:8399/arcgis/services/DD/MapServer but it says get method not allowed. I noted that the wsdl get generated but not the rest api.  Any idea why?  How do I embed the map on my website?


texts: Loading points from an array given by PHP query from MySQL db to OL

texts:
New to OpenLayers and JS in general. I have some coords stored as lat long, then PHP script which load these points from MySQL db and simply print it as follows:

var points = [[49.66356869941062, 16.993743907660246], [49.80203898920937, 17.41534424945712], [49.7259292108726, 17.078887950628996]];

Notice that these coordinates may change or there will be much more of them and so on. So I can't put them right into my code.
I have tried to load it using this script in my HTML code but did not work.

            <script>

              var map = new ol.Map({
            target: 'map',
            layers: [
              new ol.layer.Tile({
                source: new ol.source.OSM()
              })
            ],
            view: new ol.View({
              center: ol.proj.fromLonLat([17.50, 49.25]),
              zoom: 8
            })
          });

            <?php include 'pointquery.php'; ?>
            var layer = new ol.layer.vector({
                source: new.ol.source.vector({

                })
            })

            var vectorSource = new.ol.vector({
                source: // here is probably my problem, how to define the source
                        //which in my case the array given by pointquery.php
            });

            var vectorLayer = new ol.layer.vector({
                source: vectorSource
            });
            map.addLayer(vectorLayer);
            </script>

EDIT: added the map definiton

I have found similiar soutions e.g. here or here but none of these worked for me. Also, I do not want to do this stuff via GeoJSON. I want to do this as simply as possible.

EDIT2:
Improved changed code with Reference Error: points is not defined (line points.forEach...)

<script>
    var map = new ol.Map({
                target: 'map',
                layers: [
                  new ol.layer.Tile({
                    source: new ol.source.OSM()
                  })
                ],
                view: new ol.View({
                  center: ol.proj.fromLonLat([17.50, 49.25]),
                  zoom: 8
                })
              });

                "<?php include 'pointquery.php';?>"

                var vectorSource = new ol.source.Vector({

                });
                var projection = map.getView().getProjection();

                points.forEach(function(coordinate) {
                    vectorSource.addFeature(
                        new ol.Feature(
                            new ol.geom.Point(coordinate).transform('EPSG:4326', projection)
                        )
                    );
                });

                var vectorLayer = new ol.layer.Vector({
                    source: vectorSource
                });
                map.addLayer(vectorLayer);
</script>


texts: Show attribute table in the form (table of contents)

texts:
I am building a plugin in a new window with table of contents associated. I have some actions, for example, zoom to extent and others and I want also to show the attribute table.

I search and find help in here How to show attribute table in the form, with this little code:

from PyQt4.QtGui import QApplication, QTableView
from qgis.gui import QgsAttributeTableModel
cache = QgsVectorLayerCache(self.view.currentLayer(), 10000)
model = QgsAttributeTableModel(cache)
model.loadLayer()
table = QTableView()
table.setModel(model)
table.show()

But this code doesn't show anything. Could you help me to understand what is wrong?


texts: Add field in table with layername as value

texts: I need to merge multiple tables (I'm planning to use "merge" for this). I also need a new column in the original tables that show the name of the source table. Is there a way I can add a field that equals the name of the table? What is the expression I should use in “Field calculator”? I have used hours to find out, but no luck so far =( I’m sadly not familiar with Python.


texts: Migration Issue With ArcGIS 10 in ExportOperationClass

texts:
I am migrating my Application built in Version 9.3.1 to version 10
In my code I am getting execption:

{"Retrieving the COM class factory for component with CLSID {4147A362-EFC4-11D3-
 A0A4-00C04F6BC626} failed due to the following error: 80040111."}

at this place

IExportOperation pExportOp = new ExportOperationClass();

I have alreday bind my application with bindmanager code

ESRI.ArcGIS.RuntimeManager.Bind(ESRI.ArcGIS.ProductCode.Engine);

Every thing works fine like creating IWorkspaceFactory InitiallizeLicense, exception arise as this statment is executed.

IExportOperation pExportOp = new ExportOperationClass();

I found at google that beside binding , some times we have to use

ESRI.ArcGIS.VersionManager.LoadVersion(ESRI.ArcGIS.ProductCode.Engine, "10.0");

but I am not able to access ESRI.ArcGIS.VersionManager

Is there any solution to avoid exception?


texts: Why does clipping netCDF in QGIS change cell values?

texts:
My problem goes like this:

I have a gridded data in netCDF file. The cells have values in meters. For example: 0,012. I only need a small area of the grid and after clipping it in QGIS using a shapefile mask the values are multiplied with 10 000, result for the same cell being 120. I have to mention that the clipping is done on all the bands (365/366) without creating an alpha band. GDAL code looking fine.

Any ideas why it goes like this?


texts: Can't add marker saved as an EMF to the style manager

texts: I'm trying to add a marker I created in PowerPoint and saved as an EMF to ArcMap. The current questions all state that all I have to do is open the style manager and go to the marker symbols folder, and add my new symbol there. However it won't give me the option to add a new symbol as this option is un-clickable. I'm not sure why. Any thoughts?


texts: Shifted shapefiles with same coordinate system in ArcGIS

texts:
I have a strange situation with one of my GIS projects. I received from a client three shapefiles in the same coordinate system (Stereographic 1970 - a projected system). All three files should align perfectly. This happens only for two of them, the third one is shifted from the rest of the project.
I`ve tried to apply different methods of converting the coordinate system, reproject but no success till now.

Apparently, the shapefile that is not overlapping with the other two, is coming from a reconversion of a different coordinate system.
Do you have any solution that might work in ArcGIS? Or Global Mapper?


texts: Creating data driven field on attribute table in QGIS

texts: I want to create a Field on an attribute table, whose data will depend on another field within the same attribute table.


texts:
I am using Subgurim.Controles google maps in my project,i was updating every 40 with new marker to the Gmap.

I just want to update position of marker every 40 seconds, not refresh entire Gmap. Here i am using timer for every 40 sec to get new latlng from database table;

aspx code:

<asp:Timer ID="Timer1" runat="server" OnTick="Timer1_Tick">
             </asp:Timer>
             <div id="map">
             cc1:GMap ID="GMap1" runat="server" Height="100%" Width="100%" />
             </div>

C# code:

  protected void Timer1_Tick(object sender, EventArgs e)
{  GMap1.resetMarkers(); Subgurim.Controles.GLatLng gLatLng = new
    > Subgurim.Controles.GLatLng(Convert.ToDouble(status.Rows[0][0]),
    > Convert.ToDouble(status.Rows[0][1])); GMap1.setCenter(gLatLng, 16,
    > Subgurim.Controles.GMapType.GTypes.Normal); Subgurim.Controles.GMarker
    > oMarker = new Subgurim.Controles.GMarker(gLatLng);
    > GMap1.addGMarker(oMarker); }


texts: Google Earth Extension to export georeferenced image?

texts: Has somebody coded an extension to Google Earth that will export a georeferenced image from google earth? Google Earth and Google Earth Pro only export the image *.jpg *.png etc and don't export the jgw or pgw file that contains the georeferencing information. This would be a very handy extension to have since google earth pro usually has the most up to date, freely available Imagery for the regions that I work in.


texts:
After a gruesome night installing PostgreSQL 9.1, PostGIS 2 and Nominatim 2 on CentOS 6.3 using the instructions here, I am finally at the stage of importing a US planet file.

Problem: In ~/Nominatim/utils, When I attempt to run ./setup.php --osm-file us-northeast.osm.pbf --all, the output I get is:

["type"]=>
      string(2) "->"
      ["args"]=>
      array(2) {
        [0]=>
        array(9) {
          ["phptype"]=>
          string(5) "pgsql"
          ["dbsyntax"]=>
          string(5) "pgsql"
          ["username"]=>
          string(0) ""
          ["password"]=>
          bool(false)
          ["protocol"]=>
          string(3) "tcp"
          ["hostspec"]=>
          string(0) ""
          ["port"]=>
          bool(false)
          ["socket"]=>
          bool(false)
          ["database"]=>
          string(9) "nominatim"
        }
        [1]=>
        bool(false)
      }
    }
    [5]=>
    array(6) {
      ["file"]=>
      string(35) "/home/myusername/Nominatim/lib/db.php"
      ["line"]=>
      int(7)
      ["function"]=>
      string(7) "connect"
      ["class"]=>
      string(2) "DB"
      ["type"]=>
      string(2) "::"
      ["args"]=>
      array(2) {
        [0]=>
        string(19) "pgsql://@/nominatim"
        [1]=>
        bool(false)
      }
    }
    [6]=>
    array(4) {
      ["file"]=>
      string(40) "/home/myusername/Nominatim/utils/setup.php"
      ["line"]=>
      int(118)
      ["function"]=>
      string(5) "getDB"
      ["args"]=>
      array(0) {
      }
    }
  }
  ["callback"]=>
  NULL
}
DB Error: connect failed

and I am thrown back to the command line. I am logged in as a user who has a postgresql role of the same name and is a superuser. Any idea what went wrong? Do I need to create any database before doing this?


texts: Finding cell width of Raster using ArcPy with ArcGIS Pro?

texts:
I've been trying to find the cell size, or height/width/mean, of a cell by utilizing arcpy in the Python Window in ArcGIS Pro version 2.3.3., I run the following code:

import arcpy
env = arcpy.env
env.workspace = "D:\Fall 2019\GEOG 625\pr\pen"
raster = "D:\Fall 2019\GEOG 625\pr\pen\LandsatHMB201707"
description = arcpy.Describe(raster)
arcpy.env.Describe("D:\Fall 2019\GEOG 625\pr\pen\LandsatHMB201707").meanCellHeight

I've looked into the raster properties here
https://pro.arcgis.com/en/pro-app/arcpy/functions/raster-band-properties.htm

and the error code I get is:

  AttributeError: DescribeData: Method meanCellWidth does not exist


texts:
I am very new to GIS so please forgive me if the question is naive, but I haven't found a solution yet.

I generated a mock hexbin in QGIS using MMQGIS plugin, with project set to EPSG:4326 and saved to PostGIS.

When queried as Geojson via Postgis, the returned data looks stretched in Google Maps.

I have tried to save in EPSG:3857 Google Mercator but that doesn't work either.

Creating the map in 3857 was the answer, however I did try before without success, however trying via DB Manager did the trick, the Target SRID was the missing piece in the puzzle!


texts: How to proxy tile sources with reversed Y axis TMS source

texts:
I have the following URL https://htonl.dev.openstreetmap.org/ngi-tiles/tiles/50k/{z}/{x}/{-y}.png which I am trying to cache using mapproxy. I have the following config

sources:
ngi_tiles:
    type: tile
    grid: osm_grid
    url: https://htonl.dev.openstreetmap.org/ngi-tiles/tiles/50k/%(z)s/%(x)s/%(-y)s.png

grids:
   webmercator:
      base: GLOBAL_WEBMERCATOR
   osm_grid:
     base: GLOBAL_MERCATOR
     srs: 'EPSG:3857'
     tile_size: [512, 512]
     origin: nw

The format of the ULR seems to generate an error with mapproxy

File "/usr/local/lib/python3.7/site-packages/mapproxy/client/tile.py", line 81, in substitute
    return self.template % data
KeyError: '-y'

I have tried without the inverted y-axis but changing the origin of the grid from nw to other supported formats but it still does not work.


texts: Setting up PyQGIS 3 (QGIS 3.12) interpreter in PyCharm 2018.3.7

texts:
I cannot set up the correct interpreter for PyQGIS 3 in PyCharm 2018.3.7.

Using: QGIS 3.12.0, OSGeo4W64, Python 3.7 (part of the OSGeo4W64 folder), PyCharm 2018.3.7.

My .cmd files in OSGeo4W64 look like this:

pyqgis.cmd:

SET OSGEO4W_ROOT=C:\OSGeo4W64
call "%OSGEO4W_ROOT%"\bin\o4w_env.bat

@echo off
path %PATH%;%OSGEO4W_ROOT%\apps\qgis\bin
path %PATH%;%OSGEO4W_ROOT%\apps\grass\grass-78\lib
path %PATH%;C:\OSGeo4W64\apps\Qt5\bin
path %PATH%;C:\OSGeo4W64\apps\Python37\Scripts

set PYTHONPATH=%PYTHONPATH%;%OSGEO4W_ROOT%\apps\qgis\python
set PYTHONHOME=%OSGEO4W_ROOT%\apps\Python37

set PATH=C:\Program Files\Git\bin;%PATH%

cmd.exe

pycharm.cmd:

SET OSGEO4W_ROOT=C:\OSGeo4W64
call "%OSGEO4W_ROOT%"\bin\o4w_env.bat

@echo off
path %PATH%;%OSGEO4W_ROOT%\apps\qgis\bin
path %PATH%;%OSGEO4W_ROOT%\apps\grass\grass-78\lib
path %PATH%;C:\OSGeo4W64\apps\Qt5\bin
path %PATH%;C:\OSGeo4W64\apps\Python37\Scripts

set PYTHONPATH=%PYTHONPATH%;%OSGEO4W_ROOT%\apps\qgis\python
set PYTHONHOME=%OSGEO4W_ROOT%\apps\Python37

set PATH=C:\Program Files\Git\bin;%PATH%

start "PyCharm aware of QGIS" /B "C:\Program Files\JetBrains\PyCharm Community Edition 2018.3.7\bin\pycharm64.exe" %*

After running pyqgis.cmd I can import different QGIS packages so everything seems ok.
Running the pycharm.cmd file I can directly access PyCharm.

In PyCharm I´ve just defined the location of the new project and the path for the project interpreter which should be located in the OSGeo4W64 folder:

C:\OSGeo4W64\bin\python3.exe

But in the system settings the project interpreter ist still different:

The project interpreter is not the base interpreter. After trying to add a new interpreter I onl get the pip and and setuptool python packages and to not the PyQGIS packages:

After creating a new .py file in the project I can run some PyQGIS code (import some QGIS packages e.g.) but cannot run the code in the script window.

Maybe the SDK is missing but I didn´t get any error messages so I think this is not the issue.

This is really frustrating. Could you please tell me what I´m doing wrong? I´d really appreciate it.
The .cmd files are based on:

Environment Setup for Plugin Development with QGIS 3

Some GIS StackExchange references (unfortunately haven´t bring me further):

Set up PyCharm for PyQGIS 3 development

Adding QGIS-specific Python interpreter to PyCharm on Windows?


texts: How to duplicate feature-linked annotation linked to the same feature class in 9.31?

texts: I am trying to make many feature-linked annotation feature classes linked to the same (point) Feature Class, the symbology is exactly the same, but I need separate anno's for users so they can create their own queries and placements.  I don't want to have to re-create all the symbology each time.  The export function is grayed out in AM.  And copy/paste in AC creates a new FC as well (I want them all linked to the same FC), I can't change the relationships either in AC.  Help??  Thanks.


texts: What other tools can complete this exercise besides ArcGIS Spatial Analyst?

texts:
excel spread sheet here

The objective in the exercise is to determine how many public computers are available per person as a result of spatial distribution of libraries and people.

Part 1: Grid of points

Using Brooklyn, NY boundary, create a grid of points spaced at 250’.  Please report the number of grid points resulting.

Part 2:  “Smooth” population on a uniform grid of points

Using Brooklyn, NY boundary, census block population data (US Census 2000) and the resulting grid of points from Part 1, construct a “smoothed” population density for Brooklyn.  Using Spatial Analyst create a population raster from census blocks.  Note the default density and make sure it corresponds to the 250’ grid point.  Apply the population raster to the vector grid making sure that the resulting total population is equivalent to original data.  Please show results and document steps in detail.

Part 3 :  Determining Library computers per population

Using Brooklyn , NY boundary, uniform population grid resulting from Part 2 and the attached table, determine the interpolated density of publicly available computers from Brooklyn’s Public Library System per population (using the uniform grid of points population).  Please show results and document steps in detail.


texts: Text Label GeoJSON Markers Google Maps API

texts:
I am trying to create a Text Label, for each of the markers created from my GeoJSON with the following code...

function init() {
  var mapCanvas = document.getElementById('myMap');
  var map;
  var image = "http://maps.google.com/mapfiles/ms/micons/blue.png";

  var userCenter = new google.maps.LatLng(-29.789615, 30.822422);
  var mapOptions = {
    draggable: true,
    zoomControl: true,
    scrollwheel: true,
    disableDoubleClickZoom: false,
    zoom: 18,
    center: userCenter
  };

  map = new google.maps.Map(mapCanvas, mapOptions);

  var customLayer = new google.maps.Data();

  map.data.addGeoJson(jsonData);
  map.data.setStyle({
    title: '#',
    icon: image,
    map: map,
  });

  map.data.forEach(function(feature) {
    var point = new google.maps.LatLng(feature.getProperty('coordinates'))
    var name = feature.getProperty('GPSLogs_ro')
    var mark = new google.maps.Marker({
      position: point,
      title: '#',
      icon: image,
      map: map,
      draggable: false,
    });
  });
  // customLayer.setMap(map);
}
google.maps.event.addDomListener(window, 'load', init);
var jsonData = { //only showing one point from json to save space
  "type": "FeatureCollection",
  "crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
  "features": [
   { "type": "Feature",
     "properties": { "gpslogid": 106,
                     "createddat":   "2015\/07\/02 13:37:40.000",
                     "lat": -29.788272, "lon": 30.823182,
                     "GPSLogs_st": "2015\/07\/02 13:36:49.000",
                     "GPSLogs_ro": "IMPANGELE",
                     "GPSLogs_is": "False\r",
                     "id": 622,
                     "xmin": 30.821820, "xmax": 30.821920,
                     "ymin": -29.789654, "ymax": -29.789554 },
     "geometry": { "type": "Point",
                   "coordinates": [   30.821919004535157, -29.789652962174948 ] } }
  ]
}

This code is adapted from https://stackoverflow.com/questions/32893734/google-maps-javascript-v3-21-trouble-displaying-geojson-features,

Now that I have my markers adding to the map, I want to assign them a Text Box type label, containing the var name = feature.getProperty('GPSLogs_ro') content.

The closest thing I have found is https://stackoverflow.com/questions/3953922/is-it-possible-to-write-custom-text-on-google-maps-api-v3 however I try to apply it, and no labels appear, just the markers, I am sure I am doing it wrong!


texts: PostGIS: spatial analysis of point layer SRID 4326 using UTM grid

texts:
Edited, as suggested.

I've got a PostgreSQL/PostGIS table 'Datos'

ALTER TABLE datos ADD COLUMN geom GEOMETRY;
UPDATE datos SET geom = ST_MakePoint(decimallongitude,decimallatitude);
update datos set geom = ST_SETSRID(geom,4326);
CREATE INDEX idx_geom ON datos USING gist(geom);

I need to do some analysis by means of a UTM grid 10x10 km that I imported from a shapefile using QGIS. The table details are

The table SRID is 25830, so I projected the table 'Datos' to SRID 25830 using

SELECT UpdateGeometrySRID('datos','geom',25830);

Now, I need to count points within the grid cells, but I only get "0" values in the count as if they don't overlap (but they do!). I've tried several combinations of queries based on

  SELECT utmcode, count(datos.geom) as count
  FROM grid
  LEFT JOIN datos
      ON st_contains(grid.geom,datos.geom)
  GROUP BY utmcode;

What is wrong in all this? Projection issues, query syntax or both? How can it be fixed?


texts: Geoserver - cache tiles with background color

texts:
I have degines style for my osm data, but background is transparent for png tiles. Is possible to seed tiles with background color?

I want seed tiles with same result like calling wms with these parameters:

  ...&&BGCOLOR=0xfefefe&TRANSPARENT=false

current style:

black places are in tiles transparent, and my client app can not change this color, so I generate tiles with native background.


texts: Orfeo toolbox (OTB) algorithms in QGIS 3.0?

texts:
I've been trying to understand if the OTB (Orfeo toolbox) algorithms can be used in QGIS 3.0.0-Girona but found no official documentation about it. Under the 'Available QGIS Data Provider Plugins' (Help--> About), OTB is not showing.

It also doesn't help to install via OS4GEOW as mentioned in this post, as this only works for QGIS 2.X (which did work in including the OTB algorithms for me). Furthermore, OTB nowadays shows as deprecated.

Does anyone know how to include the algorithms and/or whether a future integration is planned in a future QGIS 3.X version?


texts: ArcGIS Pro Charts Auto Rounding Inconsistently

texts:
I'm making a set of maps and charts in ArcGIS Pro.  I have a scatterplot with a r2 value of 0.06 and another with a value of 0.08.  When I make chart panes, it displays the r2=0.06 as is on the chart, but it auto rounds the r2=0.08 to 0.1 on the chart within the layout.
Is there any work around to fix this glitch/inconsistent behaviour?
I'm doing this as part of a Lab assignment for my students.  I recognize the r2 values aren't really significant, I just want to show them how to do it.


texts: QGIS: Cropland Data Layer - how to add class names to .tif file?

texts:
I'm working with a clip of the Cropland Data Layer (raster), which has a number for each value coded to a crop type. I want to add these crop type / class names to the .tif file (and have them show up in the "Identify Results" pane). the CDL has a nice explanation for doing this in ArcGIS:

  f your downloaded CDL .tif file does not contain category names, then you can add them using the following instructions. Download the file: generic_cdl_attributes.tif.vat.dbf. This generic file contains all possible CDL colors and category names. As long as the .tif file and the .tif.vat.dbf file have the same file name, then the category names will load automatically in ArcMap. So, change the file name (not extension) of the generic_cdl_attributes.tif.vat.dbf to match the file name of the downloaded CDL .tif file. Then add the .tif file as a layer in ArcMap. The category names will display in the Table of Contents window.

I have tried this method in QGIS but can't get it to work. Here's a screenshot showing that only the value, not the class, is displayed in the "identify results" pane:

Has anyone encountered this before?


texts: Get exterior ring of a multipolygon using GEOSGeometry or PostGIS

texts:
I am trying to get the exterior ring of the multipolygon in the image below using GEOSGeometry (through GeoDjango) or Postgis, but I have not found an easy solution.

This is an example of a dynamically generated polygon in my django code, so the shapes can be different in other cases. The final goal is to get the length of the red line (sorry about the poor GIMP job!).

I did not find a solution using Django GEOGeometry.

As for Postgis, the best I could do was removing the holes, making a polygon from my shape, buffering a bit, and using the boundary of the final shape to calculate the length. Here's the query for this process:

SELECT st_asgeojson(
           ST_Boundary(
               ST_Buffer(
                   ST_Collect(
                       ST_MakePolygon(geom)
                     ),
                   0.0001
                 )
             )::geography
         ) As geom
FROM (
       SELECT ST_ExteriorRing(
                  (ST_Dump(
                      ST_SetSRID(
                          ST_GeomFromGeojson(
                              '{ "type": "MultiPolygon", "coordinates": [...]}'
                            ),
                          4326
                        )
                    )).geom
                ) AS geom
     ) AS foo

And here is the geojson object for the example in the screenshot:
https://gist.github.com/ka7eh/427d071c7077340cdc128f3624b30777

So my questions are:

Is there any way to do this directly in Python?
If not, is there a more straight-forward way of handling this in Postgis?
Are there any downfalls to my solution? I mean whether there might be situations that my query returns unexpected results? The only tricky part I can think of right now is the buffer value and how it might mess up my geometry.


texts: I'm new to Python but I've heard is a good way to automate tasks. I have thousands of point shapefiles and need to extract the maximum value of each one and put them in a new table. I'll appreciate your help.


texts:
I've been designing some tilesets in Mapbox Studio Classic and I keep seeing the :: notation in styles that Mapbox has created.

What does that mean? Are there only specific values that can be put after it?


texts: Help installing add-on in GRASS, I get an error message

texts: My Grass gives me this error message, Grass_info_error(4816,1):'make' required. Please install  'make' first. This happens when I attempt to fetch an add-on (r.viewshed). What am I doing wrong? and what is 'make'? Thanks for helping for a Grass newbie. I'm running on Mac Snow Leopard if that makes any difference.