texts:
Hello I'm abit new here and I hope I'll even set the question right.

Recently I've started to work on a project that will include displaying maps in qgis-web-client. At the moment I've set up the server to display the basic examples. But now I'm trying to import my own qgis project. This specific project has an custom CRS. So far I've discovered on this site that the webclient does support custom CRSs. So without hesitating I've followed the guide and did the steps. The custom CRS file I created was Proj4js.defs["EPSG:4800"] = "+title=nevem +proj=tmerc +lat_0=0.0 +lon_0=15.0 +k=0.9999 +x_0=500000 +y_0=-5000000 +units=m +no_defs";

This CRS was copied from the project in the qgis desktop application. And I edited it so it would match with the other ones. After I run the web-client I get an error from the browser console that this.resolutions is undefined. Also I should mention that the layer names load properly. For the sake of testing I tried using the default CRS just to test my qgis project with it and it loaded the layers, but they were projected incorrectly because of the CRS.

So now the question is: How to make this CRS work with qgis-web-client? And what do the a and b variables mean in:Proj4js.defs["EPSG:3857"]= "+title=WGS 84 / Pseudo Mercator +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs";


texts:
I have layers based on GeoTIFF files, provided by GeoServer, that are displayed as overlays over base maps from OpenStreetMaps, Bing, Google via OpenLayers. I want to mask out areas from the GeoTIFF overlays but enable the base maps to continue to their full extent / appear around the boundaries of my mask.

How to crop coverage with WPS chain in
sld? had the clues for creating a SLD that includes a WPS transformation that crops a layer to a shape definition. That works fine for the layer by itself, i.e., in GeoServer's own OpenLayers layer preview. But if I display that layer on top of a base map in my application, performance of the base map layer display plummets (up to minutes to display vs nearly instantaneous), CPU usage on my laptop spikes, and I end up with pink tiles where some of the base map should be. I suspect the performance hit occurs because of the (relatively simple, 10 point but continent sized, polygon used by the cropShape definition.

This GeoTools CropCoverage source code suggests that I could use a coverage definition provided by an input raster.

My questions: Can that raster be defined as part of GeoServer WPS SLD transformation? How do I define that input raster? Is there an example somewhere? Ultimately, will that improve base map display performance or is my suspicion wrong and there is some other reason for display being so slow when there is a base map in addition to the overlay with this SLD change is made for the overlay?


texts: Euclidean Distance functioning as Euclidean Direction in Arcpy

texts:
I'm trying to call Euclidean Distance in an arcpy script within a toolbox. I've used arcpy.sa.EucDistance but the output raster is a Euclidean Direction raster. When I run it in the python dialog, it works as intended and when I run the tool Euclidean Distance tool it works normally.

Does anyone have an idea of what's going on? relevant code below.

import arcpy
from arcpy.sa import *

lines = r"C:\shapes\Lines.shp
output = r"C:\rasters\EucDistLines.tif"

EucDist = arcpy.sa.EucDistance(lines, 6000, 3, output)


texts: Removing z-dimension while writing shapefile using PyQGIS

texts:
There is such a code:

from PyQt4.QtCore import QVariant
from qgis.core import (QgsFeature, QgsField, QgsFields,
QgsGeometry, QgsPoint, QgsVectorFileWriter)
from qgis.utils import QGis, iface
iface.legendInterface().layers()
sluoksniu_saras=iface.legendInterface().layers()
for cL in sluoksniu_saras:
    fn='D:\\Ko\\'+cL.name()
    ret = QgsVectorFileWriter.writeAsVectorFormat(cL, fn, "Windows-1257",\
    QgsCoordinateReferenceSystem(cL.crs()), driverName="ESRI Shapefile",\
    forceMulti =False, includeZ = False)
    if ret == QgsVectorFileWriter.NoError:
       print "Irasyta i " +fn + " !"
print 'viskas'

In theory, he should write all open layers of the project in a specific folder and reset the z-dimension. Layers are recorded, but the z-dimension is not reset. Where is the mistake? I'm using QGIS 2.16


texts: Compare list of field names located in two separate dictionaries with same key

texts:
I created two dictionaries, dictionaryone and dictionarytwo, and both contain the same key. The key is a file path to a shapefile, for example: C:\TestFolder\Folder1\A\address_pt.shp. In addition, dictionaryone has a value that is a list of field names contained within the attribute table of the shapefile.

Below is an example of dictionaryone output:

{'C:\\TestFolder\\Folder1\\A\\address_pt.shp': ['FID', 'Shape', 'FULL_ADDRE', 'STREET_ADD', 'CITY_ST_ZI', 'PREFIX_DIR', 'HALF_ADDRE', 'HOUSE_NUMB', 'STREET_NAM', 'STREET_TYP', 'SUFFIX_DIR', 'UNIT_TYPE', 'UNIT_ID', 'CITY', 'STATE', 'ZIPCODE', 'PARCEL_NUM', 'PARCEL_N_1', 'HOUSE_NU_1', 'HOUSE_NU_2', 'STATUS', 'Assembly_S', 'Senate', 'District', 'Council', 'GRIDNUM', 'PHYSICAL_A', 'PUBDATE']}

Furthermore, dictionarytwo contains a list of field names that are thought to be within the shapefile. Below is an example of dictionarytwo. Notice that both dictionaries have the same key.

{'C:\\TestFolder\\Folder1\\A\\address_pt.shp': ['house_numb', 'half_addre', 'prefix_dir', 'street_nam', 'street_typ', 'suffix_dir', 'unit_type', 'unit_id', 'street_add', 'city', 'state', 'zipcode']}

My goal is to use the list of field names in dictionarytwo and search through dictionaryone to see if field names from dictionary two exist in dictionaryone. I assume that I have to standardize the strings in both dictionary lists in order to search and if that assumption is true, I would like to capitalize the lists in both dictionaries before searching.

In addition, I would like to capture the results of the search and input the results into two new empty fields. The first field would be named something like field_found and include the matching fields found in both dictionaries. The second field would be named something like fields_not_found and include the fields from dictionarytwo that are not found in dictionaryone.

I'm not sure if this process is possible and I can't figure out the code logic. I'm not sure if using dictionaries is the best approach. Below is the code that I currently have and it creates both dictionaries.

import arcpy

arcpy.env.workspace = r"C:\TestFolder"
dbf = "ProcessControlTable"
dbf1 = "APExistOnlyEdit.dbf"
myfield = "FI_PATH"

filepathlist2 =[row[0] for row in arcpy.da.SearchCursor(dbf1,myfield)]
#print(filepathlist2)

#Dictionary holding Field Names that are within Shapefile
dictionaryone = {}

for thisFile in filepathlist2:
    ContainedFieldNames = [f.name for f in arcpy.ListFields(thisFile)]
    dictionaryone[thisFile] = ContainedFieldNames

print(dictionaryone)

##Not sure how to insert field list into dictionary two for searching
##FieldList =('FI_PATH','AP_FIPS','AP_BLDGCOM','AP_TYPE','AP_STATUS','AP_ID','AP_HOUSENU','AP_HALFADD','AP_PREDIR',
##         'AP_PRETYPE','AP_STNAME','AP_SUFTYPE','AP_SUFDIR','AP_UNITTYP','AP_UNIT','AP_BUILDIN',
##         'AP_FULLADD','AP_CITY','AP_STATE','AP_ZIP')

#Dictionary holding Field Names that are thought to exist in Shapefile
dictionarytwo = {}

with arcpy.da.SearchCursor(dbf1,[myfield,'AP_FIPS','AP_BLDGCOM','AP_TYPE','AP_STATUS','AP_ID','AP_HOUSENU','AP_HALFADD','AP_PREDIR','AP_PRETYPE','AP_STNAME','AP_SUFTYPE','AP_SUFDIR','AP_UNITTYP','AP_UNIT','AP_BUILDIN','AP_FULLADD','AP_CITY','AP_STATE','AP_ZIP']) as rows:
    for row in rows:
        dictionarytwo[row[0]] = [a for a in row[1:] if a.strip() != '']

print(dictionarytwo)


texts: Accessing individual pages within existing pdf using arcpy.mapping?

texts:
I've used the arcpy.mapping module for some time now to create mapbooks and other pdf documents, using a combination of Data Driven Pages and Python. I would like to take an existing pdf and save each individual page as separate pdfs.

I have read through the online documentation for PDFDocument, and while I can see the usefulness of the pageCount property, I don't see how I can access an individual page. The only reference I see to a page number is when using the insertPages method to insert a page before an existing page.

I know that using something like pypdf enables one to iterate over the pages in a pdf and other useful things, but I'd like to know if I can do this using arcpy before downloading it.


texts:
I have a feature with hundreds of lines and I need to extend all the lines of a given distance (for example all of 100 meter).

Is there a function in ArcGIS to do it automatically?


texts: Problem in the operation "Add xy Coordinate" in ArcGIS10

texts:
I have a problem in the operation "Add xy Coordinate" in ArcGIS10.
For a point shp theme with  many point objects I did “ArcToolboox - Data management- Features- Add xy Coordinate”,  after that the program informed that operation was created correctly, but the meanings of all point objects in the created fields POINT_X and POINT_Y – are zero (0).

What can be a reason of this error? Maybe something related with Workspace?


texts: QGIS Not Importing ExtendedData from Self-Exported KML?

texts:
I've been scouring the internet for the past few hours, and I just can't seem to figure it out!

I am using QGIS 2.4.0-Chugiak.

Here are my simple steps:

1) I create a new ShapeFile (Type: Polygon) and add 3 attributes:

2) I draw a simple polygon with 6 coordinate points and give it the following attributes:

3) I now export the Shapefile as a KML and it gives me the following text code:

    <?xml version="1.0" encoding="utf-8" ?>
<kml xmlns="http://www.opengis.net/kml/2.2">
<Document id="root_doc">
<Schema name="test1" id="test1">
    <SimpleField name="cascade_na" type="string"></SimpleField>
    <SimpleField name="sectNumb" type="int"></SimpleField>
    <SimpleField name="Market" type="string"></SimpleField>
</Schema>
<Folder><name>test1</name>
  <Placemark>
    <Style><LineStyle><color>ff0000ff</color></LineStyle><PolyStyle><fill>0</fill></PolyStyle></Style>
    <ExtendedData><SchemaData schemaUrl="#test1">
        <SimpleData name="cascade_na">test1</SimpleData>
        <SimpleData name="sectNumb">9</SimpleData>
        <SimpleData name="Market">Miami</SimpleData>
    </SchemaData></ExtendedData>
      <Polygon><altitudeMode>relativeToGround</altitudeMode><outerBoundaryIs><LinearRing><altitudeMode>relativeToGround</altitudeMode><coordinates>-73.977273863779757,40.764000085624744 -73.976720692372155,40.764083880580976 -73.976333472386827,40.763925601130516 -73.976480984762205,40.763781287185402 -73.976837473002647,40.763613696404278 -73.977390644410249,40.763674215346221 -73.977273863779757,40.764000085624744</coordinates></LinearRing></outerBoundaryIs></Polygon>
  </Placemark>
</Folder>
</Document></kml>

4) Finally I go ahead and re-import what I just exported, using "Add Vector Layer" hoping to get all the attributes back, but I get none of them!

I just want to get all the attributes imported.

What am I doing wrong?


texts:
I have some data that I collected from Trimble Juno 3D. When I download the data, I could not see them in the map. The coordinates are NOT in the correct place where it should be. Their coordinates are like this:

Long. -2.652497e+014;
 Lat. 3.569410e+014.

Is it possible change it to decimal degree?


texts:
I have got a .json file to display in Leaflet map.The points are shown with custom blue dot markers. I cannot get to change them.How do I get to change the symbol?

var greenIcon = L.icon({
  iconUrl: 'marker-shadow.png',
  iconSize: [38, 95] // size of the icon
});

var jsonTest = new L.GeoJSON.AJAX(["colleges.geojson"],{onEachFeature:popUp},{icon:greenIcon});
var jsonTest1 = jsonTest.addTo(m);


texts: Extracting distances along a line from postcode a centroid to the nearest point (bus stop)

texts:
I'm looking to measure the distance from postcode centroids (point file) in a city to the nearest bus stop (point file) along the road network (line file).

In QGIS I can calculate the distance but it is not measured along the road network. I'm sure this can be performed in GRASS GIS however, but unsure to how (v.distance maybe). Its the from and to part of the v.distance code I'm struggling with, as I can't find where to specify along the network, or does all three files need to be connected to each other and in different layers?.


texts: Multiuser Editing with ArcGIS through LAN

texts:
I am working on an ArcObjects application (C#), which is supposed to carry some customized editing operations, the database is ArcSDE based on Postgresql via (accessible via a LAN network).

One feature we're keen on is multi-user-edition (MUE) : multiple users should be able to edit the same layer (or objects in a layer) simultaneously.

The problem is that the application is frequently (once every week) overwhelmed (congested) thus disabling the MUE and leaving only one user to carry all the edition workload.

I am not very sure this is specifically related to the application (though there is a fair likelihood it is), so I had to wonder whether such an issue is experienced by users accessing a DB (via ArcSDE) and using ArcGIS desktop as an editor ; * in other words, my hope is to get an idea about how solide is ArcGIS when it comes to MUE (how well does it deal with conflicts, scalability and data size when MUEing)*. (Heuristically, if it turns out that ArcGIS is infallible or near-infallible in this regard, I'll then know that I have to optimize somehow the application)

Thanks in advance for your feedback.


texts: Removing spatially duplicate features using ModelBuilder?

texts:
I am creating a model to remove features which are duplicated within a polygon feature class. Most of the features are duplicated once, but about a quarter have up to three or four exact copies. I do not have ArcInfo, so I cannot use Delete Identical.

My model should contain the following:

Make Feature Layer on the input feature class specified by the user to output Output_FL
Get Count to count the number of rows within the feature class
Create a variable counter and initialise to 0
Add a While iterator where counter < row count
Select Layer by Attribute on Output_FL to create a NEW_SELECTION where OBJECTID = counter
Select Layer by Location on Output_FL where the relationship is ARE_IDENTICAL_TO
Select Layer by Attribute on Output_FL to REMOVE_FROM_SELECTION where OBJECTID = counter
Delete Features from Output_FL
Increment counter

I have finished all the steps except step 9. I am not sure where/how I should increment the counter. I'm also not sure if I should put a Copy Features tool somewhere before the while loop as a backup mechanism since this model will be modifying the input data, and not everyone will make a backup of their data before running the model.


texts: Export/Import CRS for use on similar/same file?

texts:
I have a slightly strange request. I have access to a set of geo-ref'd PDFs that have multiple layers within the file. As touched on in other threads, there's no good way (at this time) to turn on/off layers from a PDF in QGIS. So I end up turning off/on the layers I want (in Acrobat) and then printing a new PDF file of the same size/resolution (however, the geo-reference metadata is lost when I print the new PDF).

So what I'd like to do, is import the geo-ref'd PDF into QGIS, export the CRS data, and then import the non-ref'd file and apply the CRS from the original file.

Can this be done? If so, how?  And if you know how to control PDF layers within QGIS, let please let me know.


texts: How to undelete from ArcGIS versioned table?

texts:
We have a versioned table of polygons and one of our users have deleted every geometries.

With Oracle's SQL Developer, i can see the table still containing the missing data, along with its Axxxx and Dxxxx tables, representing the versioned actions (the Axxxx table stores the data added to our table and the Dxxxx the row number of the deleted data).

I understand that ArcGIS Desktop use an Oracle View that computes the added and deleted data to the table.

Is there a way to restore or undelete those polygons (and their attributes) ? Can I just empty that Dxxxx table, which contains row numbers of the deleted polygons ?

We have ArcGIS Server 10.2.2, ArcGIS Desktop 10.4 and Oracle 11.


texts: Shp2pgsql - exec() not allowed on external server

texts:
Me and my collegue are using the shp2pgsql and psql commands in php code (a Drupal module to be precise).

We put these commands in a php variable and use the exec() command to execute them. This works very well but now we have to use these operations on an external server which doesn't support the exec() command.

Is there any workaround or solution for this matter?

Thanks in advance.


texts:
I have a list of x,y,z points named simply "A".

I'm trying to convert them to a multipoint with Shapely.

from shapely.geometry import Point, MultiPoint

orig = Point(1.9,2.5,3.1)
destinations = MultiPoint([(A)])
print(destinations)


texts: Can you analyse .TAB files in QGIS? Or compare .shp to .tab?

texts:
I'm recently experimenting with QGIS at work (a poor charity who is struggling along with MapInfo v8).

As all our layers are in .TAB format and I was wondering if you can use the geoprocessing tools to analyse this data in QGIS or if you have to convert to .SHP first?

Also, if I was being difficult, could I use the intersect tool to see how much of a .SHP polygon was inside a .TAB polygon?

I didn't work when I tried, but this may be due to me missing something else.


texts: Colliding labels for point features in QGIS

texts:
I have a layer (point features) containing various points of interest (POIs), some of them have the exact same location. I would like to map these POIs and add a label showing their names.
However, for those POIs that share the same location, QGIS just display the overlapping labels and it is not readable.

Is there a way in QGIS to automatically place labels for point features so that they never overlap, even though they belong to the same layer?

I found this relevant thread, but no solution to the problem.


texts: postgis and st_intersects ; fail to compare a table to itself with a left join

texts:
postgresql 9.6, postgis 2.3
Following this Select features that do NOT intersect in PostGIS, I would like to compare one table to itself, to keep all records from table p1, even and above all those that don't have any record matching in table p2.
Actually I want to keep features that don't intersect any.

edit : code updated

select p1.gid, p2.gid from table p1
left join table p2
on st_intersects(p1.geom, p2.geom)
and p1.gid <> p2.gid
where p2.gid is null ;

Without the and, All that I get are matching records : the left seems to be ignored.
With the and condition, no result.
Where is (are) my mistake(s) ?


texts: Calculating longest flow path and avg. slopes for catchments with qgis?

texts:
I am a hydrologist and there are 2 important pieces of information I need for modeling. I need to know the longest flow path for each sub catchment in my area of interest and I need to know the elevations of point specific distances along those paths.

As an example, I have a shapefile with specific boundaries that represent my catchment and an elevation raster. I would like to be able to determine the longest path that water flows in each of these defined areas. Once I have a shape file containing the longest flow paths for each catchment, I would like to find the elevation at 10% and 90% the length of each flow path and the distance between each point.

I have not found either of these functions in anything except the Arc Hydro tool bar for ArcGIS. I know that the individual tools for creating these things exist in QGIS ( through SAGA and GRASS ) but I am new to python and not a programmer in any respect. I am hoping these tools are in an existing package and I am just too new to qgis/grass/saga to find them. Any help would be greatly appreciated.

I have been using QGIS for about 18 months so I have some idea what I am doing with it.


texts: Shapely buffer(0) returns an empty polygon(coordinates:[]), and discard the main polygon


texts: Changing the time of a map label in QGIS using python

texts:
I would like to adjust the "[% format_date( now(),'dd/MM/yyyy HH:mm:ss') %]", as it currently displays 07/12/2020 13:31:40 on the map, as I am in the UK right now. However, I would like to change the time to fit time zone that is 6 hours behind me (a South American country). Is there any way I can do this?
map_label = QgsLayoutItemLabel(layout)
map_label.setText("[% format_date( now(),'dd/MM/yyyy HH:mm:ss') %]")
map_label.setFont(QFont("Arial", 18))
map_label.adjustSizeToText()
layout.addLayoutItem(map_label)
map_label.attemptMove(QgsLayoutPoint(40, 33, QgsUnitTypes.LayoutMillimeters))
map_label.attemptResize(QgsLayoutSize(143, 120, QgsUnitTypes.LayoutMillimeters))


texts: Saving rasters to postgresQL using raster2pgsql, saving forecast time

texts:
I would like to save a forecast which is saved as geoTiffs into a postgresQL database. I've been able to import the rasters using:

raster2pgsql -s 0 -I -d tiff/field*.tiff -F public.latest | psql -d optigrid

THis gives me a table latest which looks like

rid  | rast  | filename
1      ...     field1.tiff
2      ...     field10.tiff

and so forth.

However I would like to add another column which stores the time the forecast, something like

rid  | rast  | filename      | forecast_date
1      ...     field1.tiff      2012-08-09
2      ...     field10.tiff     2012-08-10

Is this possible?

Thanks Matt


texts: Is there a SDE command to export/import feature class schema only without data?

texts:
I'm using ArcGIS Server 10.0, ArcSDE 10.0 and Oracle 11g.

I'd like to create a sde export of my layers' schema. The export would be used to create easily those layer in a new database instance.

Do you know of an export/import SDE command that allows that?

Or at least, is it possible to import an XML schema file using a SDE command or Python script?


texts:
I'm a totally newbie in MapServer. I'm using MapServer 6.0.3. I have managed to label some lines but I cannot create a background colour for the labels.

Using QGIS 2.10 I got something near to what I'm trying to do with MapServer. I tried to export using RT MapServer Exporter but the map file that was created didn't do the job at all.

Now, this is the map file that I have so far:

MAP
  FONTSET "H:/mapserver/data/fonts/fonts.list"
  EXTENT 19 33 30 42
  IMAGETYPE "png24"
  NAME "airways"
  SIZE 1280 1024
  STATUS ON

  LAYER
    NAME "airways"
    DATA "H:/mapserver/data/airways.shp"
    TYPE LINE
    LABELITEM "NAME"

    CLASS
      NAME "airways"
      LABEL
        COLOR 0 0 0
        TYPE TRUETYPE
        FONT Vera
        SIZE 7
        ANTIALIAS TRUE
        POSITION CC
        ANGLE AUTO
      END # end of label

      STYLE
        COLOR 0 0 0
        WIDTH 0.25
      END # STYLE
    END # CLASS

  END # LAYER

END # MAP

The output of this map file is the following:

I have added in LABEL the following lines to create background colour for the labels, changing at the same time the color in LABEL to 255 255 255 (white font, black background)

STYLE
  GEOMTRANSFORM 'labelpoly'
  COLOR 0 0 0
END # STYLE

The STYLE block was added exactly below the ANGLE AUTO.

Only one label was drawn on the map, the last entry in my shapefile.
It has white color and black background but it's only one label.

Any guidance on how to keep all my labels with background color?


texts: Filtering Polygon Data with many statements in QGIS

texts:
Working with QGIS 3.4.1 - Madeira

I am trying to perform a simple Filter, but when I include more than 31 values of a given field I get Zero returns.

Is there a limit to the number of statements that can be used in the filter?
I am using the following expression:
"Field_name" = 'Value1' OR "Field_name" = 'Value2' OR... and so on


texts: Add cookies to a WMS request with OpenLayers2

texts:
I need to add a cookie in each WMS request.  I tried to do something like below, but I don't think this is the right way:

getURL : function(bounds) {
    var requestString = OpenLayers.Layer.WMS.prototype.getURL.apply(this, arguments);
    window.document.cookie="referrer="+window.location.origin;
    return this._encode(requestString) + "&KEY=" + this.navkey;

Is it possible to do this?


texts: Specific renaming of elements in Attribute table failing until closing and reopen attribute table

texts:
Well the title tells it already. I have a column name and 100 objects. I have given them all a name. Now I mentioned that I have written some wrong of them. So I tried to correct them by double click and rename. But when I leave the name it remains the old name. I need to close the Attribute table and reopen it, then I see that it has taken the name.

This happen just to some of my entries! Most are working properly and change is done without reopening.

I just want to know if it's a bug or a feature :D I mean I know know that he just do not show it to me. But it's a bug that can be easily fixed I think.


texts: Symbols changed in conversion from QGIS to Webmap

texts: When I convert a map in QGIS to a webmap, the symbols are changed--the size is ok, but all the specific symbols-- stars, diamonds, etc.--are converted to circles. Does anyone have suggestions on how I might address this?


texts:
I have 4 tiles of worldview images. in some part of each image they overlap each other. I would like to have a single output having them all together. As far as I know I should Mosaic them, right?
I have tried to mosaic them in Envi 5 and Qgis. Its strange what I get as an output.

ENVI--> Even though I have selected all the bands but still I get single grayscale band as an output.

QGIS--> I tried to mosaic the tiles but the output is also weird.


texts: QGIS won't add layers into a saved project (but will into an unsaved project)

texts:
I have some VRT layers. If I add them into a blank unsaved project they work fine. If I save a blank project and try to add any of the files I get an error.

Invalid Layer: GDAL provider Cannot open GDAL dataset /GIS DATA/OML VRT/Buidlings.vrt: `/GIS DATA/OML VRT/Buidlings.vrt' not recognized as a supported file format. Raster layer Provider is not valid (provider: gdal, URI: /GIS DATA/OML VRT/Buidlings.vrt

I suspect this is something to do with paths but I can't figure out what it might be.

Sample of a bit of one of the VRTs

<OGRVRTDataSource>
<OGRVRTUnionLayer name="Buildings">
<OGRVRTLayer name="SD_Building">
<SrcDataSource relativeToVRT="1" shared="1">Building</SrcDataSource>
<SrcLayer>SD_Building</SrcLayer>
<GeometryType>wkbPolygon</GeometryType>
<LayerSRS>PROJCS[&quot;British_National_Grid&quot;,GEOGCS[&quot;GCS_OSGB_1936&quot;,DATUM[&quot;OSGB_1936&quot;,SPHEROID[&quot;Airy_1830&quot;,6377563.396,299.3249646]],PRIMEM[&quot;Greenwich&quot;,0.0],UNIT[&quot;Degree&quot;,0.0174532925199433]],PROJECTION[&quot;Transverse_Mercator&quot;],PARAMETER[&quot;False_Easting&quot;,400000.0],PARAMETER[&quot;False_Northing&quot;,-100000.0],PARAMETER[&quot;Central_Meridian&quot;,-2.0],PARAMETER[&quot;Scale_Factor&quot;,0.9996012717],PARAMETER[&quot;Latitude_Of_Origin&quot;,49.0],UNIT[&quot;Meter&quot;,1.0],AUTHORITY[&quot;EPSG&quot;,&quot;27700&quot;]]</LayerSRS>
<Field name="ID" type="String" src="ID" width="38"/>
<Field name="FEATCODE" type="Real" src="FEATCODE" width="38"/>
</OGRVRTLayer>
<OGRVRTLayer name="SE_Building">
<SrcDataSource relativeToVRT="1" shared="1">Building</SrcDataSource>
<SrcLayer>SE_Building</SrcLayer>
<GeometryType>wkbPolygon</GeometryType>
<LayerSRS>PROJCS[&quot;British_National_Grid&quot;,GEOGCS[&quot;GCS_OSGB_1936&quot;,DATUM[&quot;OSGB_1936&quot;,SPHEROID[&quot;Airy_1830&quot;,6377563.396,299.3249646]],PRIMEM[&quot;Greenwich&quot;,0.0],UNIT[&quot;Degree&quot;,0.0174532925199433]],PROJECTION[&quot;Transverse_Mercator&quot;],PARAMETER[&quot;False_Easting&quot;,400000.0],PARAMETER[&quot;False_Northing&quot;,-100000.0],PARAMETER[&quot;Central_Meridian&quot;,-2.0],PARAMETER[&quot;Scale_Factor&quot;,0.9996012717],PARAMETER[&quot;Latitude_Of_Origin&quot;,49.0],UNIT[&quot;Meter&quot;,1.0],AUTHORITY[&quot;EPSG&quot;,&quot;27700&quot;]]</LayerSRS>
<Field name="ID" type="String" src="ID" width="38"/>
<Field name="FEATCODE" type="Real" src="FEATCODE" width="38"/>
</OGRVRTLayer>

I really like the way the VRTs work when they work and I know they are not fully supported but this behaviour seems odd.

QGIS 2.18 Mac OSX 10.11.6 KyngChaos


texts: Using Attribute Assistant's Autonumber or GenerateID tools in a versioned environment

texts:
Is there a workaround for using Esri's Attribute Assistant's Autonumber or GenerateID in a versioned environment without creating duplicates?

re: Attribute Assistant and Versioning

could Python autoIncrement perhaps accomplish new unique IDs in a versioned environment if AA cannot?


texts: How do I link an existing install of GDAL to FileGDB in Linux?

texts: The distribution I use (CentOS) only has GDAL 1.8 available through elgis. I've downloaded the FileGDB API from ESRI and I'd like to be able to link it after an install. How can I do this?


texts:
I'm using the ModelBuilder inside ArcMap to execute a set of raster calculations. I need to execute a different calculations based on the pixel value (If the pixel value is less then 0 - execute calculation no.2, If above 0 - execute calculation no.3).
Of course in the same raster can be both cases, so I need it to indeed chech every pixel and execute the appropiate following calculation per pixel.

Any idea how to do such thing ?


texts: Serving GeoTIFFs of different resolutions as single layer

texts:
I'm trying to serve a 30 arc second DEM with global coverage combined with high resolution 1/3 arc second DEM for a very small area.

I believed at first I could combine them into a single file, however this seems impossible as the file would end up as around 45GB

Questions

Is there another way of solving this problem using geoserver?
Alternative formats which support this?
Can the Image Mosaic plugin in Geoserver deal with this? I've tried, but it doesn't appear to work. The documentation implies only a single resolution can be dealt with, however the Heterogeneous=true option implies multiple resolutions?

Details of the two files

gdalinfo GEBCO_2014_1D_EPSG4326.tif
Driver: GTiff/GeoTIFF
Files: GEBCO_2014_1D_EPSG4326.tif
Size is 43200, 21600
Coordinate System is `'
Origin = (-180.000000000000000,90.000000000000000)
Pixel Size = (0.008333333333333,-0.008333333333333)
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-180.0000000,  90.0000000)
Lower Left  (-180.0000000, -90.0000000)
Upper Right ( 180.0000000,  90.0000000)
Lower Right ( 180.0000000, -90.0000000)
Center      (   0.0000000,   0.0000000)
Band 1 Block=43200x1 Type=Int16, ColorInterp=Gray

gdalinfo asc/puerto-rico.tif
Driver: GTiff/GeoTIFF
Files: asc/puerto-rico.tif
Size is 24301, 11881
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-67.600046296295005,18.800046283695000)
Pixel Size = (0.000092592590000,-0.000092592590000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -67.6000463,  18.8000463) ( 67d36' 0.17"W, 18d48' 0.17"N)
Lower Left  ( -67.6000463,  17.6999537) ( 67d36' 0.17"W, 17d41'59.83"N)
Upper Right ( -65.3499538,  18.8000463) ( 65d20'59.83"W, 18d48' 0.17"N)
Lower Right ( -65.3499538,  17.6999537) ( 65d20'59.83"W, 17d41'59.83"N)
Center      ( -66.4750000,  18.2500000) ( 66d28'30.00"W, 18d15' 0.00"N)
Band 1 Block=24301x1 Type=Int16, ColorInterp=Gray
  NoData Value=-9999


texts: How to assign values for one column of the attribute table on QGIS

texts:
The table is below.

Then I want to assign 3 values:

[43.8,
   39.0,
   25.9]

to the X field by row.

What should I do? If there are thousands of rows in the table, can I just paste them just like in Excel?

The version of my QGIS is 3.12.0


texts: Getting WFS layer from GeoServer to OpenLayers 3?

texts:
I am creating a web map application for a website. I am using GeoServer 2.6.2 and OpenLayers 3.2. I have not problems with WMS layers, but I am not be able to visualize vector layers. This is the javascript code that I have added (I followed this OpenLayers example and this other example):

var vector2Source = new ol.source.ServerVector({
  format: new ol.format.GeoJSON(),
  loader: function(extent, resolution, projection) {
    var url = 'http://localhost/geoserver/Workspace/wfs?service=WFS&' +
        'version=1.1.0&request=GetFeature&typename=Workspace:Polygon&' +
        'outputFormat=json' + '&srsname=EPSG:6823&bbox=' + extent.join(',') + ',EPSG:3857';
    $.ajax(url).then(function(response){
        vector2Source.addFeatures(vector2Source.readFeatures(response));
    });
},
strategy: ol.loadingstrategy.createTile(new ol.tilegrid.XYZ({
    maxZoom:19
  }))
  //projection: "EPSG:6823"
});

var vector2 = new ol.layer.Vector({
  source: vector2Source,
  style: new ol.style.Style({
    stroke: new ol.style.Stroke({
      color: 'rgba(0, 0, 255, 1.0)',
      width: 2
    })
  })
});

The browser (google chrome) shows me the following error:

XMLHttpRequest cannot load http://localhost/geoserver/Workspace/wfs?service=WFS&vers…7148101,1526294.5807983987,-9920914.775189597,1536078.5204189012,EPSG:3857. No 'Access-Control-Allow-Origin' header is present on the requested resource. Origin 'null' is therefore not allowed access.

Could someone help me?


texts: Why do I get an empty shapefile when converting this KMZ in QGIS?

texts:
I have downloaded a KMZ file containing global elevation information. I have tried to convert the data to a shapefile using qgis (save as shape option) and using the zonum solutions website but both times the conversion has not worked.

When using qgis to convert, the new shape file was empty and all attributes were 'NULL'.

Any idea on how to get around this? The KML/KMZ files I am using is from the Marine Geoscience Data system and is available from their website on the following link http://www.marine-geo.org/link/download.php?data_uid=14255&client=DataLink.


texts: Stop ArcGIS 'Geodatabase Adminstration' Reporting Connection as 'DBO' Instead of My Username?

texts:
We are running ArcGIS 10.5.1 and MS SQL Server, and using an SDE connection file which is configured to use 'Operating System' credentials.

When most users connect, the 'Geodatabase Administration' dialogue shows their Active Directory account correctly.  But when I connect using the same connection file, the same dialogue shows the connection username as 'DBO' (the MSSQL database owner user).

Now, I am a privileged user, including a member of the MS SQL 'sysadmin' role, so I'm sure it must be related to this.  But shouldn't it show the Windows/AD username that I am connected as?

Why does it show my connection as 'DBO' and not as the account with which I connected?

Here is what the connection properties looks like:

Below is the Geodatabase Administration Dialogue.  In this case, the smudged-out username is of the format '"Domain\User.Name"' and this user connected using the same connection file as myself, which it shows as 'DBO'.  (The 'GISADMIN'
 and 'SDE' connections used different connection files.)


texts: Why there is a mismatch of layers in R according to how the plot panel is resized?

texts:
Background

In R, I have downloaded two datasets using the following commands (slighlty modificed from https://mgimond.github.io/Spatial/point-pattern-analysis-in-r.html):

library(maptools)
z <- gzcon(url("http://colby.edu/~mgimond/Spatial/Data/starbucks.rds"))
location <- readRDS(z)

z <- gzcon(url("http://colby.edu/~mgimond/Spatial/Data/ma.rds"))
citylimits <- readRDS(z)

The first is of class SpatialPointsDataframe, while the second is of class SpatialPolygonsDataFrame.
I have also downloaded a SpatialGridDataFrame, which I have converted to a RasterLayer:

z <- gzcon(url("http://colby.edu/~mgimond/Spatial/Data/pop_sqmile.rds"))
popdensity <- readRDS(z)
popdensity <-raster::raster(popdensity)

Problem

Depending on the way in which (in RStudio) I resize the plot panel, it seems that the raster layer doesn't match the rest of the data (see attached images).

Question

Is that an RStudio issue, or is it an issue regarding the coordinate system of the layers? As for the latter, I have checked the CRS and the three are listed below (along with other info). I do not understand whether the raster layers' coord.ref exactly matches that of the two other layers. Is what described caused by different coord. ref. embedded in the layers?

popdensity

class       : RasterLayer
dimensions  : 308, 598, 184184  (nrow, ncol, ncell)
resolution  : 500, 500  (x, y)
extent      : 623157.2, 922157.2, 4602709, 4756709  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=18 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
data source : in memory
names       : band1
values      : 5.9, 61572.7  (min, max)

locations

class       : SpatialPointsDataFrame
features    : 171
extent      : 648032.3, 917740.7, 4609785, 4748107  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=18 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0

citylimits

class       : SpatialPolygonsDataFrame
features    : 1
extent      : 623157.2, 921923.8, 4602709, 4756659  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=18 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
variables   : 4


texts: Trying to compile QGIS but it complains it can't find QGIS library

texts:
I'm trying to compile QGIS in Windows7 with a MSYS, with the objective of making a QGIS Windows App using the API.  I am using this document.

I'm using CMake-GUI and it gives the error

  Warning Target "qgis_core" requests linking to directory
  "c:\OSGeo4W\bin".  Targets may link only to libraries.  CMake is
  dropping the item.]

Screenshot

I think its looking for qgis_core.lib, but I don't think it exists because I want to compile it... so this error makes no sense to me.


texts: Web AppBuilder for ArcGIS Query Widget spatial filter relational layer

texts:
I am working on a project that has to spatially analyze a point set to a parcel layer in ArcGIS online. Inside the Query Widget, I am using the spatial filter "Only return features that have a spatial relationship with features in another layer". I want to only allow for only 1 layer to be the relational layer, not every layer in the web map.

This is to minimize confusion for the end user and to generate results with limited interaction from them. I don't know if this is possible without editing the code.


texts: How to choose raster interpolation method?

texts:
I work with a dataset of ~1.3 million residential buildings; each of them is described by a normally distributed variable ranging from 0 to 100. The dataset covers whole Switzerland so there are obviously areas with both very dense and very sparse point density.

The limitation imposed by our data provider is that we cannot release datasets more detailed than 5 buildings.

I was thinking about converting this point data into a raster surface. My main goal in this exercise would be to create the surface that minimizes the error that will arise when somebody tries to get the point value from this raster by extracting data in their's point locations of interest (as compared to extracting it from the original dataset of point locations).

I'd like to know what kind of method should be best suited to do that. There are quite few to choose from in the Spatial Analyst toolbox and even more in Geostatistical one so I'd very much appreciate some help in where to start.


texts: Saving sublayers of route layer to disk using ArcPy?

texts:
I want to calculate a route consisting of 4 points with a script tool and save the results to disk. I need the route to appear as layer on ArcCatalog for continuing working with it (modification of its attribute table).
The script tool completes successfully but no output is saved.
So, I tried to add the Save To Layer File tool.

I receive the following error message:

  File "C:\Users\Manu\Desktop\SS 2015\Analyse mit GI-Technologien\Abgabe\gruppe1\Teil1_240915.py", line 84, in
      arcpy.management.SaveToLayerFile(outNALayer,outLayerFile,"","CURRENT")
   File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\management.py", line 6640, in SaveToLayerFile
      raise e
  ExecuteError: ERROR 999999: Error executing function.
  Failed to execute (SaveToLayerFile).

This is my code:

#Create a new Route layer.
outRouteResultObject=arcpy.na.MakeRouteLayer(in_network_dataset,out_network_analysis_layer=outNALayerName,impedance_attribute=impedanceAttribute,find_best_order="USE_INPUT_ORDER",accumulate_attribute_name=["Length"],output_path_shape="TRUE_LINES_WITH_MEASURES")
#Get the layer object from the result object. The route layer can now be referenced using the layer object.
outNALayer=outRouteResultObject.getOutput(0)

#Get the names of all the sublayers within the route layer.
subLayerNames=arcpy.na.GetNAClassNames(outNALayer)
#Store the layer names that we will use later
stopsLayerName=subLayerNames["Stops"]
routeLayerName=subLayerNames["Routes"]

fieldMappings=arcpy.na.NAClassFieldMappings(outNALayer,stopsLayerName)
fieldMappings["RouteName"].mappedFieldName="stopover"

arcpy.na.AddLocations(outNALayer,stopsLayerName,"Stations",fieldMappings,"","Order","SNAP")
arcpy.na.Solve(outNALayer)

arcpy.management.SaveToLayerFile(outNALayer,outLayerFile,"","CURRENT")


texts: How to add a map created in QGIS to a Matlab plot

texts:
I am generating a number of plots in Matlab. I want to overlay these over maps I have created in QGIS. I thought the best way would be to import the maps as a geotiff into Matlab and then overlay my plots over the image. I have saved the project as an image, generating a tiff and a tiffw file.

I got stuck trying to import and display the tiff and tiffw file in Matlab. I am working with QGIS 2.10.1 and Matlab R2012b.

Here is what i have tried so far:

Matlab has a geotiffread() function to read geotiff files. Using this function to read the .tiff produces the following error.

  Error using geotiffinfo>readinfo (line 259) The file 'GeoTest.tiff'
  does not contain any GeoTIFF Tags.  Use the function IMFINFO to obtain
  information about the file.

I was able to import the .tiffw file into matlab using the code.

filename = 'GeoTest.tiff';
RGB = imread(filename);
worldFileName = getworldfilename(filename);
R = worldfileread(worldFileName, 'planar', size(RGB))

But still have to figure out how to put that into a figure together with my plot.


texts: Applying same algorithm to various input files in QGIS modeler?

texts:
I wonder what´s the best way to create a ROI (by clipping algorithm) for all bands of one LANDSAT scene in QGIS graphical modeler.

Is there a way that avoids listing all (raster) bands and applying the algorithm to each?
something similiar to batch processing?


texts: ArcMap export map function not working properly

texts: I've got 3 different data frame's in my file'; all with a standard arcgis background. Whenever I try to export my layout view; the pdf/emp/png etc only loads one out of three background maps. What can be done about this?


texts: Extracting lines from polygons/buildings based on adjacency to lines/streets

texts:
I have a building and street shapefile of a district. I am looking for a way to select and eventually extract those lines or edges of the buildings that face the street.
however the extracted edges for each building need to be somehow grouped or merged into 1 line since on the next step i need to create certain number of points on the lines based on a value given for each building.
Is there anyway to do this?


texts: Why does my arcpy script's da.UpdateCursor get slower the longer it runs?

texts:
I have a layer containing pedestrian crosswalks (line features), and two layers containing street data: street segments (line features) and street nodes (point features). Every street segment has a node at either end, with its own Node ID. Most of the crosswalks cross one or two street segments. For each point at which a crosswalk crosses a street segment, I'm trying to determine which node ID is closest. I'm using this method:

Generate the intersection points between the crosswalks and the street segments.
At each of those points, split the street segment.
Spatially join the node ID to each of the two (now split) street segments.
Find the shorter of the two segments.
Return the node ID of that shorter segment.

Because a street segment can have many crosswalks, steps 2-5 have to be done for each crosswalk/segment intersection point individually. Part of my script is below, it takes the output of step 1 above and runs a da.UpdateCursor to do steps 2-5 for each crosswalk/segment intersection point. The script works, it generates the correct output.
I'm writing and running the script in PyCharm, using the Python 3 environment that comes with ArcGIS Pro. The problem is that the longer it runs, the slower it gets. It starts out processing one row in 0.7 seconds. After 1000 records, it's up to 4 seconds, getting slower and slower until I stop it manually. Stopping and restarting the script in PyCharm resets the performance.
I've read that nesting a da.SearchCursor inside another cursor can cause this progressive slowdown, and I'm using a nested cursor here, but I'm not clear how that would be causing the slowdown in my case: it's only ever searching through 2 records and I'm deleting the SearchCursor after each iteration of the UpdateCursor.
A few other things:

I have turned off logging using arcpy.SetLogHistory(False).
PyCharm's RAM and CPU usage remain about the same throughout: about 700MB of RAM used and ~25% of the CPU, with 6GB RAM remaining free on the PC and > 30% of CPU capacity remaining. I've also tried running the script in IDLE and Spyder with the same result, so it doesn't seem specific to the IDE.
Temporary data is using the memory\ workspace. The datasets that are on disk are on a solid-state drive.

# seg_split_at_crosswalk, seg_split_with_nodes are feature classes in the "memory\" workspace.
# segment_fc, node_fc are existing data on disk.
# crosswalk_seg_crossings is a point layer of the intersections between crosswalks and segments, created earlier in the script.

# for each crossing point, find closest node which is on the segment crossed, and update row with NODE_ID
with arcpy.da.UpdateCursor(crosswalk_seg_crossings,["SEG_ID", "SHAPE@", "NODE_ID", "OBJECTID"]) as cursor:
    for row in cursor:
        if row[2] is not None:
            continue
        elif row[2] is None:
            loop_start_time = time.time()
            sql = "SEG_ID = " + str(row[0])
            segment = arcpy.management.SelectLayerByAttribute(segment_fc, "NEW_SELECTION", sql)
            arcpy.management.SplitLineAtPoint(segment, row[1], seg_split_at_crosswalk)
            arcpy.analysis.SpatialJoin(seg_split_at_crosswalk, node_fc, seg_split_with_nodes)

            # select shortest split seg
            dist_to_nodes = dict([(searchRow[0], searchRow[1]) for searchRow in arcpy.da.SearchCursor(seg_split_with_nodes,["NODE_ID", "SHAPE@LENGTH"])])
            row[2] = (min(dist_to_nodes, key=dist_to_nodes.get))
            cursor.updateRow(row)

            # print stats
            loop_execution_time = (time.time() - loop_start_time)
            print(sql, "|", row[2], "|", row[3], "|", loop_execution_time)

            # delete temp data
            arcpy.management.Delete(seg_split_at_crosswalk)
            arcpy.management.Delete(seg_split_with_nodes)
            del dist_to_nodes
            del loop_execution_time
            del loop_start_time
            del sql
            del row


texts: Specifying PROJ.4 projection with North pointing down?

texts:
I would like to create a world map "upside down", with the South on top (in Mapnik).

How can I achieve this?

Can I take an existing projection string, e.g. Robinson projection, ESRI:54030.

+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs

and flip/mirror it? So the North will point down.


texts: How do I access a feature layer's "displayField" from a REST SOE?

texts:
Understanding I can't use the 'ESRI.ArcGIS.Carto.IFeatureLayer' interface, how can I get to a Feature Layer's "displayField" from a REST server object extension using arcobjects C#?

FeatureLayer.DisplayField


texts: QGIS + PostgreSQL + row-level-security + SET ROLE in multi-user environment

texts:
Is it possible to execute SET ROLE from within QGIS?
What I want to achieve is to create individual roles for users, that should have access to the same part of data within tables, via existing RLS policies which use CURRENT_USER against a specific table column values. Then use SESSION_USER to audit the tables so that we can keep track of who exactly made some changes.
I can execute it from within DB Manager but it seems to only apply there. I need an easy for a user way to execute it so that it applies for QGIS-browser (so that a user can actually see the tables that he has access to after SET ROLE).
Unfortunetaly PostgreSQL doesn't have any built-in AFTER LOGON triggers so I can't handle it using only backend. I know that in theory I could hardcode RLS policies so that let's say roles A, B, C, D, E all see one anothers updates and inserts but that would be a nightmare to maintain.


texts: arcgis C#: restrict addin to mxd file only

texts:
I have an addin toolbar (arcmap 10.3 in C#) for a certain type of analysis for our clients. We want to restrict the availability of the addin to a mxd file. We do not want other COM or Addins to be able to 'enable' it.

Is this possible?


texts:
I have been struggling with something since yesterday and Google did not help much...

I am working with two raster images, same projection (I re-projected one of them, they were originally different), same resolution and same extent. BUT the cell grids won't align. I am attaching a picture for you...

I tried a couple of things, like the align tool on QGIS and also re-sampling. The re-sampling "worked", the cell grid got aligned but the values of my DEM raster got all messed up... Is there any way I can "force" than to align? Any advice, tips?


texts:
I have a really simple question I can't find the answer to. I'm looking to delete a range of values from an attribute table in ArcGIS 10. I'm looking for a quick way to delete about 3000 records. The idea is to have some blank records so when I do a dissolve in order to form polygons, the 3000 tiny polygons are no longer there and only larger polygons are dissolved to form bigger polygons on the basis of the same column contents.

So I basically have 2 questions:

How do I delete a range of 3000 records from the attribute table, by FID range or any other efficient method?
Will blank cell contents prevent records from being included in a dissolve?

(Or 3., is there any way of creating a conditional dissolve which avoids a range of records according to FID/other column value?)


texts: QGIS: How to set up a datasource connection to pull locally sourced layers when the name changes

texts:
I am setting up a template QGIS project which pulls layers from the folder the template is in. The template needs to connect to a raster layer that has a name that slightly changes each time. The name changes depending on the property, ie property_NCAS.tif . I want to save a connection that caters for this and connects regardless of the property name.
Currently every time the project opens a pop up window titled "handle Unavailable layers" pops up and requires you to edit the datasource "./project_NCAS.tif" and apply the changes. Is there a way to skip this step?


texts: Comparing point pattern and raster data?

texts:
I have a simulated data set showing distribution of a certain habitat type as a point pattern.

Then I have a empirical data set but there the format is raster.

The two datasets are covering the same area.

How can I compare these two data sets?

I want to know if the simulated data is a good representation of the empirical one.

I use R for the analyses.


texts: Interpolate multiple point files using the python console in QGis

texts:
Using the python console in QGis, I would like to interpolate several point files (txt x,y,z) to make rasters.

I have been trying this but it's not working and I don't know enough about python to make it work.

import qgis.core
import qgis.analysis

uri="file:///C:/test.txt?delimiter=,&xField=field_1&yField=field_2&crs=epsg:32755"
layer = QgsVectorLayer(uri, 'test', "delimitedtext")
layer_data = qgis.analysis.QgsInterpolator.LayerData()
layer_data.vectorLayer = layer
layer_data.zCoordInterpolation=False
layer_data.interpolationAttribute =2
layer_data.mInputType = 0

tin_interpolator = QgsTINInterpolator([layer_data])

export_path = "C:/test.asc"

rect = layer.extent()
res = 10
ncol = int( ( rect.xMaximum() - rect.xMinimum() ) / res )
nrows = int( (rect.yMaximum() - rect.yMinimum() ) / res)
output = QgsGridFileWriter(tin_interpolator,export_path,rect,ncol,nrows,res,res)
output.writeFile(True)

I keep getting this error: Traceback (most recent call last):
  File "", line 1, in
  File "C:/test.py", line 13, in
    tin_interpolator = QgsTINInterpolator([layer_data])
NameError: name 'QgsTINInterpolator' is not defined

Can anyone explain what this variable means layer_data.zCoordInterpolation=False, and whether I should put this to True?


texts:
I am having troubles interpretating the output of r.cost.full.

Currently i am working on a project comparing some arcgis tools with qgis tools (workflow and output). One of my examples is to make a cost distance raster with both softwares.

The input Cost Raster and the Input Start Points (point features) are the same in both softwares. My problem here is that the result values of arcgis are within 0 and ~7 which make sense in terms of interpretating. But the result values of r.cost.full are between ~0.02 and ~0.3. When i visualize the output rasters with the same amount of classes, the rasters look very similar.

I ran the calculations on arcmap 10.2 and 10.3 and on qgis 2.6 and 2.8.2

Does anyone know why the values are so different? Am i missing something obvious here?


texts: How to batch rename tables in GeoPackage with OGR

texts:
I would like to rename all tables in a GeoPackage file. The "translations" are stored in a .csv file.
I would like to do this on the command line and found the below thread that uses the following command:
ogrinfo my_geopackage.gpkg -sql "ALTER TABLE old_layer_name RENAME TO new_layer_name"

(from: Rename layer in gpkg using QGIS or PyQGIS?)
How do I modify the above so that old_layer_name and new_layer_name are fetched from the .csv file? Keep in mind that I am a noob.


texts: NetCDF projection conversion to EPSG3857 is cropped in QGIS

texts:
I am using QGIS and PyQGIS to transform NetCDF files projection to EPSG:3857. My ultimate goal is export it as raster (TIFF or PNG) to display in Mapbox.
When first imported, the NetCDF file itself does not have a projection so I set the layer's CRS to EPSG4326 while the project CRS is EPSG:3857. It seems to have converted correctly, however, when a world shape file is overlaid, the north above Greenland is cropped. I understand that the because EPSG:3857 projection outputs a square, the latitude only goes from -85 to 85.
Is there a way to manually configure so it outputs correctly like the second image? (produced from Panoply).


texts:
I have shape vector layer in ESPG:3859 CRS. I try to get layer crs

layer.crs.geographicCRSAuthId()

but method always returns EPSG:4326 no matters which CRS I chose in layer properties.

.qpj file for shape layer:

PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],AUTHORITY["EPSG","3857"]]

.prj

PROJCS["WGS_84_Pseudo_Mercator",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator"],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1],PARAMETER["standard_parallel_1",0.0]]


texts: Adding KML as operationallayers on local server

texts:
I want to add a KML file, hosted internally on a server, as operationallayers. It doesn't seem to be working.

If I put the KML file in Google Drive or Github and then call it like :

{
            type: 'kml',
            url: 'http://googlemaps.github.io/js-v2-samples/ggeoxml/cta.kml',
            title: 'Quakes2',
            options: {
                id: 'quake_kml',
                opacity: 1.0,
                visible: true,
                //imageParameters: imageParameters
                mode: 0
            },
            identifyLayerInfos: {
                //layerIds: [0]
            },
            legendLayerInfos: {
                layerInfo: {
                }
            },
            layerControlLayerInfos: {
                swipe: true,
                metadataUrl: true,
                expanded: true
            }
  }

. It works. If i change to "url: 'http://localhost/mykml/cta.kml'", it doesn't work. But I don't want to put my resource on outside premises.

How can I put this KML file on my internal network so that I can call it in my app?

I have read about KML parser called geoxml3, but I don't know how to convert the result of geoxml3 so that it can be passed to my CMV app. Any ideas?


texts: Latitude/Longitude to city and zip, batch convertion

texts:
I have ~10 mln lat/lon points, and I need to extract city and zip code from them.

I can do this with OpenStreetMap using the Nominatim API (using something like http://nominatim.openstreetmap.org/reverse?query) through geopy.

The problem is that the API takes only one lat/lon pair at a time, which, with 1 second delays between each query, would take ~3000 hours in total or maybe more.

Is there a similar API where I can submit a batch request with a few thousand pairs at the same time? Or something that I can (easily) set up locally that would give me the same information?

Also, I know that it's possible to set up Nominatim locally. However it seems quite complex, and I need only a part of the data, not everything. And it's a one-time job, so probably I won't need this in the future.

Some extra details: the lat/lon pairs are from Russia, the results I get from OSM are very good, so I'd like to have the same or something similar.


texts:
I'm using the postcode boundary polygons from here: http://www.opendoorlogistics.com/downloads/

As the data is reconstructed, it has glitches here and there.

SELECT geom FROM postcodes WHERE name = 'SW1A'; results in 01060000000300000001.

Mapped, this looks like

The top two polygons are correct but the bottom blob is completely incorrect.

What would be the most efficient method to remove this?


texts: What coordinate system is this? How can I project it to WGS84 (lat-long) system in Python?

texts:
I had an AutoCAD drawing of a floor plan. I imported it to python using Fiona and checked the metadata. The CRS is empty. When I checked the coordinates for the geometry, it was in the following form,
{
    "type": "FeatureCollection",
    "features": [{
        "type": "Feature",
        "id": "0",
        "properties": {
            "A-WALL__ID": null
        },
        "geometry": {
            "type": "LineString",
            "coordinates": [
                [5319.568084840497, 5864.945849450571, 0.0],
                [5319.568084840497, 5882.945849450571, 0.0]
            ]
        }
    }, {
        "type": "Feature",
        "id": "1",
        "properties": {
            "A-WALL__ID": null
        },
        "geometry": {
            "type": "LineString",
            "coordinates": [
                [4767.568084840497, 5863.945849450571, 0.0],
                [4767.568084840497, 5882.945849450571, 0.0]
            ]
        }
    }
]
}

I imported the AutoCAD file in following way,
c = fiona.open("level0dwg.dxf", 'r')
c.meta

Following was the output,
{'driver': 'DXF',
 'schema': {'properties': OrderedDict([('Layer', 'str'),
               ('PaperSpace', 'int'),
               ('SubClasses', 'str'),
               ('Linetype', 'str'),
               ('EntityHandle', 'str'),
               ('Text', 'str')]),
  'geometry': 'Unknown'},
 'crs': {},
 'crs_wkt': ''}

I am confused as in what is this coordinate system. Is it pixel-based or something like that?
and how can I transform it to WGS84 system with lat-longs to place it on the map?
The AutoCAD file is not georeferenced.
Is there a way to do this in Python?
I want to georefer the DWG file to place it on the map using Python. I have been stuck at this issue for quite some time now.


texts: How to set the minimal value of a layer to be displayed in google earth engine?

texts:
so i'm using OpenLandmap Sand Content and trying to show layer of an area only with average soil sand content above 50%;
here's my attempt

    var areas = ee.FeatureCollection('users/sdwijati/IDN_adm2')
    var nama_prov = ['Jawa Barat', 'Jakarta Raya', 'Jawa Tengah',
                      'Jawa Timur', 'Yogyakarta', 'Banten']// 'Sulawesi Tengah']
    var area = areas.filter(ee.Filter.inList('NAME_1', nama_prov));
    //---------------------------------------------------------------------------//
    var sand = ee.Image("OpenLandMap/SOL/SOL_SAND-WFRACTION_USDA-3A1A1A_M/v02");
    var sandavg = sand.reduce(ee.Reducer.mean());
    var ltsandavg = sandavg.mask(sandavg.lte(50));
    var sandvis= {min: 0, max: 100, palette: ['fbff00','b4ff00','00c4ff','0005ad']};
    Map.addLayer(ltsandavg.clip(area), sandvis)
    Map.centerObject(area, 7)

but instead it only shows the area with sand content below 50%, any idea on how to do the opposite?


texts:
I have a problem with ArcMap that somehow affect georeferencing.

So I have three layers, two with .lyr and one is a raster data.

The two .lyr datatypes are a border and cities with names. I georeferenced the raster
using the border and the cities and everything went well. But what I wanna do now is
to combine the two .lyr datasets with my raster so I can export the Raster with the overlayed
border and the cities as a pdf.

Is that even possible???

I already tried dissolve, merge, package etc. but none of them seems to work. I would be really grateful if someone could help me with it or at least let me know that it's not possible so I don't need to waste any more time on this topic.


texts: QGIS attribute data updating by not commiting

texts:
I am trying to change the values of an attribute in one of my layers. My code is as follows:

from qgis.utils import iface
from PyQt4.QtCore import *
layer =iface.activeLayer()
pr=layer.dataProvider()
layer.startEditing()
for f in layer.getFeatures():
    geom=f.geometry()
    xval= geom.asPoint().x()
    f['XX']=xval
    layer.commitChanges()
    print f['XX']
for f in layer.getFeatures():
    print f['XX']

For my first print statement, I am getting the values out that I want, but when I run the second print statement, the values are what they were before running the code.

Am I Not committing changes correctly?


texts: Set workspace to SDE version using ArcPy?

texts:
My company uses a versioned SDE database to perform edits, update features, etc.

I am working on a script that would inevitably update features and would like to set the environment to a specific version.

The environment set to the SDE would look like this:

arcpy.env.workspace = "C:\\Users\\MyName.CENTRAL\\AppData\\Roaming\\ESRI\\Desktop10.3\\ArcCatalog\\Production.sde

My version is called "SY.TEST"

How do I define this version (only) as the workspace?


texts: How to set a polyline length as per a given Image

texts: I have an image of a map representing the road networks. I have converted these road networks into a shapefile with only polylines. I want to set a distance to this road network. For e.g I know from the image that the distance between A and B is 10 Km. However, I want to set the same distance to the polyline also. Because. If i draw a buffer of 2km for this road network, it goes beyond the size of the image. Any idea how to do it?


texts: Strange reflectance values in Sentinel-2 images over sea ice

texts:
I am dealing with Sentinel-2 (Level 1C) images over sea ice. I have downloaded some and was wondering: When I look e.g. at Band 8A, there are some images where the reflectance over the water between the ice has realistic values of 0.4 (this is at the same time the lowest reflectance value in this image). But in other images, the whole water has a reflectance of almost zero, partially even exactly zero which is a rather unrealistic reflectance for water.

How can this be and is there a possibility to correct the low values?


texts:
I'm applying a function that distributes 500 points within the FeatureCollection.`

var getSamples = function(image, samplesFcs, samplesPoints, nsamples, label){

  var samplesCollection = ee.List(samplesFcs).iterate(
    function(sample,totalSample){
        return ee.FeatureCollection(totalSample).merge(sample);
      }
    ,ee.FeatureCollection([]));

    samplesCollection = ee.FeatureCollection(samplesCollection)
    return samplesCollection

    var randomPoints = ee.FeatureCollection.randomPoints(samplesCollection, nsamples, 1989)

      samplesPoints.forEach(function(points) {
        randomPoints = randomPoints.merge(points);
      })

     return randomPoints
}

however it is wrong, after all, when I give a print on the console to check how many points were sampled inside the Features, the number returned is about Features within the FeatureCollection, as opposed to the number of samples

var test = getSamples(addndvi_2018,class4,class4Points,nsamples)
print(teste.size(),'samples Vegetation')

Can anyone point to where I'm doing wrong here?


texts: Camera reports database online to look it up?

texts:
I received tons of aerial photos that go back to 1930's and I have the information about the camera number and lens number and focal length. They were taken in a picture of the roll that tells us what it is.  What I would like to know is what they are and what project they were flown in and I want to get more information about them.

I was wondering if there is a camera calibration report database that I could look up online. Can anyone point me to where or who does that?


texts:
In my JS web-application, I'm trying to store the tiles in my localstorage using web workers.I'm using the code from this Esri sample. I just deleted line 20 in the ioWorker.js file:

 url = "/arcgisserver/apis/javascript/proxy/proxy.ashx?" + url;

But when the web worker requests a tile, I get those three errors:

On my side, CORS should be working fine. Is this a problem on Esri's side or have i got something wrong?
edit: Or is it not possible only with CORS, do I need a proxy-page?


texts: Calculating haversine versus checking latitude and longitude digits after decimal point?

texts:
If I want to tell whether 2 points are close to each other, say at proximity of 50 meters, is it OK to compare both numbers up to a certain decimal point, for example the two points:

  36.391155, -56.095878 and 36.391187, -56.095899

then checking if both points latitude and longitude are equal at the 4th decimal point, it's safe to say they are close, or there are cases where it won't work, and I must calculate haversine?


texts: Moving between pages in QGIS Print Composer

texts:
I'm using QGIS 3.6 Print Composer, and I managed to create a layout with two pages.

I now want to delete page 1. How do I move from page 2 to page 1 so I can delete it?


texts: Creating multiple path distance outputs using ModelBuilder?

texts:
I am trying to create a model in ArcMap that will run a path distance analysis on each individual polygon in a feature class.

Below is the model I have so far but it is only outputting one path distance rather than an individual one for each of the polygons in my feature class.

I have selected 'Id' as the field in the Iterator (each polygon has a different ID number'. The output from the iterator is my input feature in my path distance tool and then the 'ReCamp2' is my input cost raster.

Do I need another tool to get my path distance to save the results from each polygon separately?


texts: Get ID for each single v.to.points point in GRASS?

texts:
How do I get for each v.to.points point an ID in the database (DB).

So far I splited the line (each 10m) but it returns only one DB entry for each original line.

v.db.select line
cat|Id|Name
1|1|Line 1
2|2|Line 2

# split the lines every 10m
v.to.points in=line out=line_points dmax=10

# check the result
v.db.select line_points
cat|Id|Name
1|1|Line 1
2|2|Line 2

Each point should have its own DB entry. The result should be:

v.db.select line_points
cat|Id|Name
1|1|Line 1
2|1|Line 1
3|2|Line 2
4|2|Line 2

Workaround: So far my workaround is, to export (v.out.ogr) the points as shape and import (v.in.ogr) it again. Then each point has his own ID. But there must be a smarter way to do that!


texts: PyQGIS script crashes QGIS 3 when remove a custom pluginLayer which has custom legend

texts:
I upgrade a 2.18 plugin to 3.4.

PluginLayer has been refactored and I could not create a custom legend as simple as previous major version of QGIS 2 :

  QgsDefaultPluginLayerLegend was removed. Use QgsMapLayer::setLegend()
  to provide legend nodes for plugin layers.
  as it's wrote here https://qgis.org/api/api_break.html

Well I followed steps :

create my customed pluginlayer
add to mapLayer

QgsProject.instance().addMapLayer(layer)
legend = layer.makeLegend()
layer.setLegend(legend)

Here my pluginLayer class :

class PSudLayer(QgsPluginLayer):

def __init__(self,session, type, name):
    QgsPluginLayer.__init__(self, type, name)
    self.layerType = type
    self.session=session
    self.conf = None
    self.transparency = 0
    self.setValid(True)

def makeLegend(self):
    return PSudLegend(self.session,self.conf['uuid'])

and here may customed legend

class PSudLegend(QgsMapLayerLegend):

    def __init__(self, session, uuid):
        QgsMapLayerLegend.__init__(self)
        self.nodes = []
        url = session.baseUrl+'/geocat/legend/'+uuid
        r = session.get(url)
        self.legends = json.loads(r.readAll().data().decode())
        self.globalcount = 0
        self.countBoucle = 0
        self.countRecursive = 0
        self.icon = QIcon()

    def recursiveLoad(self,legends, nodeLayer):
        self.countRecursive +=1
        self.countBoucle = 0
        size = legends['children'].__len__()
        for child in legends['children']:
            if 'children' in child:
                self.recursiveLoad(child, nodeLayer)
            else:
                if size > self.nodes.__len__(): # hack car la boucle passe plus d'une fois ....
                    image = QPixmap()
                    iStr= child["icon"].split(",")[1]
                    i64 = base64.b64encode(iStr.encode()) # hack parsing error (multiple de 4)
                    image.loadFromData(QByteArray.fromBase64(base64.b64decode(i64)))
                    self.icon = QIcon(image)
                    node = QgsSimpleLegendNode(nodeLayer,child['title'],self.icon)
                    self.nodes.append(node)

    def createLayerTreeModelLegendNodes(self, nodeLayer):
        if 'children' in self.legends:
            self.recursiveLoad(self.legends, nodeLayer)

        return self.nodes

    def removeLegend(self):
        for node in  self.nodes:
            node = None
        self.nodes = None
        self.icon =None
        return True

All work fine as I remove the layer from the treeLayerPanel.  QGIS simply crash.

I tried to remove nodes in my legend with no success and call it when removeLayer is called.

QgsProject.instance().layerWillBeRemoved.connect(self.PSudData.removeLayer)

def removeLayer(self,layerid):

    for layer in self.parent.service_layers:
        if layer.id() == layerid:
            layer.legend.removeLegend()
            QgsProject.instance().removeMapLayer(layer)

I get a typeError when I tried the filtering option of the layerTree.
PSudLegend.createLayerTreeModelLegendNodes(), index 0 has type 'QgsSimpleLegendNode' but 'QgsLayerTreeModelLegendNode' is expected
Which I don't understand cause QgsLayerTreeModelLegendNode is an abstract class and QgsSimpleLegendNode extended from QgsLayerTreeModelLegendNode.

When I tried to refreshLayerLegend() from the LayerTreeModel and I also get the same TypeError.

Very confused.


texts: Using QGIS to extract elevations of specific points from a DEM layer

texts: With QGIS My objective is to extract Elevation Values from a DEM dataset. The locations to extract elevations come from a segp1 file. My DEM values are georeferenced to a UTM, my segp1 is georeferenced to a state plane system, I want output of extracted Elevations to the State Plane system?


texts: Encoding error between maptools and ArcGIS

texts:
I am having a trouble of encoding at writing shapefile by maptools (0.8-34) of R by MacOSX and reading that by ArcMap 10.2.1 in Windows 8.1 through parallels desktop. Using an original shapefile, I wrote an edited file by writeLinesShape as follows:

library(maptools)
original <- readShapeLines("original.shp")
writeLinesShape(original, "edit.shp")

After this, I try to read "edit.shp" with ArcMap but cannot read Japanese strings properly in "dbf" file, which I could read in the original file. (Of course I did some edit before writing "edit.shp" but I have the trouble even in this simplest code.)

The original shapefile seems utf-8 because

iconv(original$japanese,from='utf-8',to='shift_jis')

worked in R. Also, I tried the following suspecting that ArcGIS can read shift_jis or cp932,

iconv(original$japanese,from='utf-8',to='shift_jis')

and

iconv(original$japanese,from='utf-8',to='cp932')

before writing the file, but Japanese strings were not shown properly in ArcMAP again.

Is there any way to solve this problem? I suspect that maptools is not writing utf-8 file because ArcMap is able to read utf-8.


texts: Splitting lines at polygon intersection with attributes saving using PostGIS

texts:
This relates to Splitting line layer at polygon's boundary of polygon layer using PostGIS but what about keeping the attributes.

Having

line_layer (id, geom, attribute1, attribute2)
poly_layer (id, geom, attribute3, attribute4)

How do I split the lines after polygons and preserve the attributes for each segment?

My current code is:

SELECT row_number() OVER () AS id, *
   FROM ( WITH lines_in_polygons AS (
                 SELECT st_intersection(ln.geom, a.geom) AS inter
               FROM line_layer ln,
                polygon_layer a
            ), diff AS (
             SELECT st_difference(ln.geom, st_union(pl.geom)) AS geom
               FROM line_layer ln
          JOIN polygon_layer pl ON st_intersects(ln.geom, pl.geom)
         GROUP BY ln.geom
            )
             SELECT st_geometryn(d.geom, n.n) AS geom
               FROM diff d
         CROSS JOIN ( SELECT generate_series(1, st_numgeometries(diff.geom)) AS n
                       FROM diff) n
    UNION
             SELECT lines_in_polygons.inter
               FROM lines_in_polygons) coco
  WHERE geometrytype(coco.geom) = 'LINESTRING';


texts: Using Fontawesome 5.7 icons with OpenLayers 3

texts:
I'd like to use this Fontawesome version with Openlayers 3. I tried changing the font value from FontAwesome to Font Awesome Free, but nothing happens. On the map appears a square (it means that the unicode cannot be found inside the font).
return new ol.style.Style({
  geometry: geom,
  text: new ol.style.Text({
    text: '\uf041',
    font: '20px FontAwesome',
    // font: 'Font Awesome Free',
   })
});


texts: Merging agricultural rasters using ArcGIS Spatial Analyst?

texts:
I am working in ArcMap 10.3 and I would like to merge two rasters, one of wheat harvest area, the other of corn harvest area. Looking at the corn rasters table shows the amount of corn harvested in that pixel. I would like to combine the two and get a value for the highest.
I am thinking something like set corn to 1 and wheat to 2. Whenever corn is higher than wheat assign a 1 and so on.  The output raster then have every pixel showing a 1 or 2 to represent which of the two was higher in that location. Both rasters have the same spatial extent

Can this be done?

I have been reading through the tool list and have not found any tools that look like they will do the job.


texts: Extracting all field names of a featureclass using query filter programmatically

texts: I want to get all field names present in a feature class using queryfilter programmatically.


texts: Add an attribute name as parameter - QGIS 3 modeler

texts: I'd like make a model where I can change an attribute name to a predefined one using refact tool. How can I get that "attribute" by a parameter from a vectoral field ?


texts: What is the difference between geocoding and georeferencing?

texts:
After some googling around both seem to be interchangable words, where we speak of tying the spatial info of a real place and storing in database.

But still I am not satisfied, there should be some major differences. Also are there any tools in Oracle Spatial suite for this to carry out.


texts: Modifying DEM elevation data in QGIS using a polygon

texts:
I've tried searching through the existing questions but am struggling to find a succinct answer to what I am trying to achieve.

For background, I am running a flood model in TUFLOW to help support a residential development on greenfield land. Part of the site has been shown to flood in a given event and I now want to reprofile areas of the land to allow additional flood storage. The hope is that with the reprofiled land, the flood extents will be greatly diminished and so there will be more land available to the development.

I have drawn polygons in QGIS over the exact areas that will need to be lowered by specific amounts i.e. the land in one area varies from 66.42-66.73m AOD. This will now be a uniform 66.10m AOD.

Is there a way of simply modifying the existing raster DEM file using the polygons?


texts: QGIS 2.14 Couldn't load plugin 'processing' due to an error when calling its classFactory() method

texts:
I understand that my problem originates from the letters in my username which apparently are non-ascii and very Nordic. I just have absolutely no clue how to get around this. I'm on Windows 10 64bit and QGIS is 2.14 Wien, not 2.14.1.

I don't really want to create a separate user account only for QGIS. Is there a way I can make QGIS "understand" my filepath without changing my username's lettering. What's weird is that I've not had any problems using QGIS with this error message popping up every time I start the program. Every analysis I've churned out has checked out accurately. The headache this error has given me has been purely cosmetic as far as I know.

So is there a pain-free workaround?

Here's my error message poem:

2016-05-02T16:19:19 1   Traceback (most recent call last):
              File "C:/PROGRA~1/QGISES~1/apps/qgis/./python\qgis\utils.py", line 306, in startPlugin
                plugins[packageName] = package.classFactory(iface)
              File "C:/Users/V�iski/.qgis2/python/plugins\processing\__init__.py", line 37, in classFactory
                from processing.ProcessingPlugin import ProcessingPlugin
              File "C:/PROGRA~1/QGISES~1/apps/qgis/./python\qgis\utils.py", line 572, in _import
                mod = _builtin_import(name, globals, locals, fromlist, level)
              File "C:/Users/V�iski/.qgis2/python/plugins\processing\ProcessingPlugin.py", line 47, in <module>
                if cmd_folder not in sys.path:
              File "C:/PROGRA~1/QGISES~1/apps/qgis/./python\qgis\utils.py", line 58, in showWarning
                "warning:%s\ntraceback:%s" % (warnings.formatwarning(message, category, filename, lineno), stk),
            UnicodeDecodeError: 'ascii' codec can't decode byte 0xe4 in position 18: ordinal not in range(128)


texts: How to convert an insert SQL query to a PostGIS layer?


texts: How to use ICalculator with a pre-expression in C#

texts:
Another wonderful predicament. I am trying to use an instance of ICalculator to update a field based on another field. Textbook stuff. Unfortunately, the example in the SDK for C# is a bit lacking.

Here is what I have puzzled out so far from the C# GeoDatabaseUI 10.x and VB.NET ICalculator 9.3 examples/documentation (NOTE: I am doing this in C#, in case I was not clear).

If you want to use a string in the Expression property, do not use triple quotes (must be legacy information), but escaped quotes for normal string-within-string notation.
You can use VB.NET code in the pre-expression when working in VB.NET (this confuses me).

With these assumptions in mind, I have been attempting to use python or C# in the pre-expression property, but I am having no luck.

snip....

Has anyone tried this yet? If so, what is the proper way to use a python function and/or a C# method when working with the ICalculcator?

UPDATE: Okay, trying the method kenbuja suggests below, I changed it to an IGeoProcessor method with the following code:

public void GPCalculate(ITable table)
    {
        IGeoProcessor gp = new GeoProcessorClass();
        IVariantArray varray = new VarArrayClass();
        string input = workspacepath+"\\"+targetfile;

        string codeblock = "states = {'2':'AK', '1':'AL', '5':'AR', '60':'AS', '4':'AZ', '6':'CA', '8':'CO', '9':'CT', '11':'DC', '10':'DE', '12':'FL',"+
            "'13':'GA', '66':'GU', '15':'HI', '19':'IA', '16':'ID', '17':'IL', '18':'IN', '20':'KS', '21':'KY', '22':'LA', '25':'MA', '24':'MD', '23':'ME',"+
            "'26':'MI', '27':'MN', '29':'MO', '28':'MS', '30':'MT', '37':'NC', '38':'ND', '31':'NE', '33':'NH', '34':'NJ', '35':'NM', '32':'NV', '36':'NY',"+
            "'39':'OH', '40':'OK', '41':'OR', '42':'PA', '72':'PR', '44':'RI', '45':'SC', '46':'SD', '47':'TN', '48':'TX', '49':'UT', '51':'VA', '78':'VI',"+
            "'50':'VT', '53':'WA', '55':'WI', '54':'WV', '56':'WY'}\n" +
            "def getstate(id):\n" +
            "  id = str(id)\n"+
            "  if (len(id) == 14):\n" +
            "    id = id[0]\n" +
            "  if (len(id) == 15):\n" +
            "    tid = id[0:2]\n" +
            "  return states[id]";
        // First Parameter in_table
        varray.Add(input);
        // Second Parameter calculate field
        varray.Add("PSTATE");
        // Third Parameter expression
        varray.Add("getstate(!ID!)");
        // Fourth Parameter expression_type (optional)
        varray.Add("PYTHON_9.3");
        // Fifth Parameter Codeblock (optional)
        varray.Add(codeblock);

        IGeoProcessorResult result = new GeoProcessorResultClass();
        settesttext("Starting Calculation");
        result = gp.Execute("CalculateField_management", varray, null);
    }

The above code works perfectly now (thanks to an edit @blah238).


texts:
Sometimes I think about possible analysis using GIS models (I use QGIS 1.8), and I always think about few types of visualization techniques: demographics, traffic flows, distance ones...

Thanks to you and reading those pages I discover many possible criteria (and I begin to use it) analysis and I like to collect the most important one in a single item to help the consultation.

Is this a good idea or not?

I try with some examples of a project where I worked on

generated traffic by zone by magnitude

traffic IN/OUT zone by zone by bar diagram

journey time on different routes

OD from a zone (thanks to Flow Mapper plugin my question on it )

quantity and quality parks on links


texts: How to reverse the MultiLineString command to a list of arrays


texts: Adding local GeoJSON files to Leaflet web page

texts:
I'm quite new to JavaScript and I have been trying to add my GeoJSON files to a web page. I can't seem to figure out how to get the GeoJSON to actually load. I've checked it using GeoJson.io and it draws correctly there. I just can't get it to load on my page. The OSM portion works just fine no issues showing that.

L.tileLayer('https://api.tiles.mapbox.com/v4/{id}/{z}/{x}/{y}.png?access_token=pk.eyJ1IjoibWFwYm94IiwiYSI6ImNpejY4NXVycTA2emYycXBndHRqcmZ3N3gifQ.rJcFIG214AriISLbB6B5aw', {
    maxZoom: 18,
    attribution: 'Map data &copy; OpenStreetMap contributors, ' +
        'CC-BY-SA, ' +
        'Imagery © Mapbox',
    id: 'mapbox.streets'
}).addTo(mymap);

var StratStyle = {
    "color": "#d3d3d3",
    "weight": 5,
    "opacity": 0.65
};

var geojson = new L.GeoJSON.AJAX("StratfordBoundary", {style:StratStyle});
geojson.addTo(map);

here is the entire code


texts: Shapefile for India which is up-to-date to display regional borders

texts:
I need a shapefile for India which is up-to-date to display the regional borders. By up-to-date (in 2014), I mean that it includes Telangana as a state. Many shapefiles that I downloaded do not have Telangana as a state as it was partitioned a month ago, it is displayed as part of Andhrapradesh.

So I need the latest shapefile.


texts: How can I merge/collapse nearby and parallel road lines (eg. a dual-carriageway) into one single line?

texts:
I have a shapefile from OSM which contains all the roads within a city. Larger roads (like dual-carriageways) seem to have 2 parallel lines. I'd like to collapse / combine those into 1 line so I can export and render them nicely in an external 3D application (2 overlapping roads will be complicated and will look weird in 3D).

How can I achieve this with either QGIS or PostGIS? I am not worried about a slight loss of accuracy (within a few metres) and would like the resulting, single line to be between (ideally the centre of) the current, parallel lines.

Thank you.

(here is an example of the dual road lines that I want to combine)


texts: gdal image quality problem when using gdalbuildvrt?

texts:
I have ~1000 1GB images that I want to make into a COG w/ JPEG compression.  Any time I try to use gdalbuildvrt to create a vrt of my entire directory of images then I get big quality loss.
Can someone show me what I'm doing wrong?
Left image is a sample of how I'd like to make my entire directory of images into a single COG.  This VRT has only 2 images, N00-E000_ur and N00-E000_ul.

gdalbuildvrt test.vrt *.tif

gdal_translate test.vrt testjpg90vrt.tif -of COG -co COMPRESS=JPEG -co QUALITY=90

Right image is just a comparison of making a COG from just 1 of the images:

gdal_translate N00-E000_ur.tif N00-E000_ur_cog_jpeg.tif -of COG -co COMPRESS=JPEG -co QUALITY=90


texts: compare three rasters and three conditions with overlay in r

texts:
I'm trying to run an overlay function with rasters where I want to meet all of 3 different conditions at each cell and produce a single raster as output.

Running ifelse with the & operator seems to look at the conditions in a linear fashion from left to right - If the first two conditions are met then it will produce the if condition as output, regardless of the third condition.  I think that to meet all three conditions I would need something equivalent to &&, but this can't be used here because it is not vectorized.

You can see this with this example below.  Since the third raster has NA values and the result raster has fewer NA values, it is clear that it is not evaluating all three arguments.

Hoping someone can point me in the right direction.

library(raster)
fn <- system.file("external/test.grd", package="raster")
  s <- stack(fn, fn,fn)
  s[[1]] <- round(runif(ncell(s), 1, 2))
  s[[2]] <- round(runif(ncell(s), 1, 2))
  s[[3]] <- round(runif(ncell(s), 1, 2))
  #convert some values in s[[3]] to NA
  s[[3]][s[[3]] == 1]<- NA

  result.rast <- overlay(s[[1]], s[[2]], s[[3]], fun =
      function(x,y,z) {
            ifelse( x == 2 & y == 1 & z ==2, 1, 0)
      } )

I copied some of the building code from this post Compare two rasters with raster overlay and replace cell value with variable


texts:
Which algorithm, that is based on open Street map takes the current GPS position, a duration and an average speed or way of transport(by foot, by bike or by car) and outputs possible routes, that start and end at the current GPS position and fulfill all criteria? That means, you start from a GPS position, you move around till the time has elapsed and when the time has elapsed, you should be at the same place where you started.

Usecase:

You are in a new country. You just arrived at your hotel. In one hour you are going to meet a friend of yours. The two of you agreed to meet at your hotel. You decide to get to know the close surrounding a little bit by walking. At the same time you want to be back in time. You would like to not see the same road again when you come back. So you grant GPS access in your mobile device, tell the amount of time that you have and that you don't want to see the same road again if possible. Then the program should output a series of routes that you can walk that fulfill those criteria.


texts: Write shp into python and fill in the attribute table

texts:
I am struggling to create a shp and I want it to writes records on the attribute table. I can create the shp but I can not fill in with the correct records.
I import an Excel file, that contains the fields "code", "x_Start", "y_Start", "x_Finish", "y_Finish" and I want the field 'code' on the shp table.

Here is the code:

import pandas as pd
import geopandas as gpd
import requests
from shapely.geometry import Point

#import excel
import requests
df = pd.read_excel("NI.xls")
df.head()

#then I create trajectories through an API and I get a list.

a = list(zip(*[list(df['x_S']), list(df['y_S']), list(df['x_F']), list(df['y_F'])]))
#print (a)
b = [bb for bb in [create_trajectory(coords) for coords in a] if not bb is None]
dg = pd.DataFrame(b)
dg

#I can create the shp from dg but I can't know how to fill in the attribute with th df. Ν.Β. that I do not get the same number routes as I give on excel file.

import shapefile

# Create a polyline shapefile writer
with shapefile.Writer("shapefiles/test/routeNI6") as w:
    #w = shapefile.Writer(shapefile.POLYLINE)

    # Add our fields
    w.field("NAME", "C", "40")

    #lines = []
    for i, row in df.iterrows():
        w.record(df.code)

        w.line([list(row.dropna())])

I want to add only the 'code'  records. and here are the errors that I get

TypeError                                 Traceback (most recent call last)
<ipython-input-115-04aa2ea41a9f> in <module>
     21
---> 22         w.line([list(row.dropna())])

~\Anaconda3\lib\site-packages\shapefile.py in line(self, lines)
   1670         shapeType = POLYLINE
-> 1671         self._shapeparts(parts=lines, shapeType=shapeType)
   1672

~\Anaconda3\lib\site-packages\shapefile.py in _shapeparts(self, parts, shapeType)
   1756                 if not isinstance(point, list):
-> 1757                     point = list(point)
   1758                 polyShape.points.append(point)

TypeError: 'float' object is not iterable

During handling of the above exception, another exception occurred:

ShapefileException                        Traceback (most recent call last)
<ipython-input-115-04aa2ea41a9f> in <module>
     20         #lines.append(list(row))
     21
---> 22         w.line([list(row.dropna())])

~\Anaconda3\lib\site-packages\shapefile.py in __exit__(self, exc_type, exc_val, exc_tb)
   1099         Exit phase of context manager, finish writing and close the files.
   1100         """
-> 1101         self.close()
   1102
   1103     def __del__(self):

~\Anaconda3\lib\site-packages\shapefile.py in close(self)
   1120                 raise ShapefileException("When saving both the dbf and shp file, "
   1121                                          "the number of records (%s) must correspond "
-> 1122                                          "with the number of shapes (%s)" % (self.recNum, self.shpNum))
   1123         # Fill in the blank headers
   1124         if self.shp and shp_open:

ShapefileException: When saving both the dbf and shp file, the number of records (1) must correspond with the number of shapes (0)

I have added the libraries that I use. Actually, I have some couple of thousand polylines. This code is a part of the whole code.

I import the excel file as df and then I create trajectories as dg. So I want to create a shp with features from the dg and to fill in the attribute from df.
I hope it helps.
I am sorry for the inconvenience but I am quite fresh.


texts: Where can I get (USA) county neighbour relationship data?

texts:
I need a graph where each node is a US county and each edge represents a border shared between counties.  I do not particularly care about the absolute position or shape of each county (though that would be a plus).

Where can I find that information (for free)? Army Corp of Engineers? US Geo Survey?

Ideally, I could just get a csv list something like

FL-Polk, FL-Lake
FL-Polk, FL-Orange
FL-Polk, FL-Osceola
etc for every county, including "kitty-corner" borders, borders across water, or across state borders

I'm a programmer so expect I can handle any standard exchange format.  (That could just be hubris though). (I don't have any cool GIS apps.  Just humble scripting in perl and python.)


texts: geotools parsing kml with multigeometry only gives me first geometry

texts:
I'm trying to parse the states.kml file (http://docs.geotools.org/latest/userguide/_downloads/states.kml) on the geotools website using the following code:

Parser parser = new Parser(new KMLConfiguration());
SimpleFeature f = (SimpleFeature) parser.parse( is );
@SuppressWarnings("unchecked")
Collection<SimpleFeature> placemarks = (Collection<SimpleFeature>) f.getAttribute("Feature");

for (SimpleFeature fe : placemarks) {
    Geometry g = (Geometry) fe.getAttribute("Geometry");
    for (int i = 0; i < g.getNumGeometries(); i++) {
        System.out.println(g.getGeometryN(i));
    }
}

Each of the geometry fields for each placemark in that file is a MultiGeometry made up of a Point and a Polygon. When I run through my for loop and print out though, only the points are printed. Does anyone know how I can access the Polygon also? Is it accessible, or did it get dropped in parsing?

Thanks


texts: Closest projection to GeographicLib's Local Cartesian

texts:
We are using the excellent GeographicLib https://sourceforge.net/projects/geographiclib/ in our application, a professional ship simulator. Our 3D scene is rendered using a Local Cartesian projection provided by this library.
But we need to use some GIS tools, including QGIS, to edit the scenario, but we don't known if there is an equivalent projection based on GDAL/PROJ.4. The best results was obtained with +proj=tmerc +ellips=wsg84 +k=0.99997 but It doesn't matches exactly. I tried ortho too, but it looks far than tmerc.

Does anyone know which GDAL/PROJ.4  is equivalent to that Local Cartesian?


texts: Converting multiple polygons into rasters using ArcGIS Desktop

texts:
I want to convert a couple of thousands of polygons (fish species data, each fish species shown as one polygon) into raster data (size of 10*10km), the conversion tool (polygon to raster) in ArcGIS only allows me to do each polygon individually which is impractical given the data size. I'm considering using R to rasterize the data or using iteration in Modelbuilder (using ArcGIS); but I'm not very familiar with both methods, I was wondering if there's any other way doing this efficiently or elaborate these two methods.
My goal is to show the number of fish species per grid (or species richness per grid; I will make each raster size to be 10 by 10km for all the polygons). I've already added a new field called presence (with value of 1) to each species polygon. My next step (I think) is to rasterize the polygon data using field value of "presence" and set raster size to 10000; in the end, I can tell how many species per pixel cell (size of 10 by 10km) by looking at the "presence" value.
Correct me if I'm wrong, thanks in advance.
edit:


texts: How to change the default gray color of the QgsComposition background?

texts:
I'm building a map, with elements like QgsComposerMap of course, but also QgsComposerLabel, QgsComposerLegend ,etc...
Each of these elements has a boolean "setBackgroundEnabled" property, to make them white backgrounded or not.
But when aggregating these white backgrounded elements in the QgsComposition, we can see gray bands between them, caused by the deepest gray background. And it's printed when saved as image !

As for canvas, ( i mean canvas.setCanvasColor(Qt.white)), is there a way to set the color of the overall QgsComposition background (the deepest one) ?

Is there something like :

comp = QgsComposition(mapRenderer)
comp.setCompositionColor(Qt.white)


texts: DC-GIS WMS server. Metadata shows polygon layer, why am I getting a raster

texts:
This page gives the description of the layers available from the DC_GIS mapserver. When I look up the metadata for various layers, say record lots, the description says
           Geometry Type: esriGeometryPolygon
but when I fetch the layer using either qgis or ArcMap, I get a raster layer. I've tried this for several different layers which describe themselves a polygons, but always get a raster.

In ArcMap I download the .lyr file. In qgis, I copy the WMS URL, open it in the New WMS servers dialog and select the layer to add.

Is it me or is it the metadata?


texts: Is it possible to download all the contents of OpenRailwayMap.com as a shapefile?

texts: OpenRailMap.com is a free website. You can use its contents as a layer in Google Maps, Open Layers and leaflet. Is it possible to download its contents as a shapefile?


texts: Limit to number of overviews in QGIS print composer?

texts:
I am creating a general map that has seven (7) detailed map frames in QGIS 2.18.10. I want to show the overviews of all these frames on my general map. Unfortunately, not all map frames show in the combo box.

Is there a limit to the number of overviews you could show?