texts:
So I am trying to display the appropriate legend for my first map. This is my javascript for my first_map.

function firstmap() {
    var first_map = new ol.Map({
        target:'FirstMap',
        layers: [
            new ol.layer.Tile({
                title: 'OSM',
                source: new ol.source.OSM()
            }),
            new ol.layer.Group({
                title: 'Overlays',
                layers:[new ol.layer.Tile({
                    title: '- NDVI',
                    source: new ol.source.TileWMS({
                        url: 'http://mygeonodeserver.org/geoserver/wms',
                        params: {layers: 'geonode:xyz',
                                      format: 'image/png'},
                        serverType: 'geoserver' }),
                    showLegend:true
                })]
            })
        ],
        view: new ol.View({
            projection: 'EPSG:900913',
            center:[3930427, 359171],
            zoom:6
            })
    });
    first_map.addControl(new ol3_legend({
    map: first_map,
    class: 'ol_legend' }));
}

The get legend java script i got online is this one below

<script type = 'text/javascript'>
function ol3_legend(e) {
    var options = e || {};
    var wmsVersion = options.wmsVersion || '1.1.1';
    var format = options.format || 'image/png';
    var legendP = document.createElement('p');
    legendP.innerHTML = 'Legend:';
    var legendDiv = document.createElement('div');
    legendDiv.className = options.class + ' ol-unselectable';
    legendDiv.appendChild(legendP);
    var layers = options.map.getLayers().getArray();
    for (var i = 0; i < layers.length; i++) {
        if (layers[i].get('showLegend') === true) {
            try {
                var url = layers[i].getSource().getUrls()[0];
            } catch (err) {
                var url = layers[i].getSource().getUrl();
            }
            var legendImg = document.createElement('img');
            legendImg.src = url + '?REQUEST=GetLegendGraphic&sld_version=1.0.0&layer=' + layers[i].getSource().getParams().layers + '&format=' + format;
            legendDiv.appendChild(legendImg);
        }
    }
    ol.control.Control.call(this, {
        element: legendDiv
    });
}
ol.inherits(ol3_legend, ol.control.Control);

proj4.defs("EPSG:900913", "+proj=somerc +lat_0=47.14439372222222 +lon_0=19.04857177777778 +k_0=0.99993 +x_0=650000 +y_0=200000 +ellps=GRS67 +towgs84=52.17,-71.82,-14.9,0,0,0,0 +units=m +no_defs");
</script>

The above script is not working because of my layers are grouped i.e new ol.layer Group


texts: Spatial Query plugin broken in MacOS version?

texts:
I'm running QGIS 1.7.4-4 on a MacOS 1.6.8.
I'm trying to use the "Spatial Query" plugin but it's not loading.
The plugin is activated and I can see it in the main menu Plugins > Spatial Query > Spatial Query but when I click on it no window appears, nothing happens.
I've tried to deactivate it, restart QGIS, reactivate, but nothing changes.

Any help?


texts: Calculating percentage of polygons within specific areas using ArcGIS for Desktop?

texts:
I have a parcel layer and I am trying to determine the percentage of each parcel that falls within specific zones.  I tried tabular intersect but it just returned an empty output.

I'm using ArcGIS 10.2 for Desktop.


texts:
I would like to export a shapefile out of the Graphical Modeler.
When I export the data as a file, I only can give a fixed name. Is there a possibility that I can give the export a variable name? I don't want use the Python script.


texts: GeoTools: SimpleFeatureCollection to KML error CoordinateReferenceSystem

texts:
I want to export a SimpleFeatureCollection (Read from a MySQL DB) to a KML file, here is my code used to read the infos into the console:

private void featureCollectionToKML(SimpleFeatureCollection features)
        throws Exception {

    Encoder encoder = new Encoder(new KMLConfiguration());
    encoder.setIndenting(true);
    encoder.encode(features, KML.kml, System.out);
}

When I run the program there are first some WARNINGS when I get the Features and then the program crashes with this error log:

  PM org.geotools.referencing.factory.PropertyAuthorityFactory

  WARNING: Axis elements found in a wkt definition, the force longitude
  first axis order hint might not be respected: PROJCS["Google
  Projection",GEOGCS["WGS 84",DATUM["World Geodetic System
  1984",SPHEROID["WGS 84", 6378137.0, 298.257223563,
  AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",
  0.0, AUTHORITY["EPSG","8901"]],UNIT["degree", 0.017453292519943295],AXIS["Geodetic longitude", EAST],AXIS["Geodetic latitude", NORTH],AUTHORITY["EPSG","4326"]],PROJECTION["Popular
  Visualisation Pseudo Mercator",
  AUTHORITY["EPSG","1024"]],PARAMETER["semi_minor",
  6378137.0],PARAMETER["latitude_of_origin", 0.0],PARAMETER["central_meridian", 0.0],PARAMETER["scale_factor", 1.0],PARAMETER["false_easting", 0.0],PARAMETER["false_northing", 0.0],UNIT["m", 1.0]]

  ott 23, 2013 4:06:49 PM
  org.geotools.feature.simple.SimpleFeatureTypeBuilder add

  WARNING: Creating LookAt with null CoordinateReferenceSystem - did you
  mean to setCRS?

  ott 23, 2013 4:06:49 PM
  org.geotools.feature.simple.SimpleFeatureTypeBuilder add

  WARNING: Creating Region with null CoordinateReferenceSystem - did you
  mean to setCRS? java.lang.RuntimeException: Failed to get property: id
    at
  org.geotools.xml.impl.GetPropertyExecutor.visit(GetPropertyExecutor.java:105)
    at
  org.geotools.xml.impl.BindingWalker$BindingExecutionChain.execute(BindingWalker.java:220)
    at org.geotools.xml.impl.BindingWalker.walk(BindingWalker.java:186)
    at
  org.geotools.xml.impl.BindingVisitorDispatch.walk(BindingVisitorDispatch.java:47)
    at
  org.geotools.xml.impl.BindingVisitorDispatch.walk(BindingVisitorDispatch.java:32)
    at org.geotools.xml.Encoder.encode(Encoder.java:902)    at
  org.geotools.xml.Encoder.encode(Encoder.java:603)     at
  shapefileManager.ShapefileManager.featureCollectionToKML(ShapefileManager.java:173)
    at
  shapefileManager.ShapefileManager.convertArchiToFeatuesKML(ShapefileManager.java:114)
    at Main.Test.main(Test.java:19) Caused by:
  java.lang.UnsupportedOperationException   at
  org.geotools.data.store.ContentFeatureCollection.getID(ContentFeatureCollection.java:492)
    at
  org.geotools.kml.bindings.FeatureTypeBinding.getProperty(FeatureTypeBinding.java:284)
    at
  org.geotools.xml.impl.GetPropertyExecutor.visit(GetPropertyExecutor.java:103)
    ... 9 more

I can't understand which is the problem...

Here is how I get the SimpleFeaturesCollection from MySQL DB

 public void convertArcsToFeatuesKML(String[] ids) {

    try {

        FilterFactory2 ff = CommonFactoryFinder.getFilterFactory2();

        Set<FeatureId> fids = new HashSet<FeatureId>();

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

            fids.add(ff.featureId(ids[i]));

        }

        Filter filter = ff.id(fids);

        SimpleFeatureCollection features = featureSource
                .getFeatures(filter);

        this.featureCollectionToKML(features);

    } catch (Exception e) {             e.printStackTrace();        }
}

FeatureSource is the MySQL DB.

PS I've checked, with the getId() method, and all features have an unique ID.


texts: How to hide/suppress and rename/alias fields when exporting to KML?

texts: I have a feature class that needs to be exported to a KML file. I wonder if there is a way for me to export only some of the columns of the feature class's attribute table; some attributes are not needed in the KML file. Also, some of the column headings of the attribute table need to be edited to make them more understandable by a layperson. Is there a way of doing this when exporting to KML? Or, is there a way to make these edits in a KML file once it already exists? I'm exporting from ArcGIS 10.0, in case this helps.


texts: How can I show a text column on MapServer?

texts:
I'm working on a project with MapServer with PostGIS and I need to retrieve from my table and display a layer that shows some information like a street name.

What I have? I have my map rendered perfectly, showing all streets of a city, and I need to define a layer to show on the rendered street their names.

I've searched in MapServer Documentation about text layers but I did not found anything that actually help to solve my problem.

I've read the link suggested in the comments of this question and make this code, but it does not work.

  LAYER
    NAME face_guarulhos
    TYPE line
    CONNECTIONTYPE postgis
    CONNECTION "db stuff config"
    DATA "geomlatlon from ctageo.face_guarulhos"
    STATUS default

        METADATA
            "wms_title" "face_guarulhos"
        END

    CLASS
      NAME "Faces Guarulhos"
            OUTLINECOLOR 139 134 130
      COLOR 255 145 145

            LABEL
                COLOR 255 255 255
                OUTLINECOLOR 0 0 0
                SIZE 10
                ANGLE 0
                POSITION AUTO
                PARTIALS false
            END
  END

END


texts: Renaming attribute via pyQGis and data provider

texts:
I'm looking for a solution to rename a field/attribute via dataprovider (see the code below):

def modifierChamp(self,file_table,nom_champ,nom_champ_new):
    #vÃ©rifie que le champ existe
    dp = file_table.dataProvider()
    ind = dp.fieldNameIndex(nom_champ)
    if ind==-1:
        self.logFile.write("Le champ "+nom_champ+u" Nexiste PAS dans la table "+file_table.name()+" \n")
    else:
        # renameAttributes exists but I don't find how to use it !
        dp.renameAttributes(QgsFieldNameMap(ind,nom_champ_new)) # it does not work
        dp.updateFields()
        self.logFile.write("Correction du nom du champ "+nom_champ+" par "+nom_champ_new+" dans la table "+str(file_table.name())+" \n")


texts: Selecting Parcels by Watershed boundary - QGIS 3.4

texts: I am using QGIS 3.4. I would like to select all the property parcels withing a watershed boundary, easily done with Clip. But I want to retain the complete parcels that cross the watershed boundary and not have them clipped to the boundary.


texts: Are any fictional worlds available in some GIS format?

texts:
Is there any fictional world (e.g. Middle-earth, Discworld etc.) which could be used in a GIS context?

For example, an MBTiles archive and a shapefile or any other format usable with OpenLayers / Polymaps / Leaflet?


texts: Add points at specific gps coordinates in qgis

texts:
I created a layer in order to draw points with the azimut and distance plugin. IT didnt work at all. There was not even 1 point created. I followed the tutorial here: http://www.teuderun.de/workshops/qgis/aufschlusspunkte-hinzuf%C3%BCgen/

I entered the gps data in decimal format as suggested in the tutorial.
There is no draw button. No points are drawn. The layer is in active mode.

Are there any other options to draw a few points using their coordinates?


texts: Getting points that overlap with polygon using PostGIS?

texts:
I have two geometries:

multi point A

multi polygon B

I want to get all the points of A that overlap with the polygon B

I am not sure how to run this query. I have tried the following query:

SELECT ST_Intersection(A.geom, B.geom)
FROM point AS A, poly AS B

But with little succes. How can I change it so that i get only the points of A that are overlap with B?


texts: Calculating correct longitude when it's over |180|?

texts:
I'm trying to develop a "formula" to correct the lat-lng values.

I'm using vue-leaflet but when you pan outside the "first" world you get big numbers. Over +180 or under -180.

For example: when I pan to America to the right (east direction), I get as lng 215.
In my mind, I would just correct it with 215-360=-145

The same is for when I pan to east russia to the left (west direction) and I get for example -222. Now I need to calculate -222+360=138

However, since the world is indefinite the user could pan to the 8th world and I had to adjust the values.

Is it possible to calculate the right longitude? (and another requirement is when the user is in the first world, 24 lng should still be 24 lng.


texts: Where i will get English setup for gvSIG?

texts: I am trying to install gvSIG. It is asking me which setup language, i selected English. But when installation is complete it is giving me interface with Russian language. Can any one help me to install or download the English version setup?


texts: How can I export ArcScene animation to .mov quicktime file?

texts:
ArcGIS 10 documentation says that an animation can be exported to either .avi or .mov (quicktime).  I only see the .avi option in my drop-down list?  What am I missing.

EDIT: Quicktime is installed on my computer


texts: Add closest facility as an attribute to a point in ArcMap?

texts:
Is there a way, after performing a closest facility analysis, to add/save the closest facility number/name/unique ID to the incident?

I plan to then use this attribute to create a route going from one incident to another incident via this closest facility, which means I have to somehow save it / keep track of which facility is closest for all incidents.


texts: Is there a way to enable CSV or JSON output formats for MapServer WFS GetPropertyValue requests?

texts:
In a MapServer map file, it is possible to set up a variety of output formats for WFS GetFeature requests.

However, it doesn't seem possible to do the same for a WFS GetPropertyValue request, which only seems to work with these output formats:

application/gml+xml; version=3.2
text/xml; subtype=gml/3.2.1
text/xml; subtype=gml/3.1.1
text/xml; subtype=gml/2.1.2

Is there any way to enable CSV or JSON output for GetPropertyValue requests?


texts: How to do spatial adjustment of Line features by edge snapping in ArcGIS?

texts:
I have two polyline features of same area which are prepared based on two different raster imagery. Hence the two features are not exactly matched and has some shift.

How can I adjust them to overlay exactly one over the other? I attached the snapshot of the two polyline features.

I want Rd2 in the below imagery to adjust with Rd1 while Rd1 remains unchanged.


texts: Creating snow cover frequency map in Google Earth Engine

texts:
I would like to create a probability map from NDSI_Snow_Cover for a school project. I am VERY new to object-oriented programing and remote sensing data.

I used the code below, but I don't know how to deal with each image in the returned search individually (if I need to). My thoughts are to create an image of snow / no-snow for each image returned in the search. Then for a given pixel in all the snow/no-snow images, determine the number of times there was snow divided by the number of samples (excluding results where the pixel is masked).

Does that make any sense?

var MOD10A1_Terra = ee.ImageCollection("MODIS/006/MOD10A1");

// Setting date range, the month number, and the day number
var start = ee.Date('2000-1-01');
var end = ee.Date('2017-1-02');
var month = 1;           // The number of the starting month; Jan = 1
var numberOfMonths = 1;  // The number of the last month in search
var dayOfMonth = 1;
var numberOfDays =1;

// Filtering and Sorting MOD10A1_Terra
// the variable filteredTerra can contain multiple images
var filteredTerra = ee.ImageCollection(MOD10A1_Terra)
//.filterBounds(region)
.filterDate(start, end)
.filter(ee.Filter.calendarRange(month,numberOfMonths,'month'))
.filter(ee.Filter.calendarRange(dayOfMonth,numberOfDays,'day_of_month'))
.select('NDSI_Snow_Cover');

print('NDSI_Snow_Cover images returned from search parameters:', filteredTerra);

The code above as written will filter the image collection to 17 images. What I think I want to next is to assign each image returned by the code above to a variable.

For example:

var image1 = the first image from the image collection filteredTerra
var image2 = the second image from the image collection filteredTerra
var image3 = the second image from the image collection filteredTerra…and so on.

How would this be done in google-earth-engine?


texts: How to create a Spatial Reference object from Calculate UTM Zone output?

texts:
I have a function that extracts the well known text representation of a UTM zone:

def find_utm_zone(shp):
    """
    -Generates a UTM zone for the feature
    """
    arcpy.env.overwriteOutput = 1
    gdb = arcpy.env.scratchGDB
    fc = os.path.join(gdb, "temp_fc")
    arcpy.CopyFeatures_management(shp, fc)
    arcpy.AddField_management(fc, "utmzone", field_type = "TEXT", field_length = 1000)

    # Execute CalculateUTMZone
    arcpy.CalculateUTMZone_cartography(fc, "utmzone")

    prj = [row[0] for row in arcpy.da.SearchCursor(fc, "utmzone")][0]
    return prj

The output looks like the following (i.e. the same as a .prj file):

PROJCS["WGS 84 UTM Zone 18L (Calculated)",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",10000000.0],PARAMETER["Central_Meridian",-75.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]

When I try to create a Spatial Reference object from the output I get a generic RuntimeError: ERROR 999999: Error executing function. error:

shp = r'C:\path\to\myshapefile.shp'
sr = find_utm_zone(shp)
print arcpy.SpatialReference(sr)

What is the preferred method of creating a Spatial Reference object from the Calculate UTM Zone output?


texts: Intersection of touching geometries returns a polygon

texts:
I'm using Postgis 2.3.
I have two touching geometries that I got from an overlay (ST_Difference) operation:

id |                                              st_astext
 3 | MULTIPOLYGON(((0 -2,0 4,3 4,6 1,0 -2)),((6 6,8 8,8 4,6 6)))
 4 | POLYGON((2 6,7 9,7 7,6 6,7 5,7 0,4 0,6 1,4.33333333333333 2.66666666666667,5 4,2 6))

Yet when I check their intersection, I get not only linestrings, but also a polygon.

select st_AsText(ST_Intersection(a.the_geom, b.the_geom)) from renovation a, renovation b where a.id=4 and b.id=3;

GEOMETRYCOLLECTION(LINESTRING(7 7,6 6),LINESTRING(6 6,7 5),LINESTRING(4 0,6 1),POLYGON((6 1,4.33333333333333 2.66666666666667,4.33333333333333 2.66666666666667,6 1)))

This polygon is very weird; it's got two equal vertices, yet it is valid, and it has a non-zero area.

How does this happen, and how can I prevent that from happening or filter this unnecessary polygon out?


texts: How to iterate over features and abort iteration with keypress in PyQGIS?

texts:
I would like to run a loop with PyQGIS, it should wait for the user to press a key either SpaceBar or Escape to continue.

canvas = qgis.utils.iface.mapCanvas()
aLayer = canvas.currentLayer()

feat = QgsFeature()
iter = aLayer.getFeatures()

while iter.nextFeature(feat):
    aLayer.removeSelection()
    aLayer.select(feat.id())
    qgis.utils.iface.actionZoomToSelected().trigger()

    <wait for user input>

How can I do that?


texts: Cartodb - How to insert geojson data to table using SQL api?

texts:
I am currently using leaflet.js to generate polygons and I would like to add them to my table in cartodb. I.e. a geojson object :

{
  "type": "Feature",
  "properties": {
    "name": "wing1"
  },
  "geometry": {
    "type": "Polygon",
    "coordinates": [
      [
        [
          -60.1171875,
          54.16243396806781
        ],
        [
          5.9765625,
          70.0205873017406
        ],
        [
          41.484375,
          49.15296965617039
        ],
        [
          11.25,
          1.0546279422758869
        ],
        [
          -60.1171875,
          54.16243396806781
        ]
      ]
    ]
  }
},

I am aware I can use the following api to insert new data to cartodb table.

https://{username}.carto.com/api/v2/sql?q=INSERT INTO test_table (column_name, column_name_2, the_geom) VALUES ('this is a string', 11, ST_SetSRID(ST_Point(-110, 43),4326))&api_key={api_key}

But can someone explain to me how I can send a geojson object using the existing api?


texts: Pathfinder Differential Correction: "Use reference position from" who?

texts: How do you choose either "base files" or "base provider"?


texts: CARTO Builder, make calculation based on map point set by user

texts:
I'm have recent looking into carto in my search for a nice and powerful map analysis service.

I have set a few map nice, tidy and easy, but I'm missing an essential feature.
I'm looking on a way to make the data analysis based on input from the a users and where they click or placed a pin somewhere in the map.

It could be something like calculate the population with a 30 min isoline, with origin from where the user click on the map.

I have been unable to find this feature in the carto builder.

Can someone confirm if there is no such possibilities or point out to me where I can find it ?


texts:
I have a points layer and I want to create an interpolation grid with QGIS 2.0.
If I use Raster -> Analysis -> Interpolation the extent of the output raster is good (rectangular shape) but it doesn't reflect the original values of the points layer in input. The values of the raster are the result of interpolation calculations.

If I use Raster -> Interpolation the output file reflects the original values of the input points layer, but its shape is not what I'm looking for.

I don't see options to get the "original values" with the first method or to choose the extent of the area with the second one.
So, how can I choose the shape-extent of the interpolated raster and, at the same time, be sure that it will reflects the values of the input layer?


texts:
I want to convert spatial data from postgres to some basic formats (shp, geopackage, csv, xlsx, gml... ) using python gdal functionality. I read the answer to this question: Using ogr2ogr to convert GML to shapefile in Python?

I downloaded ogr2ogr.py and tried to convert data from postgres to some formats but I get an error saying:

KML driver failed to create /home/.../out.kml

This is the method I am using:

 ogr2ogr.main(["", "-f", "KML", "/home/.../out.kml",
              "PG:dbname='dbname' host='host' port=5432 user='username' password='password'",
              "-sql", "SELECT column1, column2 FROM schema.table"])


texts: R -plot 12 images per plot or save 12 images per image file in R

texts:
Here is R script and my question is
How can I plot 12 images per plot or save 12 images per image file because I have total 864 images.

RAD1998 <- raster(paste(getwd(), "/1998bil /1998ASC5min_ppt_spas1214_0001_19980202_0810_UTC.asc.bil", sep = ""))

list.ras <- mixedsort(list.files(paste(getwd(), "/1998bil/", sep = ""), full.names = T, pattern = ".asc.bil"))
list.ras
png('RAD1998.png')
RAD1998.all <- stack(list.ras)
plot(RAD1998.all)
dev.off()


texts: I have a DEM of a region and want to delineate about a half-mile region at the base of hills (e.g. an area where the slope changes from steep to shallow). I ran the slope analysis tool and it kinda works if I classify it with low slope values (0-10) and high (10+) (see image), but I am still left with the tops of ridges and high elevation areas also showing shallow slope. How might I just show a region at the base of hills?


texts:
I have a situation like this:

What I need to do is to connect each point to every line that is at most, let's say 200 m, away from the point. In other words, I need to draw a perpendicular line from each point to every line that is in the buffer.

Is there a way to do this in PyQGIS?


texts: Importing all layers open in QGIS to PostGIS

texts:
I'm looking for write a PyQGIS code to import all layers open in QGIS to PostGIS.

import os
#Code from  https://gis.stackexchange.com/questions/26257/iterating-over-map-layers-in-qgis-python
layers = QgsProject.instance().mapLayers().values()

for layer in layers:
  mytable=layer.name()
  path=layer.source()
  #Code from https://gis.stackexchange.com/questions/254182/saving-qgis-layer-to-postgis-using-pyqgis
  con_string = """dbname='Robanostra' host='localhost' port='5432' user='Daniele' password='******' key=gid type=MULTIPOLYGON table="s_500_patprova."mytable (geom)"""
  err = QgsVectorLayerExporter.exportLayer(path, con_string, 'postgres', QgsCoordinateReferenceSystem(3003), False)

The result is:

TypeError: QgsVectorLayerExporter.exportLayer(): argument 1 has unexpected type 'str'

I tried to convert str to path but I doesn't work. What am I doing wrong?


texts: How to include filename when importing txt / csv file in postgres?

texts:
When I import a csv-file in postgres using psql I would like to add the filename in a column (before each row). Is this possible?

The command I am using:

psql.exe -h host-U user -d database -c "\COPY schema.table FROM filename.csv with csv


texts:
I have recently switched to QGIS 3.4.0 and currently, I can't open projects saved in .qgz-format.

No error message, no crash, they simply don't open.

When I try to load the project a second time, I get a small tooltip at the bottom of the window "Project Loaded" - though that definitely is not the case.

On the other hand, .qgs-files can be opened without any hassle.

Is this a known error, maybe already solved?


texts:
Whilst working with OGR, I came across the following field used define a geographic dataset:

COORD_DIMENSION

The spatial dimension (2, 3 or 4 dimensional) of the column.

Whilst I've used 2D and 3D features, are there any examples of datasets in GIS that use the fourth dimension, and do any GIS systems actually handle these?

Update:

The 4d with regards to OGR / PostGIS is likely to be used for M-values (M-Aware in ArcGIS) (though I've yet to find this in a tech doc). However the quesion of (real?!) 4d in GIS remains open.


texts: Raster domain output polygon has is an empty feature class. Shape_Length and Shape_Area = 0

texts:
I am trying to convert a 32-bit floating point raster with 205289 columns and 120615 rows to a polygon using Raster Domain. It creates a feature class that does not contain any polygon.  There is a record in the attribute table but it says zero in the Shape_Length and Shape_Area fields.

Any ideas why I am getting this output? Or other suggestions how to convert this raster to a polygon? I just need a simple polygon to delineate the boundary of my raster.


texts: QSpatialite : Import multiple TXT / CSV?

texts:
I'm using QSpatialite only. I don't have spatialite-gui, neither the shell.

I want to import multiple TXT/CSV at one time, but I found nothing on the internet.

How is it possible to do this?


texts:
I am working on 3.12.0-București.
I have created buffers (radius=300 meters) on a set of sampling points.
I do not remember the number of segments that I have used. However, I have set a high number of segments in order to have a buffer resembling to a circular shape. I am working in the CRS WGS 84/UTM zone 32N - EPSG: 32632
I have tried to calculate the area of these buffers (in squared meters) using the field calculator $area (see figure).
I have noticed:

The value of the calculated area is different for each buffer even if I have created all the buffers using the same settings;

The calculated area is always slightly different from the area of a circle with radius = 300 meters: probably this is because the shape is not "strictly" a circle

See the coloumn of the area buffer (each row represents a different buffer).

Could you clarify these two points?


texts:
Using google maps I have located 4 points and have created a polygon with them in this way

insert into dummy_zip values (4,'prueba_directa','2016-04-13 08:08:15.973731','prueba_directa',1,1,100,'2016-04-13 08:08:15.973731',(SELECT ST_Transform(ST_GeomFromText('MULTIPOLYGON(((37.986504 -1.121951,37.986741 -1.121243,37.985794 -1.120792,37.985608 -1.121479,37.986504 -1.121951)))',4326),900913)));

I use the camptocamp example base_geoengine_demo for odoo (openERP 8), only can use the 900913, but i don't use google maps i use OpenStreetMap, if i use distinct 900913 i obtain a error "ERROR: Geometry SRID (0) does not match column SRID (900913)"

but when displayed on the map they are in a completely different location than where they should be?


texts: Computing area of different pixels classes in a multiple areas raster?

texts:
I work on forest data. My raster is a binary file : 0 is forest and 1 is non-deforest. I have multiple ROIs in which I would like to know the area of forest (so pixels of value 0 in the ROIs).
I managed to do it with r.report for a single ROI but I have a lot of ROIs and I would like to obtain it for all of them. Here is an example :

In the example I show, I would like to know area of black pixels in each of the ROIs, and not the sum of them as I obtain with r.report.


texts: Projection with Transformation using ogr2ogr GDAL

texts:
I am trying to reproject a shape from EPSG:3857 to EPSG:102008 (North America Equal Area Conic Albers) using GDAL.

I need to use the "NAD 83 to WGS 1984 5 (WKID 1515)" transformation since the datum are different.

My current command is using the -s_srs and -t_srs parameters.
i.e. -s_srs EPSG:3857 -t_srs EPSG:102008

How do I inject this particular transformation into the above command?


texts: How does one access a featurelayer in SDE via Python?

texts:
I am trying to use Arcpy to run the CopyFeatures_management script so that I can copy a featurelayer in SDE.

What do I use for the input (and output, for that matter, since I'll be copying the layer back to SDE) to access the layer?


texts: Using ArcGIS Raster Calculator to calculate MSAVI2 based on Landsat 8 OLI bands 4 and 5

texts:
I would like to calculate a Modified Soil-Adjusted Vegetation Index 2 on the basis of Landsat 8 OLI imagery, bands 5 (NIR) and 4 (RED), using the Raster Calculator. The formula should be:

I have tried to convert this using the following map algebra expression in the Raster Calculator:

MSAVI2 = (2 * "NIR layer" + 1 - SquareRoot((2 * "NIR layer" + 1) ^ 2  - 8 * ("NIR layer" - "RED layer"))) / 2

The returned values are between 0.5 and 58113.5. They should be between -1 and 1 though, so something is not right.

Any ideas about what is going wrong here?


texts: How to interpolate line Z from sparse attribute values

texts:
I have a series of end-to-end connected 3d lines which I'd like to add elevation values for all vertices on the lines which are currently null or 0. How to accomplish this?

In the following image there are 3 lines, with 1 selected and vertices marked (green squares). Of these 3 lines there 2 vertices which have a Z value. Each Z-bearing vertex is on a different line. How to calculate and assign the intermediate vertex Z values between A and B? Again referring to the image, the resultant point A.15 would be 661, A.14 is 662, ...

(For simplicity's sake this answer doesn't need to take X,Y distance between vertices into account. It's okay to simply divide A.Z & B.Z difference by number of points. Bonus kudos if that's included anyway ;-)

Both ArcGIS and QGIS/OGR solutions welcome; ditto ready to run tools and/or python.


texts: Intersection failing because of vertices of polygons

texts:
I've written some code to intersect features of one layer with another's, but one layer has polygons (the squares in the image) and the other has circles (not filled).

The problem is that I'm using IPointCollection interface to get the features' points and I'm only getting the border vertices of the squares

I tried getting the convex hull of one single feature but I failed because to use the intersect method I still have to convert to a IPointCollection

Here's the code I am using

// Get first layer's points
IGeometry multipoint = new MultipointClass();
multipoint.SpatialReference = m_firstLayerFirstfeatureIGeometry.SpatialReference;
IPointCollection multipointPoints = (IPointCollection)multipoint;

// Cycle all features of first layer and add them to the collection
foreach (DataRow m_dr in m_firstLayerfeatureGraphics.Rows)
{
    // Takes this feature's geometry
    ESRI.ArcGIS.ADF.Web.Geometry.Geometry adfFeature =
        m_firstLayerfeatureGraphics.GeometryFromRow(m_dr) as ESRI.ArcGIS.ADF.Web.Geometry.Geometry;
    // convert to a igeometry object
    ESRI.ArcGIS.Geometry.IGeometry m_tempfeatureIGeometry =
    ESRI.ArcGIS.ADF.Web.DataSources.ArcGISServer.Local.Converter.ToIGeometry(adfFeature, ctx);

    IPointCollection mtempfet = (IPointCollection)m_tempfeatureIGeometry;

    //Add copies of the polyline vertices to the multipoint.
    multipointPoints.AddPointCollection((IPointCollection)mtempfet);
}

ESRI.ArcGIS.Geometry.ITopologicalOperator6 mtopo1 = (ESRI.ArcGIS.Geometry.ITopologicalOperator6)multipointPoints;

// Get second layer's points
IGeometry multipoint2 = new MultipointClass();
multipoint2.SpatialReference = m_masterLayerFirstfeatureIGeometry.SpatialReference;
IPointCollection multipointPoints2 = (IPointCollection)multipoint2;

// Cycle all features of first layer and add them to the collection
foreach (DataRow m_dr in m_masterLayerfeatureGraphics.Rows)
{
    // Takes this feature's geometry
    ESRI.ArcGIS.ADF.Web.Geometry.Geometry adfFeature =
        m_firstLayerfeatureGraphics.GeometryFromRow(m_dr) as ESRI.ArcGIS.ADF.Web.Geometry.Geometry;
    // convert to a igeometry object
    ESRI.ArcGIS.Geometry.IGeometry m_tempfeatureIGeometry =
    ESRI.ArcGIS.ADF.Web.DataSources.ArcGISServer.Local.Converter.ToIGeometry(adfFeature, ctx);

    IPointCollection mtempfet = (IPointCollection)m_tempfeatureIGeometry;

    //Add copies of the polyline vertices to the multipoint.
    multipointPoints2.AddPointCollection(mtempfet);
}

ESRI.ArcGIS.Geometry.ITopologicalOperator6 mtopo2 = (ESRI.ArcGIS.Geometry.ITopologicalOperator6)multipointPoints2;

// INTERSECTION
ESRI.ArcGIS.Geometry.IGeometry mpointCol = (ESRI.ArcGIS.Geometry.IGeometry)mtopo2.IntersectEx((IGeometry)mtopo1, false, esriGeometryDimension.esriGeometryNoDimension);

// Draw the result
DrawIGeometryOnMap(m_Map, mpointCol);


texts:
My ArcMap 10.6 Python code for creating a specific raster isn't working. The code is fine but the output keeps saying failed to draw.

import arcpy
from arcpy.sa import *
arcpy.env.workspace = r'C:\Users\npenme2\Desktop\Faunalyze\Right_Femur.gdb'
rasters = arcpy.ListRasters()
(u'Feature_hul',u'Feature_KB3_Left_Humerus_14502_RC1',u'Feature_KB4_Left_Humerus_14515_RC1',u'KB1_Left_Humerus_14476_RC',u'KB2_Left_Humerus_14489_RC',u'KB5_Left_Humerus_14528_RC',u'KB5_Left_Humerus_14528_RC_zeronew1')
r1 = sum(Raster(r) for r in rasters)
r2 = sum(Con(r, 1, 0, 'Value > 0') for r in rasters)
avg_rast = FloatDivide(r1, r2)
avg_rast.save('Total_Left_HUmerus')

Also it's taking an extremely long time to process.


texts: Setting raster visibility on/off in QGIS Python API 3

texts:
I'm having a bit of a trouble migrating an algorithm from pyqgis 2 to pyqgis 3. The thing is I want to make visible and invisible a series of layers as I iterate through a list and render some images.

A little more explanation: I have this big python list which contains several dictionaries. Each item of the list contains a series of "actions" I need to perform in order to render some images with a specific configuration. One of this possible actions is turning on or off a layer's visibility (and eventually the Google Maps Raster). So, in my code what I'm trying to implement is some function like this:

def __layer_actions(actions):
    for action in actions:
        if action["type"] == "layer_opacity":
            get_layer(action["layer"]).setOpacity(action["opacity"] / 100)

        if action["type"] == "layer_visibility":
            if(action["visibility"] is True):
                get_layer(action["layer"]).draw()       # fails when the layer is actually a raster

Since the Raster's draw function is not the same as the layer's one, when I'm processing a Raster I get this error:

TypeError: QgsRasterLayer.draw(): not enough arguments

How can I rearrange things to contemplate this case in which my code is failing?
For example, if I could just manage to turn the Raster's visibility on/off somehow, I could use a try/except sentence.


texts: Update postgres/postgis point data attribute from another enclosed polygon shapefile in QGIS

texts:
I have

a placename point table in postgis with 'district' column.
a districts polygon shapefile with district name in 'Name' column.
See the image below.

I would like to update the 'district' columns(in placename psql table) with the corresponding district name from the shapefile polygons in which each placename point feature reside.

How can I do this in QGIS?
Or Is there any other way?


texts:
Our company just switched over to a new server setup and at the same time moved from ArcMap 10.1 to 10.2. Everything seems to be working ok but when I am working with feature classes I sporadically get this error (in this case alt2zon is my feature class)

One or more layers failed to draw:
alt2zon:  General function failure

It doesn't seem to be tied to any one action I am performing, although it seems to occur when I am running larger datasets. If I remove the layer and re-add it still doesn't work. To solve the problem I have to save the MXD, close it, reopen, and the issue is resolved.

Any ideas as to what is happening?


texts: what is the difference between topo to raster and points to raster tools in arcgis

texts:
I am working with point data, and I've found two ways to create raster data.

1 Topo to raster

2 points to raster

the values of each result are different, and here are some pics for reference.
I want to know what are the differences for those two tools and their applications.


texts: Populate legend items with values from another layer

texts:
I'm using QGIS 3.10.
I'm trying to label my map with values taken from a different layer (meaning not my atlas coverage layer), which is intersecting the atlas cover layer.
This layer named "Colorisé" has a categorized style using the field "color_id". I want my legend items to be labelled with the field "idprop".
I have tried using this expression (1), which returns the featured items not individually but as a list :
aggregate('Colorisé', 'concatenate_unique', "idprop", intersects($geometry,@atlas_geometry), ',', "idprop")

And I have tried using the expression below (2), which returns the first value it can find (I guess) and applies it to all my legend items :
attribute(get_feature('Colorisé', 'color_id', intersects($geometry, @atlas_geometry)), 'idprop')

Can you guys tell me if I'm anywhere close to getting the desired results ?
You can find below an excerpt from the table of attributes of the layer "Colorisé" and an image summing up the situation.

EDIT
Using this expression seems to be working in the preview, but the labels end up being empty :
attribute(get_feature(@layer,'color_id',
attribute($currentfeature, 'color_id')),'idprop')


texts: Adding MySQL to Geoserver for Non-Spatial Data

texts:
I need to connect to a MySql table, which has all the attribute values and the lat & long. I am using a GeoServer Suite to host my raster and so felt the best way to go about it would be to add the MySQL to my Geoserver, but it seems, geoserver comes with Postgres by default.

Is there any My SQL extension to geoserver?

Is this the best approach to obtain the info from a MySql database?

UPDATE:
Ok so what I am trying to do is eventually migrate my Database to PostGIS. So I am trying to create a shapefile from MySQL database. The database has non-spatial data with the lat and long associated with each attribute data. I found pyshp, but i think that is reading the shape file. Can we create shapefile using the pyshp (meaning can we create shapefile from MySQL database using pyshp)? Becoz if i am able to create the shapefile then I can easily insert it into PostGIS using geoServer and it can be used in other softwares as well.


texts: Flow Direction tool outputs strange raster that doesn't work with subsequent hydrology tools

texts:
My main goal here is to end up with a shapefile that simply defines the ridgeline in a given area. I followed the directions at How To: Identify ridgelines from a DEM and am having some trouble. Firstly I went through every step exactly as it is shown on that page, and ended up with a shapefile that defined almost every cell in the raster as "ridgeline." I went back through in order to diagnose where I went wrong. What I found was that in the step where you invert your filled DEM, the resulting raster has all negative values for elevation. I modified the equation in the raster calculator to be:

(("Rasterlayer" - Max_value) * -1) + Min_value

This made it so that my inverted and filled DEM had positive values. I threw the two DEMs (inverted and not inverted) into arcscene just to check my sanity and yes, the new inverted DEM showed valleys where previously there were ridges. All made sense.

However, when I then move on to the flow direction tool, I get weird results. I specify the new inverted DEM as the input, and choose to save the output as the same name plus "flow_dir" somewhere in the name. I checked the box to make the edge cells flow outward, and I left it at D8, which is what I wanted. In theory, this should give me a raster where the cells are defined as one of 8 discrete values, 1, 2, 4, 8, 16, 32, 64, & 128. In practice, however, I get a raster that has not 8 discrete cell values but 208 discrete values, from 1-208. What have  done wrong here and how do I fix it?

It should be of note that when I run the original filled DEM (non-inverted) through the flow direction tool with the exact same parameters, it spits out a raster that has the expected 8 discrete values. The screen grab here shows the raster that I get with 208 values (Flow_dir_7) which is not what I want but was resultant of using the tool on the inverted DEM. As you can see, right below that is the Flow_Dir_5 raster which is defined by the 8 values as it should be for D8, however this raster was produced using the non-inverted DEM as the input.

Do you know where I went wrong or know another way to define ridgelines other than inverting my DEM and pretending they are streams?

I'm running Windows 10 and ArcMap 10.7.1.


texts:
I was wondering if it's possible to use the standard auto increment code to work on groups of records:

FeatureClass = "newFC"
recs = collections.defaultdict(int) autoIncrementEXP = "autoIncrement()"
expCalc = """def autoIncrement(!Group!): recs[!Group!] += 1 return recs[!Group!]"""
arcpy.CalculateField_management(FeatureClass, "Number",autoIncrementEXP,"PYTHON",expCalc)

So instead of having the auto increment portion of the script stop when it reaches the end of the records, have it create sequential numbers using field value groups if that makes sense. You would start with something like the Group field below and end with the Number field containing the increment values:

Group   Number
A       1
A       2
A       3
B       1
B       2
C       1
C       2
C       3
C       4
C       5

Thanks for any insight you can provide.


texts: Rasterize shapefile according to another big raster using gdal_rasterize

texts:
I want to rasterize a shapefile into a raster according to another raster (97284, 154208, 15001971072  (nrow, ncol, ncell)).

Since the raster is very big, I am using gdal_rasterize following the code below

    rp<-gdal_rasterize(shapefile, big raster)

However, I am getting this error

Error in as.character.default() :
  no method for coercing this S4 class to a vector

1) What format does the shapefile need to be in to rasterize?
2) Is there another way to rasterize a shapefile? I have tried rasterize() and it is taking very long and I have also read posts that say that gdal_rasterize is much faster for very big raster


texts:
I want to implement day/night functionality on world map i.e there will be
  dark layer showing countires having night now.

This world map could be creaeted through ESRI or SVG map.I want to understand what is the
  logic of showing behind implementing this.Any known algorithms/formulas.

I have posted this question in stackoverlfow (https://stackoverflow.com/questions/19685010/implementing-day-night-functionality-on-map).

But thought this is better forum.


texts: How to create a routable network with PGRouting

texts:
I have two sets of data

Table A: Stations - contains postal code (which I can convert to lat/long), name and code.
Table B: Connections - contains a possible list of connections including source, target (both expressed as three letter codes), distance and maximum speed (I can calculate a cost from these two parameters).

I want to create a routable network - from looking, I have found PGRouting may be able to do what I need but almost all the tutorials I've found are based on road networks.

What is the best way to convert the data I have into a routable network? I'm not set on PostGIS/PGRouting but I haven't found a better solution yet - eventually there'll be around 3000 stations and 5000+ connections.


texts: Values disappear after joining .csv file to shp layer in QGIS

texts:
I am joining a .csv file to shp layer and it goes well except all the values from the .csv file disappear (there's NULL in the attribute after the join)?

I don't have X & Y coordinates, but I'm combining it with a place name.

For example, shp files contains text columns (example made in excel on purpose):

And the .csv file contains also additional info

I also tried to add these columns straight into the shp layer attribute table to fill in the info there from the csv file by hand, but I couldn't fill in any information. There's also NULL (in both form view and table view options), but I can't edit the cell at all. I can only edit the cells that already existed in the shp file itself.

What could be wrong with this as well, should I enable editing separately for these new columns?

I'm quite a beginner with QGIS, but this is really something I just can't figure out.


texts: CartoDB Visualization refresh in legends does not work until you switch template

texts:
When changing the settings through the wizard, for example of the Bubble visuatlization, the legend does not update automatically.

I managed to update the legend by changing the template to 'custom' then back to 'bubble'

Is there another way to make it update?


texts: Connecting to ArcGIS Desktop PostgreSQL database (running locally)

texts:
I know this is a total noob question but then again, when it comes to PostgreSQL I am a total noob...

I have installed the OpenGeo suite on my computer, running Windows 7, 64-bit.  I have downloaded the tutorial data and have been working through it.  I would like to access the database that I have made with ArcGIS so I can start learning about how that works.  However, when I try to connect to the database using the "Add Database Connection..." I can't seem to figure out what to put in for 'Instance' - everything I've tried (localhost, localhost,54321, my IP address with and without the port, my computer name) doesn't seem to work.

I installed and created the nyc test database according to all the instructions in the OpenGeo tutorial pages but am at a total loss here.  I know this must be an easy thing to deal with, I just can't get my head around it.  Can someone with more experience throw me a bone on this one?

Frustratingly enough, QGIS connects to the database with almost no effort and everything works fine - but my organization is heavily ESRI focused and switching to QGIS is a long shot...

EDIT
Thanks everyone for pitching in on this.  I really appreciate the advice, it's what makes this site great.  I have installed

OpenGeo Suite 3.0.1 with no extensions
ArcGIS Desktop 10.1, no license for server

I just went to the Opengeo website and downloaded the Windows installer and installed the default configuration. I have tried every permutation for the instance name that I can think of - using colons and commas, my machine name, localhost, postgresql, server, my IP address, random curse words, anything I can think of.  I have been using port 54321 instead of 5432 because that is what the OpenGeo workshop told me to set up my 'nyc' practice database to use.

It seems that localhost should work, here is a shot of the server properties from pgAdmin:


texts: Disabling annoying Vertex Editor Panel in QGIS

texts: Is there a way to disable the Vertex Editor panel popping up, saying "Right click on an editable feature to show its table of vertices"? This panel should be visible on demand, i.e. when you actually need to display and edit coordinates in a table but not by default. It's really annoying.


texts: How to cluster data (and not only the number of markers)

texts:
I wanted to display some data (a number, different for each point) on a map (either Google maps or OpenStreetMap). Each point has an associated value (the number above).

For markers which are close to each other, the display is messy and the natural solution for this is to cluster the points.

While there several libraries or built-in mechanisms to cluster data points (including OpenLayers which seems to be very complete), all the ones I found display the number of markers in a cluster, without the possibility to display a sum of values carried by each point (or an average, or other functions which takes values an an input and transform them into one value).

Is there a library which would help with this problem? (in JavaScript for pure client-side rendering and ideally open source)


texts: Setting simple definition query using ArcPy?

texts:
I am trying to set a definition query on a Application Lands layer in my project.mxd using arcpy. But it's not working for me. When I launch the project the query is in there but it doesn't get applied.

Can anyone help me out here?

import arcpy
mxd = arcpy.mapping.MapDocument(r'C:\Users\Desktop\GIS PROGRAMMING\Project.mxd')
Dockets = arcpy.mapping.ListLayers(mxd, "Application Lands")
Applications= Dockets[0]
Applications.definitionQuery = '"Docket"' + " = " + "'1254-2017 '"
mxd.save()
del mxd


texts: How to handle 500 Internal Server Error with ESRI Silverlight API

texts:
I have a simple AGS map service. I accidentally specified an invalid SRID explicitly in my Silverlight XAML:

<esri:Map x:Name="Map" IsLogoVisible="False">
    <esri:Map.Extent>
        <esri:Envelope XMin="526801.81" YMin="489775.0655" XMax="593105.41" YMax="523904.2045">
            <esri:Envelope.SpatialReference>
                <esri:SpatialReference WKID="3807" />
            </esri:Envelope.SpatialReference>
        </esri:Envelope>
    </esri:Map.Extent>
    <esri:Map.Layers>
        <esri:ArcGISDynamicMapServiceLayer x:Name="MyLayer" ID="MyService" Url="http://localhost/ArcGIS/rest/services/MyService/MapServer" Initialized="MyLayer_Initialized" InitializationFailed="Layer_InitializationFailed" />
    </esri:Map.Layers>
</esri:Map>

This results in a RESTful request that looks like this:

http://localhost/ArcGIS/rest/services/MyService/MapServer/export?bbox=526801.8125,479254.134305427,593105.4375,534425.146944573&size=1078,897&format=png24&transparent=true&layers=show:0,1&f=image&imageSR=3807&bboxSR=3807

This request fails with the following IIS error:

Error Summary

HTTP Error 500.0 - Internal Server Error
The page cannot be displayed because an internal server error has occurred.

Detailed Error Information

Module  ESRI.RestModule-Integrated
Notification    ResolveRequestCache
Handler StaticFile
Error Code  0x00000000

If I replace the invalid SRID, everything works fine.

Is there any way to get the Silverlight API to report exceptions like this? As it is I have to use Fiddler to find it.


texts: Seeking documentation for {z}/{x}/{y} syntax in tilesource urls?

texts:
Speaking of tilesource urls such as

http://{s}.tile.stamen.com/watercolor/{z}/{x}/{y}.png

is there any technical documentation or canonical reference for the curly braces syntax?

I'd specifically like to know what other things can be used in curly braces, and in which places, and how they work.


texts: Building Custom Mosaic Dataset Overviews in ArcGIS

texts:
We'd like to resample our mosaic dataset using the MAJORITY option, but it it isn't given in the Build Overviews tool. Our plan is to generate the rasters manually using Resample and Split GP tools, and then register them with the mosaic dataset.

Once we've created the TIFFs that we're going to use as overviews, how do we register them as such within the mosaic dataset? How do we let ArcGIS know that these rasters should be used as overviews?


texts: Storing metadata of various services according to ISO 19119?

texts:
I need to store the metadata of web services according to ISO 19119 standard.

Hence I am looking for some guidelines for designing the database for storage.

Can any one suggest any document to help in writing the SQL script


texts: Ask for polygons but getting points with QuickOSM in QGIS Modeler

texts:
I am attempting to download both points and polygons using QuickOSM+modeler in QGIS 3.10 (also tried QGIS 3.8). It's downloading the points fine, but when I attempt the same thing with multipolygons it is still giving me points.

Here's a look at the full model:

And here's the String Concatenation tool:

The log output gives some clues as to what might be happening:

URL: 'http://www.overpass-api.de/api/interpreter?.....&info=QgisQuickOSMPlugin\'|layername=multipolygons\''

This made me think that single quotes were being included as part of the string, so I tried again removing the single quotes. The output log now read:

URL: 'http://www.overpass-api.de/api/interpreter?.....&info=QgisQuickOSMPlugin|layername=multipolygons'

But, the result was the same I am still only getting points and not the building outlines of the shops that I seek.


texts: Loading a qgsproject when another one is loaded

texts:
I am loading a project with:

QgsProject.instance().read(QFileInfo(path))

When i try to load another project it subtracts the layers and shows nothing. But if I try doing it again it works.


texts:
Can someone please tell me how will I discard the clouds from Landsat 8 images using R? I have used the code

LC08_20151128_new <- cloudMask(LC08_20151128, threshold = 0.8, blue = "layer.1", tir = "layer.9", buffer = NULL, plot = FALSE)

I read that for landsat 8 layer 1 and layer 9 might be useful to detect the clouds. Is it correct? I got the image below after runninng the code

And then I masked the image with the cloudMask

LC08_20151128_crop <- mask(LC08_20151128,LC08_20151128_new, maskvalue = NA)

But it seems that the code retains the area with the cloud cover instead of removing it. Please see the RGB image below


texts:
I have two shape files: road layer and a district layer. Is it possible to know which roads belongs to which districts?

I'm currently using QGIS.


texts: How to use GDAL RasterizeLayer (python binding) -where flag

texts:
I am trying to use GDAL RasterizeLayer python binding, but I can not find a way to use the -where flag.

The 'pure' GDAL command I am trying to run is:
gdal_rasterize -b 1 -burn 0.00001 -where 'DN=255' -l pilot_1class.lazground_poly pilot_1class.lazground_poly.shp testdiff_again.tif

I am trying to use this in python:
gdal.RasterizeLayer(target_ds, [1], layerPoly, None, None, burn_values=[0.0001], options = ['ALL_TOUCHED'])


texts:
I am wanting to get a copy of the OS 1km Grid layer.

I have found these previous questions:

Ordnance Survey grids - download or generate?
and
Ordnance Survey Grids

But neither provide the information I required.

The first question does provide some links to the OS web site but these now return the 404 error. A subsequent comment mentions its half way down this page:

https://www.ordnancesurvey.co.uk/opendatadownload/products.html

But which specific item should I be downloading? It is not obvious!


texts:
I am using QGIS 3.6.3.

I am trying to rasterize 5 different maps from the same attribute table but it seems 3 of the maps produced have the same value even though they have different attribute values.


texts: Setting waste water gravitational flow direction using DEM automatically?

texts:
I have a Geometric Network of Waste Water and I need to set line digitize direction according to the gravitational flow direction using a DEM.

I'm trying to avoid flipping manually 1000 Km of pipes...

I'm well acquainted with ModelBuilder but a beginner with python so any hint will help!


texts: Get list of fields in empty table using PyQGIS

texts:
I have seen many examples of how to get the list of fields, but almost all using the layer that is already on the map, what I need is to create a new layer but from the structure of one that exists in PostGIS that when it has elements there are no problems, but when this is empty it is not possible to obtain the list of fields of its structure.

This is part of the code

def AddPostgisLayer_to_gpk(self,esquema,esquemagrupo,loadlayer, tipo,strtipo,nomenclatura):
    uri = QgsDataSourceUri()
    urilistacampos = QgsDataSourceUri()
    GpName = 'd:\\temp\\' + nomenclatura + '.gpkg'
    options = QgsVectorFileWriter.SaveVectorOptions()
    urilistacampos.setConnection("localhost", "5432", "bdsapromat", "postgres", "postgres",
                                 QgsDataSourceUri.SslDisable)
    urilistacampos.setDataSource(esquema, loadlayer, "geom",  "geocodigo like " + tipo)
    urilistacampos.setSrid("4326")

   vlayerlistfields = QgsVectorLayer(urilistacampos.uri(), loadlayer, "postgres")
    vlayercreada = QgsVectorLayer()

    # Get its list of fields
    inFields = vlayerlistfields.dataProvider().fields()

        if strtipo == "_areal":
            vlayercreada = QgsVectorLayer('MultiPolygon?crs=EPSG:4326', loadlayer, 'memory')
        if strtipo == "_lineal":
            vlayercreada = QgsVectorLayer("MultiLinestring?crs=epsg:4326", loadlayer, "memory")
        if strtipo == "_puntual":
            vlayercreada = QgsVectorLayer("MultiPoint?crs=epsg:4326", loadlayer, "memory")
        if strtipo == "_textos":
            vlayercreada = QgsVectorLayer("MultiPoint?crs=epsg:4326", loadlayer, "memory")

        # Copy the fields  into the new layer

        adiciono = vlayercreada.dataProvider().addAttributes(inFields.toList())
        vlayercreada.updateFields()
        QgsProject.instance().addMapLayer(vlayercreada, False)

    if vlayercreada.isValid():
        options.layerName = esquemagrupo + "." + loadlayer + strtipo
        if os.path.isfile(GpName):
            options.actionOnExistingFile = QgsVectorFileWriter.CreateOrOverwriteLayer

        options.EditionCapability = QgsVectorFileWriter.CanAddNewLayer
        error = QgsVectorFileWriter.writeAsVectorFormat(vlayercreada, GpName, options)

I am using QGIS 3.10.3.


texts: Getting base GIS files for use with pgRouting library

texts:
I'm just starting to read the tutorial for pgRouting and the first issue that troubles me is the origin of the "base map" to feed pgRouting with. I'm aware that there's available a tool called osm2pgrouting, which converts OpenStreetMaps data into pgRouting compatible ones. But I'm not sure if it is the correct/best source for this info, although it clearly might be a starting point and is fairly easy to get it.

What are the other possible source options for this purpose?

I work in an AVL company so, we gather second by second lots of GPS data from vehicles, which we don't actually process to "connect the dots" and turn them into streets, but I'm not sure if this actually helps or not.

Data needed for this project would be for Chile and, as starting point only for the city of Santiago.


texts: How to create and display a shakemap in openlayers?

texts:
OpenLayers is a library that develops tiled maps and shows heatmap on map, it also provides clustering. The examples are here.

I tried with this example and tried to customize it for Shakemap. Below is example of expected output.

This above map is USGS map.

I used following function to customize the example given by OpenLayers.

   function get_gradient() {
    var gradient = {};
        //colors selection
        var colors = ["#FFFFFF", "#BFCCFF", "#A0E6FF", "#80FFFF", "#93FF7A", "#FFFF00", "#FFC800", "#FF9100", "#FF0000", "#c80000"];
       //some logic of selection of colors
     //  gradient[1] = colors[1];
        return colors;
        //return gradient;
    }

I gave this as input to vector layer's gradient attribute

  var vector = new ol.layer.Heatmap({
       source: new ol.source.Vector({
         url: 'data/kml/2012_Earthquakes_Mag5.kml',
         format: new ol.format.KML({
           extractStyles: false
         })
       }),
       blur: parseInt(blur.value, 10),
       radius: parseInt(radius.value, 10),
       //gradient : ["#FFFFFF", "#BFCCFF", "#A0E6FF", "#80FFFF", "#93FF7A", "#FFFF00", "#FFC800", "#FF9100", "#FF0000", "#c80000"]
       gradient : get_gradient() //here
     });

I want some business logic that gives me idea of what filter should be applied to colors.

This is done using OpenLayers 4.


texts: How to get the center line between two linestrings

texts:
Where is the track of a ship when it is shipping always in the center of the English Channel  between the coasts of Great Britain and France?

I started with 2 linestrings and I created a lot of points on the linestring with a postgis query.
After that I created a distance matrix with Raster/Analysis/Proximity in QGIS.
Finally I created an aspect raster from the output raster.

Is the red line really the "center line" between these two linestrings?

Is there a more efficient way to do this?


texts: Rotate a Leaflet FeatureGroup with finger touch

texts:
I have some rectangles that I store in a leaflet FeatureGroup. Currently, if I define the FeatureGroup, it appears parallely to the screen. I would like to only rotate this FeatureGroup on top of the Leaflet map, ideally by touch commands. For example, you can touchselect that FeatureGroup and then rotate it with 2 fingers. Otherwise, rotation by a number value like a rotate-button with 360° is fine as well.

Any ideas?


texts: Dynamic Label or Annotation in QGIS canvas

texts:
I'm creating an animation from several Raster images, and I'd like to add a label or an annotation to the map canvas, not the print composer, that will enable me to show the layer's name, which holds the date.

I'm assuming that I'll need to use the Form Annotation, but all examples that I've seen are used in a vector layer.

The end result should be a label that shows the visible raster layer(s) in the TOC.

How do I go about doing that?


texts: QGIS 2.0 Plugin error: Unable to load GdalTools plugin?

texts:
I just used the osgeo4w exe to install QGIS and all the other software included in the basic setup.

When I open QGIS, I get the following error:

Plugin error
Unable to load GdalTools plugin.
The required "osgeo [python-gdal]" module is missing.
Install it and try again.

I opened the Python.exe in the osgeo4w bin folder and am able to import gdal just fine.

Any idea what's wrong here?


texts: Install unstable Qgis-master alongside release version

texts:
Good day
I have the release version of Qgis (2.16) installed into Mint and it can be updated just with apt install qgis whenever a notification comes out that a new version is released.

I want to also install the latest master (2.99) alongside 2.16 for testing but I don't want 2.16 to be overwritten in case there is a problem with master.

Unlike the Windows scenario where you can install different versions of the osGeo4W package to have several versions of Qgis, I haven't been able to work out a way of just creating a 2nd Qgis installation alongside the existing one. I haven't been able to find a simple install package I can download, just repositories, but I haven't been able to work out how to get the master version from them.


texts: Pyqgis: Some way to import shp(shapefile) to plugin(to create a toolbar of shpfileS), when we want a use a shpfile to active or deactive

texts: Some way to import shpfile(one shp file) in plugin, not go over the standard procedure(Select shp file, Toogle editing and Add features)-to do it with python, because I am trying a make a custom toolbar. Any solutions?


texts:
EDIT: OK, added two lines to the code below including a projectAs().

When I start my mxd and run this code it selects a polygon. However, if I then pan to a new location - still within the bounds of the target FC - it doesn't work. Is there some update command I need to run after panning? It's as though the dataframe extent is lost after panning.

====

I'm trying to selectByLocation using a pointGeometry. The code below creates a pointGeometry from the center of the current dataframe extent. It works if the dataframe and the target feature class (polygon features) are in the same coordinate system but not when they're in different systems.

My usage model is that the target will usually be in a continental projected system or geographic coordinates while the dataframe will be in NAD 1983 UTM 10N to 20N.

import arcpy

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]

arcpy.env.outputCoordinateSystem = df.spatialReference
**sr = arcpy.Describe('targetFC').spatialReference**

myPt = arcpy.Point((df.extent.XMax + df.extent.XMin)/2.0, (df.extent.YMax + df.extent.YMin)/2.0)
myPtGeometry = arcpy.PointGeometry(myPt, df.spatialReference)
**myPtGeometry.projectAs(sr)**
arcpy.SelectLayerByLocation_management('targetFC', 'INTERSECT', myPtGeometry)

del df, mxd

Thank you,
Scott


texts: OSM Hebrew caption are left to right instead of right to left in QGIS

texts:
Whew I am usieng OSM maps in QGIS the Hebrew caption are in the right order, right to left
However, when I am downloading the streets with "OSMDownloader InaSAFE" plugin the words are left to right. I am using UTF-8 encoding and the caption in the Attribute Table are in the right order.

How can I solved this and display this Hebrew caption?


texts: I am trying to create a polygon feature from an enclosed line feature using the Geometry Tool: Lines to Polygons. When I do this, I get several spikes with the newly created polygon. I then run the "Check Geometry Validity" tool and it notes that I have about seven duplicate nodes. I  am trying to delete the duplicate node using the Edit tools: Node Tool, select the node, but have not been successful with deleting the extra (or duplicate) node. I was told to use Fn + Delete on the keyboard (Mac OS), but nothing seems to happen. Please tell me the best way to get rid of these duplicate nodes.


texts: Polygon intersection with bounding box with polygon's holes intact

texts:
I am trying to serve GeoJSON FeatureCollection tiles with each feature's geometries trimmed to the tile boundary.

I am working on a Leaflet plugin to re-assemble the geometries. It works by building up a pool of lines for a unique feature, removing any equal lines in opposing directions, then walking any unwalked lines to recreate the linear rings.

This currently works, except I am having trouble dealing with polygon holes. When I have a polygon with a hole and generate the GeoJSON to intersect it with the tile boundary - I lose the hole information.

Does anyone have any ideas on how I could keep the hole information intact on the server side, or alternatively detect this happening on the client-side?


texts: Importing GPS bird tracker text file and i want to display attributes

texts:
I am importing GPS bird tracker data with a cut down file consisting of

Date/time.
Lat/Long.
Speed.

I would like the points to display with the speed as text.  Is this possible, if so how?

I'm using QGIS Madeira (3.4).


texts:
I'm trying to use ogr2ogr to extract some attributes from a Shapefile to a netCDF. I can limit which attributes are read using the -select argument, but it seems that the geometry is written to the output file no matter what. Is there some way to exclude the geometry using ogr2ogr ?

ogr2ogr -f "NetCDF" /tmp/ogrout.nc \
        -dsco "FORMAT=NC4C" \
        -select "HYBAS_ID, NEXT_DOWN" \
        /tmp/basins.shp


texts: Circles projected from AEQD to WGS84 too elongated

texts:
I'm trying to draw circles around lat-long (WGS84) points by:

transforming the point into AEQD space
creating a shapely buffer on the projected point
projecting the shapely buffer back to WGS space

```

import pyproj
from shapely.ops import transform as shapely_transform
from shapely.geometry import Point

from functools import partial

def non_euclidean_circle(center, radius):
    lat, lon = center
    # proj4str = '+proj=aeqd +lat_0=%s +lon_0=%s +x_0=0 +y_0=0' % (lat, lon)
    AEQD = pyproj.Proj(proj='aeqd', lat_0=lat, lon_0=lon, x_0=0, y_0=0)
    WGS84 = pyproj.Proj(init='epsg:4326')

    # transform the given lat-long onto the flat AEQD plane
    tx_lat, tx_lon = pyproj.transform(WGS84, AEQD, lat, lon)
    circle = Point(tx_lat, tx_lon).buffer(radius)

    def inverse_tx(x, y, z=None):
      return pyproj.transform(AEQD, WGS84, x, y)

    # inverse projection from AEQD to EPSG4326-WGS84
    return shapely_transform(inverse_tx, circle)

c = non_euclidean_circle((42, -72), 100)
print(c, c.area)

When I plot the polygon c on a map, it looks like:

Complete plotting code: https://codepen.io/anon/pen/yEoJrG?editors=0010

My question is, is this correct? Are these points on the manifold of the Earth actually equidistant from the center, or am I doing something wrong here?


texts: Cut one polygon into smaller polygon with other polygons

texts:
I want to use one layer with numerous polygons in it (layer 2 in attached image) to cut a big polygon on another layer (layer 1 in image), so that said big polygon is separated into smaller polygons.

For example, these are drainage areas that run across a landscape, but the different sections in between these drainage areas have different sensitivities and must be classified differently. That is the aim of this exercise. Hope it makes sense!

I have been unable to figure how to do this with either standard vector tools or with plugins (like clipper or digitize tools).

Any ideas?


texts: Merging multiple attributes of a raster into a single attribute using R?

texts:
I have a raster that has 50 land cover classes. I am trying to combine multiple attributes from the attribute table into a single attribute. For example, I would like "Evergreen", "Semi-Evergreen" as - EVERGREEN.

I used the Editor Toolbar in ArcGIS and selected these features, but the Merge button on the Editor toolbar is greyed out. I read another solution that asked me to "delete the Shared Features tool on the Topology toolbar", but I was unable to do so. Apparently, deleting the "Shared Features" tool will activate the merge tool?

Either way, I would like to combine the attributes into a single attribute in R.

What I have tried in R so far (just loading and viewing the attributes)

> lulc_2015 <- raster("E:\\MasterThesis\\Absence Approach\\Data\\LandCover_Data\\2015_lulc")
> lulc_2015
class       : RasterLayer
dimensions  : 26782, 11768, 315170576  (nrow, ncol, ncell)
resolution  : 60, 60  (x, y)
extent      : 163125, 869205, 855611, 2462531  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=43 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
data source : E:\2015_lulc
names       : X2015_lulc
values      : 1, 134  (min, max)
attributes  :
    ID     COUNT CLASS_NAME
from:   1 118944856
to  : 134   3589862 Settlement

> lulc_2015@data@attributes
[[1]]
ID     COUNT                       CLASS_NAME
1    1 118944856
2    2   2397462                        Evergreen
3    6     39397         Sub-tropical broadleaved
4   11   2457007                   Semi-evergreen
5   12   8798892                  Moist deciduous
6   13       184        Sal mixed moist deciduous
7   14   1076081       Teak mixed moist deciduous
8   15   6675337                    Dry deciduous
9   17   1790087         Teak mixed dry deciduous
10  18    476694                     Thorn forest
11  19    247793                     Bamboo mixed
12  20        16             Temperate coniferous
13  22      1109                              Sal
14  23    708833                             Teak
15  25    138858                           Bamboo
16  26         4                             Pine
17  30     85417                       Hardwickia
18  33      2179                        Boswellia
19  35       210               Anogeissus pendula
20  41    107796                         Mangrove

How do I combine these classes?

For those providing a possible solution, please note I am currently using Reclassify tool from the Spatial Analyst extension in ArcGIS to achieve the same, but if there is an alternative, please provide the same below.


texts: Removing black region in published ECW file in GeoServer?

texts:
I had published an ECW file in geoserver 2.11.2 with SLD generated from QGIS. After publishing the layer i could see some black region in layers preview with the SLD generated from QGIS as below. This is the SLD generated from QGIS

<?xml version="1.0" ?>
<sld:StyledLayerDescriptor version="1.0.0" xmlns="http://www.opengis.net/sld" xmlns:gml="http://www.opengis.net/gml" xmlns:ogc="http://www.opengis.net/ogc" xmlns:sld="http://www.opengis.net/sld">
<sld:UserLayer>
    <sld:LayerFeatureConstraints>
        <sld:FeatureTypeConstraint/>
    </sld:LayerFeatureConstraints>
    <sld:UserStyle>
        <sld:Name>SAD61_Tamale_Satellite_ECW25_20170114</sld:Name>
        <sld:Title/>
        <sld:FeatureTypeStyle>
            <sld:Name/>
            <sld:Rule>
                <sld:RasterSymbolizer>
                    <sld:Geometry>
                        <ogc:PropertyName>grid</ogc:PropertyName>
                    </sld:Geometry>
                    <sld:Opacity>1</sld:Opacity>
                    <sld:ChannelSelection>
                        <sld:RedChannel>
                            <sld:SourceChannelName>1</sld:SourceChannelName>
                        </sld:RedChannel>
                        <sld:GreenChannel>
                            <sld:SourceChannelName>2</sld:SourceChannelName>
                        </sld:GreenChannel>
                        <sld:BlueChannel>
                            <sld:SourceChannelName>3</sld:SourceChannelName>
                        </sld:BlueChannel>
                    </sld:ChannelSelection>
                </sld:RasterSymbolizer>
            </sld:Rule>
        </sld:FeatureTypeStyle>
    </sld:UserStyle>
</sld:UserLayer>

If i publish the layer with default raster style, i could see some black speckles around the boundary as below

How to remove those black region after it gets published in geoserver?


texts:
A have a raster image in .tif format (WGS 1984). I want to retrieve the cell value when input the pixel row and column. I found the most suitable tool
https://pro.arcgis.com/en/pro-app/tool-reference/data-management/get-cell-value.htm

I tried the following  code

result = arcpy.GetCellValue_management ("G:/06_DROUGHT_RISK/VCI/GeoTif/001_2001.tif", "309 320")

Even though it seems quite straightforward, every time i try to return cell value, I get NoData as the output

I also looked into the following link
https://community.esri.com/thread/69009

However I could not find the solution from given answers as well


texts: Google Earth Engine: How to automate calculating area for each class for each image in an imageCollection?

texts:
My goal is to be able to use the JRC Yearly Water classification History imageCollection to calculate the area of seasonal and permanent water for each year available in a specified extent. I am able to do this for a specific image but I am hoping to find a way to automate the process for all the images in the collection. I know that I can map over an imageCollection to analyze it, but I can't print a chart within that function or find a way for the function to return the values of the reduceRegion. The questions I have are:

1) Is there a way to map over the imageCollection and get the values of the area for each class for each image? Or should I be doing something else?

or

2) Is there a way to get a graph with the area of each class for each of the images? Or do a time series graph by class?

Here is my code for a single image. From here I would just export the csv and then combine them in a different program, but can I do it all in Google Earth Engine?

//import needed images and collections
var img1984 = ee.Image("JRC/GSW1_1/YearlyHistory/1984");
var gsw = ee.Image("JRC/GSW1_1/GlobalSurfaceWater");
var yearSeasonality = ee.ImageCollection("JRC/GSW1_1/YearlyHistory");

var roi = ee.Geometry.Polygon([[-113.31109206384143,38.70660386074948], [-113.09136550134143,38.69588646672195], [-112.93481032556018,39.112673412644575], [-113.20672194665393,39.12332806661866], [-113.31109206384143,38.70660386074948]]);

Map.addLayer(roi);

//These are the parameters if you want to map the extent or seasonality layers
var extentParams = {min: 0, max: 1};
var seasonalityParams = {min:1, max: 12, palette: ["#99D9EA", "#0000AA" ]};

//use the max extent of the water body to mask the seasonality layer
var mask = gsw.select('max_extent') .eq(1) .clip(roi);
var y1984 = img1984.updateMask(mask);

//calculate the area of the pixels
var y1984Area= ee.Image.pixelArea().addBands(y1984);

var classNames = ee.List(["No Data", "Not Water", "Seasonal Water", "Permanent Water"]);

//create a chart that adds up the areas of each of the classes

var TimeChart = ui.Chart.image.byClass({ image: y1984Area, classBand: 1, region: roi, reducer: ee.Reducer.sum(), scale:30, classLabels: classNames, }); print(TimeChart);

//another way to get the total areas for each class

var reduced1984 = y1984Area.reduceRegion({ reducer: ee.Reducer.sum().group({ groupField: 1, groupName: 'waterClass', }), geometry: roi, scale:30, bestEffort: true, });

print("reduced values",reduced1984);

Here is the code for me starting to automate the process (I'm not including the many things I tried that didn't work):

var seasonality = yearSeasonality.map (function(image){ var mask = gsw.select('max_extent') .eq(1) .clip(roi);

return image.updateMask(mask); });

var seasonalityArea = seasonality.map (function(image){ var seasons = ee.Image.pixelArea(); return ee.Image.pixelArea().addBands(image).copyProperties(image,['system:time_start','system:time_end']); }); print (seasonalityArea,"seasonality with area added,not summed");


texts: How to create multi criteria cost surface - Slope reclassification?

texts:
I´m working on thesis aimed at investigating of possible routes and paths between Mayans cities in Guatemala. My thesis should help archaeologist to locate and reveal old paths in jungle and speed up all works in field with some scientific background.

I will use 30m DEM and anisotropic LCP approach with Tobler´s hiking equation. Problem:

I need to create multi criteria cost surface like here: http://resources.arcgis.com/en/help/main/10.2/index.html#//009z0000001z000000

I don´t know how to create cost raster with correct parameters. I mean I don´t know how to reclassify slope and add costs to it. I tried to use slope

how it was created (just derived from DEM using "Slope" function
reclassify slope from 1. and, of course, results are different.

I can´t find correct range for intervals and its costs. I need the same, exact range like 0-5° (cost 1), 5-8° (cost 2),... based on scientific data.

And is this correct way of solving LCP problem?

Should I create multi criteria cost surface and use it as "in_cost_raster" and Tobler equation as vertical factor in Path Distance Tool or use slope (without reclassification) as "in_cost_raster", calculate Path Distance and then add cost of Land use?


texts: QGIS - points at the North and South Pole disappear in the composer

texts:
I am using QGIS 3.4 to produce maps of the Arctic and Antarctic with weighted points from a csv file.
I noticed that I lose the points at the North and South Pole when I add the map in the composer. They are present on the map, they just disappear in the composer (see pictures below).
Is there a way to avoid this?

The map in QGIS:

The map in QGIS's composer:

(in case you were wondering: I change the projection when making maps of the Antarctic. Projections used are Lambert Azimuthal Equal Area for the North and South Pole. I tried playing with other pole-centered projections and the problem remained)


texts: How to to get current legend column count using arcpy?

texts:
I can change the numbers of columns in the legend using arcpy.adjustColumnCount(column_count).

However, in my particular case I would like to know how many colums there are in the legend in the beginning and save that information as a variable, then change the number of columns to a specific value for a map and last change the number of colums back to the value of my variable.

Is there a possibility to do this?


texts: Creating geodjango model with geometry field that can accept different projections

texts:
I am working on web application with GeoDjango and within the model (db table) RoadAxis I have geometry field or more precisely LineStringField where I need to store different polylines (road axis) for different projects (each polyline has FK to project ID).
The problem is that projects will be done in different coordinate systems (i can't know them in advance) so the polylines that I will have to store will have different coordinate systems too.
I know that in PostGIS database I can have geometry column with rows in different projections but when I define my Geodjango model I need to define unique SRID for all them, or he takes 4326 by default.

There are few options that I considered but none of them really solve my problem.
I will try to go through some of them and I hope some of you could tell me if I got something wrong or if there is some much better idea

1) transform data just before saving to one single coordinate system (e.g. WGS84)

Problem with this approach:
I can't just tansform my geometries before saving them to database by using geodjango built in functions because I would lose survey accuracy that I need for some later calculations.

2) Dynamically create new table in runtime
Idea: First having only one table where geometries would be stored and it would be all in one projection system (e.g. wgs84) and then when the new project is added - if it has different coordinate system I dynamically create new table with geom field for that projection system and then I store data for that project in that new table

Problem: dynamically creating tables in django doesn't seem like happiest and easiest solution, and it could later cause lot of additional problems.

3) Using noSQL db so instead of dynamically creating new tables in relational database I could use something like MongoDB and have it document oriented
Problem: this add more complexity to application and probably it's not really necessary, beside I lose possibilities I have with PostGIS (spatial queries, eventual transformation functions,...)

4) creating geometry field with some fake epsg and put all geometries in it even though they have different projection and then create additional field where i would save real srid of geometry ..
Problem: not really correct and clean solution but it could help this work

5) ?

Maybe I am just missing something obvious but I really don't have clear idea what is the best way to do this.


texts:
Can anyone tell me where I can find a file like this on Pakistan?
I found it on ochaonline.un.org I cannot find it when I search on there.
http://www.scribd.com/doc/38714310/Afghanistan-Settlements-Index-A3


texts: How to set image size with Image Map Plugin on QGIS 2.4?

texts:
I use QGIS 2.4, but I have noticed that the "Set Map/Image Size" feature on the HTML Image Map Plugin no longer seems to work. No matter what you set the "Set Map/Image Size" for, your image will always be 928 x 725. This feature worked fine with QGIS 2.2, but I don't want to go back to 2.2, for a variety of reasons.

How do I get this feature to work with 2.4?


texts: Preserve datatypes when using SQLite and QGIS

texts:
Following Nathan's SpatiaLite/DBManager solution in this answer, I find that once the layer is loaded from DB Manager into QGIS (2.6.1), any calculated fields are appearing as text values instead of numeric types.  This is especially a problem because I would like to use those fields as the basis of style using graduated colors.

In the original table, the "pop" column is an integer datatype.  But when I look at the Layer Properties > Fields tab in QGIS, it clearly shows the field having type=QString, "Type name"=TEXT for these fields that were computed in the query, such as:

COALESCE(SUM(a2.pop),0) AS pop_neighbours

This seems to happen anytime there is a calculation in the SELECT statement, even for something as simple as:

a.pop + 1 AS pop1

I get the same result when trying to force the datatype:

CAST (a1.pop + 1 AS INTEGER) AS pop1

The datatype of any unmodified columns (such as a.pop) seems to be preserved without any problem.  The QGIS Layer Properties > Fields tab shows type=qlonglong, "Type name"=INTEGER.

Is there a way to preserve the datatypes of calculated fields?

UPDATE:

For what it's worth, I see the same issue sometimes when using QSpatiaLite, depending on the method.  These work fine in QSpatiaLite:

Create Table & Load in QGIS
Create Spatial Table & Load in QGIS
Create Spatial View & Load in QGIS

But these don't:

Load in QGIS as Spatial Layer
Create View & Load in QGIS
"CREATE TABLE AS" within a query


texts: Is there a way to get gdal and arcpy to work in the same python script?

texts:
I have done a little searching into this topic, but I really don't understand virtual environments, and setting up paths, using visual studio or anaconda all that well, so the when I read any advice I get overwhelmed and confused by it.

But basically I have ArcGIS Pro but also want to run some gdal specific tools in the same script. I was able to create a virtual environment with gdal in it using anaconda and I know how to get the the ArcPro's virtual environment/python but I can't download one library into the other's environment.

usually when I'm running scripts I am using pycharm and set the python interpreter to one of the two environment's and switch between the two depending on what script I'm using, but this isn't efficient or really working anymore.

Any ideas? I feel like I need an "Idiot's guide to installing gdal and arcpy" or something. I hope I made sense.


texts: Running into endless loop when adding feature to layer with pyqgis

texts:
I am working on a script where I want the user to add feature to a layer in edit mode, and then decide whether to keep the newly added feature in the layer or not. I am doing it the following way:

layer.featureAdded.connect(store_feature_id)
def store_feature_id(feature_id):
    featureIter =
    layer.getFeatures(QgsFeatureRequest().setFilterFid(feature_id))
    for feat in featureIter:
       feature = feat
    reply = QMessageBox.question(iface.mainWindow(),
                                    'Save Changes?',
                                    "Save your changes?",
                                    QMessageBox.Yes, QMessageBox.No)
    if reply == QMessageBox.Yes:
        layer.updateFeature(feature)
        layer.commitChanges()

The problem I have is that the command layer.commitChanges will trigger the layer.featureAdded condition again, thus this program will run in an endless loop. I am not sure why, because I expected that the feature add event to only emit once when the feature is initially added to the layer.

This problem can probably be solved by using an extra bool variable as flag, but I am wondering if there is a more elegant solution using the available pyqgis functions?


texts: Value Table in Python Script not displayed correctly?

texts:
I am trying to implement a Python Script tool where I need a table to input shapefiles and corresponding names for new created fields.

I found out that there is the Value Table provided by ArcGIS, but using the "Geostatistical Value Table" as data type in the properties of the script, it does not appear correctly in the user interface, see attachment where I just made a small test.

The tool is a python script tool (*.tbx) and the geostatistical analyst extension is enabled. My code so far includes only the Import Statements and Statements to read out the Parameter.

import arcpy
inputShp = arcpy.GetParameterAsText(0)

Then I passed on the the script to a custom toolbox and added the Parameter where I chose Geostatistical Value Table as data type.

After finishing the wizard, I can open the tool, but as you can see from the figure, it seems that there is too few space below test to enter a shapefile.

It should look somehow like in this example, right: Creating a ValueTable with a column that has dependency


texts: Calculating length and area of habitat polygons within overlapping buffers in ArcGIS

texts:
I have a created 500 m buffers around a series of sampling points and want to know how much of a certain habitat is contained within each buffer. I used intersect and summary statistics to calculate length and area of habitat a within one example buffer x and got a value y.

When using  the entire set of buffers, the result y for habitat a in buffer x is larger than before. Buffer x overlaps with several other buffers. I am interested in the total amount of habitat within each buffer, regardless of the overlap with another buffer. Hence the result should not be different.

It seems that when several buffers overlap the intersect result consists of multiple sections so that the total length calculated using summary statistics is different.

How can I avoid this problem and get a correct output ?


texts: Changing from data view to layout view using ArcPy?

texts:
Changing from layout view to data view using ArcPy?, has been asked and answered, but I can't figure out how to switch to the layout view from data view using ArcPy and ArcMap 10.

Can it be done?


texts: Changing the coordinate value of a Pointz object with python

texts:
How can I change the X,Y or Z value of a Pointz feature with python in QGIS?
I am trying to change the Z values of Pointz with a python plugin.   I am first trying to pass the selected feature to a variable in the program.  I can get to the geometry in the program and on the console.

selection = layer.selectedFeatures()[0]
geom = selection.geometry()
print(geom)
<QgsGeometry: PointZ (26141648.89000000059604645 642064.42000000004190952
767.54999999999995453)>

Then I get stuck trying to get the Z coordinates out.   I have tried the geom.z() and geom.id() and I get errors. geom.asPoint()[1] does give me the X and Y values but not the Z.   Is there a way to get and change the "Z"

<QgsGeometry: PointZ (26141648.89000000059604645 642064.42000000004190952
767.54999999999995453)>
geom.x()
Traceback (most recent call last):
File "C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\code.py", line 90, in runcode
exec(code, self.locals)
File "<input>", line 1, in <module>
AttributeError: 'QgsGeometry' object has no attribute 'x'

UPDATE

This is similar to what I ended up with.   The lines with wkt in them generate the geometry from "well known text".   This did work but I took it out for the better .setX() method.  You would have to change the feature.setGeometry(geom) to feature.setGeometry(new_pnt) to use that method.

layer = iface.activeLayer()
    print(layer)
    all_features = layer.getFeatures()
    layer.startEditing()
    for feature in all_features:
        print(feature['NUMBER'])
        if feature['NUMBER'] == '1003':
            print('found it ********\n  ')
            feature['SUMPDIST'] = '18.55'
            geom = feature.geometry()
            print(geom)
            new_feat = QgsFeature()
            point1 = QgsPoint(26132109.5, 644411.8, 940.600)

            #new_pnt = QgsGeometry.fromWkt("POINTZ(26132109.5 644411.8 940.48)")
            wkt = "POINTZ(" + str(point1.x()) + " " + str(point1.y()) + " " + str(point1.z()) +")"
            new_pnt = QgsGeometry.fromWkt(wkt) #This well known text works also But I like the .setX, .setY, .setZ method much better.
            geom.get().setX(point1.x())
            geom.get().setY(point1.y())
            geom.get().setZ(point1.z())
            feature.setGeometry(geom)

            layer.updateFeature(feature)
            print('feature '+feature['SUMPDIST'])
            print(str(point1.z()))
            print(str(feature.geometry().get().z()))

            break
layer.commitChanges()
print(str(feature.geometry().get().z()))


texts: Creating Shapefile of subset of input features in PyQGIS?

texts:
I'm having some trouble applying my Python skills to develop some scripts in PyQGIS. For me it's far simpler to write it in ArcPy, but then again I wouldn't get that much satisfaction after finishing my project. Trying to work with open source tools as much as possible.

In short, for now I want to write a Python script that takes in a shapefile and from that generates another one with some specific attribute data related to buildings. If I were doing it in ArcPy, my algorithm would be the following: use the MakeFeature Layer_management, then perform a SelectLayerByLocation_management, within the selection select the data you need, and then create another layer with desired data.

How could I this in PyQGIS?

One idea, reading the PyQGIS Cookbook, to use .getFeatures(), QgsVectorFileWriter to write the shapefile, but how do I perform the query?


texts:
I Have code

   ISelectionSet pSelectionSet;
   layer = AxMapControl1.get_Layer(0);
   pFeatSel = (IFeatureSelection)layer
   pSelectionSet = pFeatSel.SelectionSet;
   ITable table = pSelectionSet  as ITable;

How to convert selection to a table?


texts: Open/Free Area Code shapefiles for the US/Canada

texts:
I would like to map magicJack's areacodes and exchanges against known area codes and exchanges to create a coverage map, just to see what it would look like as compared with, say, a mobile phone coverage map.

Having read this previous post:
Shapefiles for U.S. Area Codes (NANP)
I see the only available databases seem to cost more than I'm willing to pay. And the maxmind information is point data rather than polygon data. I'm curious if any new free/open resources have become available since the previous posting.


texts: Deciding whether to cast or Project when meters needed for ST_Distance?

texts:
I have a database and I am using geoms to build polygons around points.

I need to do st_distance, st_area, st_overlaps, st_union, and st_instersects. also using st_x and st_y and st_project, so pretty sure i have to use geometry data type for these operations.

My issue is that I need meters, not degrees to get st_distance between points, and and st_area of polygons.

Should i do a cast with the :: operator to convert to geography or should i project?

Also, how do I choose which srid to pick?

My input could be anywhere in the US (including AL, HI).


texts: Tool to find the nearest points in one layer, ArcGIS

texts: How can I find nearest points in the same layer in ArcMap? I would like to find points that are in the distance of less then 0,5m.


texts: QGIS Layer 1 Point. Layer 2 Line. Calculate distance along the line that the Point sits

texts:
As per the title really.
I have a layer consisting of a network of Lines. I have another layer of point items 'snapped' to the lines. I'd like to attach a length value to the point based on how far along the line it is.


texts: Automatically moving along columns within a layer as moved to the next item in a QGIS Atlas?

texts:
My problem is quite hard to explain, I have attached some images to help illustrate my issue.

I want to produce a large number of maps which overlap in geographic area. Currently my data is arranged as shown in image 1.

I have a number of sites in each column, which colour UK wards based on differing values.

I want to set up a system in QGIS which will map Site 1, and then move along the colons to map site 2 and so on.

Unfortunately each ward is affected by multiple sites. Ward E050001 may have a value for site 1 of 0.75, and then a value of 0.25 for site 2 etc.

I currently manually change the values to move along the columns in the layer properties as shown in image 2, and it's this manual step that prevents me from using a QGIS Atlas.

If Image 3 was my map, I'd like it to display the area around site 1, export a map, and then move on to site 2 and export a map. This is a simplification however, as the areas often overlap for sites, and I'd like the colour to vary based on the value in the cells shown in image one.

I feel the best way would be to use an Atlas and a bit of python to move each column along as it moves to the next item in the Atlas, but an unsure how to do this or if it would even work.

Does anybody have any idea how I would go about achieving this?

It's trivial to do manually, but I have potentially thousands of maps that I need to do!


texts: Pansharpened Raster Color Changes After Clipping

texts:
I have a mosaic raster of an agricultural area. In the mosaic raster, a specific agricultural field (field no. 5) looks just like this

When I clip the raster with the field's shapefile (shapefile of the field no. 5) the color changes a bit and looks like the following image.

The following codes have been used to clip the mosaic raster -

#clipping the raster
arcpy.Clip_management(in_raster= mosaic_pan_layer, rectangle=str(extent_of_shape_file),
out_raster=pan_save, in_template_dataset=shape_file, nodata_value="65536",
clipping_geometry="ClippingGeometry", maintain_clipping_extent="NO_MAINTAIN_EXTENT")

#saving the clipped raster
arcpy.CopyRaster_management (in_raster = pan_save, out_rasterdataset = Clipped_pan_save,
nodata_value = "9999")

Is there any way to get the clipped raster to look exactly the same as the mosaic layer?


texts: Calculating nearest distance between 2 points and creating line between them avoiding barrier using ArcGIS Desktop?

texts:
I have 2 point datasets.
One shows a number of offshore locations.
The other is the location of a number of ports in the UK.

I want to get the name and distance to the nearest port for each offshore location.  I am not allowed to cross over land (I have a polygon of the UK).

I have access to ArcMap 10.2.2 advanced and AGOL.

I am very experienced with ModelBuilder but have little to no experience with Python.

I have looked at the Near tool, but this does not allow for barriers.

The only method I can think of is:

Create a raster that incorporates the area that all my points fall within
Assign a value of 1 to the land and a value of 0 to the sea
Run some sort of cost analysis

Using this method, I can't think of a way that would let me generate lines showing the "route" between my offshore points and the port locations.

I imagine a Python solution would be best, but I am open to suggestions.