texts:
I've written a PyQGIS script based on this answer to the following question:

Writing numpy array to raster file

My data shows up at the correct, expected location but the raster cells, which are 25m x 25m squares, are still X/Y aligned when I would like them to be rotated as well (aligned with the red polygon limit, not the X/Y axis).

Is there a way to do this with QGIS / PyQGIS / GDAL ?


texts: Shapefile disappearing in exported map image from ArcGIS Pro

texts:
I'm exporting a map image from ArcGIS Pro, which is made of three overlayed shapefiles.
Everything looks dandy in the layout, but after exporting, one of the shapefiles is not on the map. The legend contains the symbology for the layer, but it is missing in the map frame.
What could be the problem?
I tried exporting in different formats, but it doesn't make any difference.
I'm using one Map and one Layout with one map frame.


texts: Listing Unique Pixel Values From Raster in Google Earth Engine

texts:
I'm having an issue with what I hope is a rather simple task. I have a raster of worldwide extent and I simply want a list of the unique pixel values in a given band.
Example:
//load in fire data from the last two years as dummy data
var burn = ee.ImageCollection('ESA/CCI/FireCCI/5_1')
                  .filterDate('2018-11-06', '2020-11-06');

How do I grab a list of unique values from this?
My real dataset is pretty sparse on a worldwide level, so ~300m is where I'd like to set my scale.
Furthermore, can I use this list to filter out entries in a FeatureCollection that have the unique values I find?
I've tried a variety of methods and a number of S/O threads... so far either reduceRegions has produced a memory error, or a toDictionary/frequencyHistogram method has turned up with empty outputs.


texts: R: Display single band raster in true color

texts:
I have a static image of a Bing maps tile (Mercator EPSG: 3857) with satellite overlay in PNG format. I read the PNG as a raster and georeferenced it in R by setting the corner coordinates using the raster package.

Here is the original PNG.

r <- raster("mymap.png")
# add coordinate system
crs(r) <- "+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"
# plot raster
plot(r)
# try some different colors, not true color though
colortable(r) <- rainbow(256)
plot(r)

However, when I try to plot the raster, I get some kooky colors. I can change the colortable, but it doesn't match the true color of the PNG. What is happening to the "color values" when I read a PNG with the raster package?

Here is the raster:

Formal class 'RasterLayer' [package "raster"] with 12 slots
  ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
  .. .. ..@ name        : chr "C:\\Users\\me\\Documents\\traffic\\bing\\images\\mymap.png"
  .. .. ..@ datanotation: chr "INT1U"
  .. .. ..@ byteorder   : chr "little"
  .. .. ..@ nodatavalue : num -Inf
  .. .. ..@ NAchanged   : logi FALSE
  .. .. ..@ nbands      : int 3
  .. .. ..@ bandorder   : chr "BIL"
  .. .. ..@ offset      : int 0
  .. .. ..@ toptobottom : logi TRUE
  .. .. ..@ blockrows   : int 1
  .. .. ..@ blockcols   : int 2560
  .. .. ..@ driver      : chr "gdal"
  .. .. ..@ open        : logi FALSE
  ..@ data    :Formal class '.SingleLayerData' [package "raster"] with 13 slots
  .. .. ..@ values    : logi(0)
  .. .. ..@ offset    : num 0
  .. .. ..@ gain      : num 1
  .. .. ..@ inmemory  : logi FALSE
  .. .. ..@ fromdisk  : logi TRUE
  .. .. ..@ isfactor  : logi FALSE
  .. .. ..@ attributes: list()
  .. .. ..@ haveminmax: logi TRUE
  .. .. ..@ min       : num 0
  .. .. ..@ max       : num 255
  .. .. ..@ band      : int 1
  .. .. ..@ unit      : chr ""
  .. .. ..@ names     : chr "mymap"
  ..@ legend  :Formal class '.RasterLegend' [package "raster"] with 5 slots
  .. .. ..@ type      : chr(0)
  .. .. ..@ values    : logi(0)
  .. .. ..@ color     : logi(0)
  .. .. ..@ names     : logi(0)
  .. .. ..@ colortable: logi(0)
  ..@ title   : chr(0)
  ..@ extent  :Formal class 'Extent' [package "raster"] with 4 slots
  .. .. ..@ xmin: num 0
  .. .. ..@ xmax: num 2560
  .. .. ..@ ymin: num 0
  .. .. ..@ ymax: num 1280
  ..@ rotated : logi FALSE
  ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
  .. .. ..@ geotrans: num(0)
  .. .. ..@ transfun:function ()
  ..@ ncols   : int 2560
  ..@ nrows   : int 1280
  ..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr "+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +no_defs"
  ..@ history : list()
  ..@ z       : list()

I'm not entirely sure why it says there are 3 bands when this is just a single band PNG image.


texts: I am working in Qgis. I have a vector layer called Kimball Building with about 5,000 records in it. I have multiple fields including one called PermitType. Most of the values in this field are NULL. Some records have another value in the PermitType field. I would like to delete all of the records that have NULL in the PermitType field. I am thinking this could be done in Python, but not sure where to begin. I can delete items individually, but deleting 4,000 records one at a time is not a practical solution.


texts: Using label expression to label highway with last group of numbers in field

texts:
I'm using Tiger road files and I want to label the highways with their shield number. Tiger has an inconsistent naming convention for road features; I am specifically concerned with the highways. I have defined 3 seperate label classes: interstate, US, and state route. For each I applied an SQL query of their corresponding route types.

The part I am struggling with is creating a label expression for each that would include only the last SET of numbers. I say "set" because there are varying numbers of characters annotating the highway number.

Example:
US Hwy 24,
US Hwy 6,
State Hwy 78 N

Is there a way to include only the highway number for the label?


texts: Loading NAD83 and WGS84 layers together in QGIS?

texts:
I have downloaded parcel files for several counties in Washington State, and am attempting to overlay them on the OpenStreetMap Standard base layer (which is in WGS84 Pseudo Mercator). The majority of the county layers were in WGS84 and automatically loaded into the correct extent/location, but some (such as King County) were in NAD83 (ESPG: 4269) is placed over West Africa in the OSM map.

What do I need to do to get this King County layer to load in the correct place on OpenLayers map? Or alternatively, how do I convert the file to the correct CRS using ogr2ogr?

Here are my project CRS settings:

If I set the CRS for the layer to the correct settings NAD83 (EPSG:4269), the layer is not visible. (but doing zoom to layer still takes me to west Africa). If I change the CRS in layer properties to EPSG:3857, it becomes visible but is still in the wrong location.

I also tried saving the King County NAD83 layer as a new file with CRS set to WGS84 PseudoMercator, and using the extents laid out in the metadata files:

... but when I loaded it, I can't see the layer or make it visible no matter what CRS I choose for it. So then I gave up trying to do it in QGIS, and attempted to use some command line tools to try to convert/reproject the file.

Here is the output of ogrinfo for the King County NAD83 shapefile:

ogrinfo -so -al parcel.shp
INFO: Open of `parcel.shp'
      using driver `ESRI Shapefile' successful.

Layer name: parcel
Metadata:
  DBF_DATE_LAST_UPDATE=2018-04-06
Geometry: Polygon
Feature Count: 617501
Extent: (1218256.223500, 30932.680000) - (1577865.300000, 287704.244000)
Layer SRS WKT:
PROJCS["NAD_1983_HARN_StatePlane_Washington_North_FIPS_4601_Feet",
    GEOGCS["GCS_North_American_1983_HARN",
        DATUM["NAD83_High_Accuracy_Reference_Network",
            SPHEROID["GRS_1980",6378137.0,298.257222101]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.0174532925199433],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["False_Easting",1640416.666666667],
    PARAMETER["False_Northing",0.0],
    PARAMETER["Central_Meridian",-120.8333333333333],
    PARAMETER["Standard_Parallel_1",47.5],
    PARAMETER["Standard_Parallel_2",48.73333333333333],
    PARAMETER["Latitude_Of_Origin",47.0],
    UNIT["Foot_US",0.3048006096012192]]
MAJOR: String (6.0)
MINOR: String (4.0)
PIN: String (10.0)
Shape_area: Real (19.11)
Shape_len: Real (19.11)

I tried using ogr2ogr to convert it to WGS84, but got an error message:

ogr2ogr -f "ESRI Shapefile" -t_srs EPSG:3857 -s_srs EPSG:4269 'King County Parcels (WGS84, EPSG 3857).shp' parcel.shp
ERROR 1: Failed to reproject feature 0 (geometry probably out of source or destination SRS).
ERROR 1: Terminating translation prematurely after failed
translation of layer parcel (use -skipfailures to skip errors)

... I'm at a complete loss for what to do. I've looked through tons of tutorials and GIS Stackexchange posts on projections in QGIS, and seem to be missing something, but can't figure out what.

What am I doing wrong here, and how can I make this layer work with OpenStreetMap, and the other WGS84 layers?


texts: CartoDB: how to let users update data from the map infowindow

texts:
In cartoDB, I have a map showing public lands in a German city. I want to integrate a form into the cartoDB infowindow to let the user update the data table.

How can I pass the values from the form to the SQL statement? Should I call a php script or is it better to use javascript?


texts:
how can I create a new pen style in qgis 2.16?

I want to create a pen style that looks like that - - . . - - . . (dash dash dot dot). This is not defined yet.

I tried to use two simple lines and two marker lines (1st option) or four marker lines (2nd option) or use custom dash pattern (3rd option) but all does not work fine. The problems are:

1st Option:
This would work fine, if it is possible to use offest in x direction (offset is only defined in y direction) for simple lines. Is this possible?

2nd Option:
This would work fine too, but there is no way to use an user defined rectangle as dash.

3rd Option I tried was use custom dash pattern. Is it possible to define dash (10px) space (5px) dash (10px) space (20px)?

Is there a way to define an dash dash dot dot line?


texts: PostgreSQL : question about SELECT .. WHERE .. = max()

texts:
If you have, for example;

SELECT * FROM data
WHERE line = min( (SELECT length FROM lines) )
      AND number < 500;

(It's just an example, not a real code)

How does PostgreSQL react when the result from the MIN has a number > 500 ? Is he then going after the second MIN or what happens ?


texts: How to set up a background map in WGS84 projection for interactive map?

texts:
I've started developing my own interactive map and I have some things which are not clear to me.

I heard that OSM, Google maps and all similar maps use Google Mercator projection. But if I use that map as my background map can I transform it to WGS84 because my other data is in WGS84 projection. Is it better to do it immediately in QGis or is it possible to change it later in some other step of creating web gis? (is it possible to do it using postgis, geoserver or openlayers?).

Is it even possible to transform it to WGS84 or it has to stay at Google Mercator?

And one more thing. What is actually the best way to add background map to my future web GIS application? Should I do it immediately in QGis (if so, is it possible to import it in GeoServer later because I can't have .shp file of it because its not shp file, it's just satellite or other type of data provided with QGis plugins.. and as far as I know, I can import only shp files in geoserver, or tables from postgis) .. or can I do it some other way?


texts: How can I index USGS Landsat Collection 2 Level 2 datasets into the Open Data Cube?

texts:
I'm trying to ingest data using Landsat 8 Collection 2 Level 2 data into the Open Data Cube, but I couldn't find any examples for data preparation.
What is the best way to index USGS Landsat Collection 2 datasets into the Open Data Cube?


texts: Using Orthophotos - not multispectral - to generate Object-oriented image classification vectors?

texts:
I'm using  QGIS version 2.14.3-Essen on a Macbook pro. I'm sort of new to this and feeling my way through this.

I have high-resolution Orthophotos officially sourced from local country in Sub-saharan Africa. These are already composite TIFFS and not multi-band.

I want to be able to visualise the different land coverage i.e agriculture, buidlings, forestry etc from these images, and to have this in a vector file format.

I'm not sure how to do this.

Exploring the web I came across Semi-Automatic Classification Plugin for QGIS. However I am not getting very far and I'm wondering is it because they are Flattened Tiffs and not multispectral? I keep getting errors when I try to draw the ROI's.


texts: ArcGIS Portal use server-side obtained authorization in the client / JavaScript

texts:
I have an ASP.NET MVC web application that is authenticating users server-side against ArcGIS Portal successfully using Owin.Security.Providers.ArcGISPortal.

I would like to use the granted authorization client-side in the browser without requiring the user to login a second time via ArcGIS's JavaScript API. Does anyone have any pointers regarding how to obtain a relevant token and pass this through client-side to access portal secured web maps, etc.?

There is a context.AccessToken that I obtain via ArcGISPortalAuthenticationHandler and make available as a claim in order to share the access token.

Any ideas? If you know of a working example for ArcGIS Online then this could also help me as ArcGIS Online and ArcGIS Portal seem to be identical in their mechanics.

Further investigation...

Following the linked example found through the documentation for esri.IdentityManager.registerToken(), I've tried the following method to share the access token:

var credentialsJSON = {
    serverInfos: [{
        server: "https://[HOST]",
        tokenServiceUrl: "https://[HOST]/arcgis/tokens/",
        adminTokenServiceUrl: "https://[HOST]/arcgis/admin/generateToken",
        shortLivedTokenValidity: 1800,
        currentVersion: 10.5,
        hasServer: true
    }],
    oAuthInfos: [],
    credentials: [{
        userId: user.userId,
        server: "https://[HOST]/arcgis",
        token: user.userAccessToken,
        expires: user.userAccessTokenExpiry,
        validity: 1800,
        isAdmin: false,
        ssl: false,
        creationTime: user.userAccessTokenIssued,
        scope: "server"
    }]
};
esriId.initialize(credentialsJSON);

But I'm still not getting access to the webmap. I see a html login modal over the empty map div with the title "Please sign in to access the item on https://[HOST]/arcgis (b11824af61df463586dad40d1df7abbd)".

In the console log I see the following message logged:

dojo.io.script error Error: You do not have permissions to access this resource or perform this operation.
    at Object.g.load (init.js:984)
    at init.js:87
    at c (init.js:103)
    at d (init.js:103)
    at a.Deferred.resolve.callback (init.js:105)
    at c (init.js:104)
    at d (init.js:103)
    at a.Deferred.resolve.callback (init.js:105)
    at init.js:999
    at n (init.js:107)

and for the network request to

https://[HOST]/arcgis/sharing/rest/content/items/b11824af61df463586dad40d1df7abbd?f=json&callback=dojo.io.script.jsonp_dojoIoScript1._jsonpCallback

I see the following response:

dojo.io.script.jsonp_dojoIoScript1._jsonpCallback({"error":{"code":403,"messageCode":"GWM_0003","message":"You do not have permissions to access this resource or perform this operation.","details":[]}});


texts: How to add GeoServer vector tiles in ArcGIS using ArcGIS JS API

texts:
I am trying to display vector tile layers from GeoServer in ArcGIS using JavaScript API. I have referred this.here they have used tegola service but i want to add "application/x-protobuf;type=mapbox-vector" which is published in GeoServer.

In OpenLayers I am able to display using Mapbox vector tiles service is "http://localhost:8080/geoserver/gwc/service/tms/1.0.0/Cadastral:CADASTRAL@EPSG:900913@pbf/{z}/{x}/{y}.pbf"

Same thing I want to add in ArcGIS. How to add it in ArcGIS using JavaScript API?


texts: Calculating Center (GPS Point) and Radius (km) of n GPS points

texts:
I'm new to GPS and position mapping but i have a background in programing and some of the math related to it.

I've been tasked with writing a small golang program that takes n GPS points (DMS) and returns a GPS point and Radius that makes a circle that contains all n GPS points.

if n <= 1, Not Valid.

if n = 2, The GPS point returned will be the midway point and the radius will the distance divided by 2.

if n >= 3, The GPS point and Radius that makes a circle that contains all n GPS points.

Calculate the Centroid of a polygon of n GPS points, returns center GPS point in (DMS) format.
Calculate the distance from the above centroid to each coordinate from the n set using Pythagoras. The biggest value will be the radius of the circle.

This is as far as my thought process goes. Is there any GPS specific math I need to find points in (DMS)? Are there golang packages that could help with this type of math? Can someone give me an example of this with n >= 3 How can I drop pins and a radius to check the results, google maps?


texts: Tiled Map Service (TMS) in QGIS not showing highest resolution?

texts:
I've added a Tiled Map Service (TMS) in QGIS, but when i try and look at the highest zoom levels, the highest resolution I can see is level 19. When I add the service as a WMS, I can get the full resolution, level 21, but I can no longer call a specific image date from the server.

Does any one know how to call a higher resolution using a TMS in QGIS?

Is this a minidriver issue?

I don't have any code to post, but I've been fighting with this for days and have just run out of ideas on how to fix this.


texts: Change language from Swedish to English QGIS pisa

texts:
Can I change language from Swedish to English in my QGIS Pisa?

I don't remember even choosing languages when I installed it.


texts:
I have a large number of ECW format raster images and would like to use with GDAL (and ultimately in QGIS) on Mac OSX.

Having seen that there is no official ESRI support for the ECW library on any other platform other than Windows, I managed to get hold of the old version of the libecwj2-3.3 library. There are a few versions knocking around, but none of them will successfully compile.

Anyone had any success compiling the ECW library on OSX (10.6.7) and then compiling GDAL with ECW support.


texts: Leaflet Geodata Marker clustering and intersection

texts:
I'm having a problem with grouping(clustering) markers
I have a list geosjon features that I display as markers in the leaflet. And it looks like this:

These points are saved in my Database and are located on my Geoserver.
And when fetch them from DB I get it like this, then I loop over them and add them as markers.
const points = [
  { lat: 48.371768, lng: 16.723857 },
  { lat: 48.354206, lng: 16.703258 },
  { lat: 48.401474, lng: 16.727156 },
  { lat: 48.283497, lng: 16.545787 },
  { lat: 48.361275, lng: 16.699054 },
  { lat: 48.283497, lng: 16.545787 },
];

The problem that I'm having is that I have to summerize points on a within a specific radius and grouped(cluster) them automatically.
I have to expand points for some specific radius say 100 meters
Like in Turf js buffer option:
https://turfjs.org/
Like here:

And in Points where those circles intersect with each other, add an additional marker with those i already have
I have tried with Leaflet Clustering marker
https://github.com/Leaflet/Leaflet.markercluster
But haven't had any luck, markers are not displayed at all.
Here is my code so far:
var markers = L.markerClusterGroup({
  maxClusterRadius: function (zoom) {
    return 10;
  },
});

points.forEach((point) => {
  const location = { latlng: { lat: point, point: event.lng } };
  var marker = L.marker(location.latlng);
  markers.addLayer(marker);
  this.layers.markersLayer.instance.addLayer(markers);
});

And that works fine but I would like to cluster points, by distance in meters instead of pixels like in maxClusterRadius


texts: Layers are not visible after restarting GeoServer

texts:
I have ArcSDE connection on GeoServer. In order to make GeoServer faster I have restarted it. After restarting my published layers (datas) are not visible. When I publish this layer again, GeoServer says that "you already published this layer"

How can I see my published layers after restarting GeoServer.


texts: Permanently Reordering Subtypes in ArcGIS Pro?

texts:
I have a feature class in a file geodatabase with subtypes (designed so when you select a country, you are presented with every province in that country as a series of domains).

I recently updated the name of a country from Saudia Arabia to Kingdom of Saudia Arabia. Everything is functioning normally, however, my subtypes are now out of order.

I was unable to find a subtype manager in ArcCatalog or ArcMap. I used the one in ArcGIS Pro 2.0 and reordered the subtypes, and clicked Save. It appears that this adjustment was only temporary.

Is there a tool that sorts subtypes, somewhere?

Ideally, something like the Sort Coded Value Domain tool.


texts: Defining and changing unknown CRS into WGS84

texts:
I came across a shapefile (no information available about the original CRS) with some unknown Reference System (City of Stuttgart in Germany). Example of coordinates looks like 3465878, 5482952. I assume it is lat/lon.

How can I change the CRS to the decimal system of WGS84 with QGIS 3?

Is it possible to overlay it to another shapefile with the borders and match the CRS?


texts: ListFeatureClasses brings back more than expected

texts:
Why, when I issue the function:

fcList = arcpy.ListFeatureClasses(feature_type='point')

does it bring back a list of grid/rasters as well as all my point featureclasses from ArcSDE.  All I want is my point featureclasses, this doesn't make any sense!  Anyone had this happen to them? Plus it is incredibly slow!

Any help much appreciated :)


texts: How to load Interlis1 files directly into QGIS 2.0.1 Dufour 64Bit (W7Pro)

texts:
With earlier version of QGIS (1.8) it was possible to load data in Interlis 1 format directly.

If I try now it says something like unknown format. Is this because I'm using the 64Bit version of the installer?

A possible workaround could be using the examples at http://giswiki.hsr.ch/HowTo_OGR2OGR "interlis to shp".

It would be very nice to be able to open Interlis Files directly again.


texts: Move points along a line, 100m left or right

texts:
In ArcMap 10, I have several hundred points taken along a zig-zag transect (line) but I need to move these points 100m to the left or right, perpendicular to the line. The points are stored as a single shape file with a 'SIDE' field, which is either 'left' or 'right' indicating which side of the line to move the point.

I have tried linear referencing but got totally lost. Is that the tool to use or is there another simple method of achieving this?


texts: How to rotate a geometry created by the geometry builder in QGIS

texts:
This topic goes into making polygons based on point locations;
Using QGIS Geometry Generator to get rectangle from point?
Is there a way to apply a rotation to the geometry generated? I have a field that contains a 0-360 degree rotation value and would like to rotate my built geometries based on that value.


texts: qgis raster style 'no data value' checkbox is ignored - what's a good workaround?

texts:
See bug report filed here:

http://hub.qgis.org/issues/10236

Note the initial bug report got 'fixed' but that fix introduced this new error, so the issue was reopened a few days ago.

Anyway, it looks like the team might have bigger fish to fry at the moment - thanks to anyone who does find time to work on it - but in the meantime, can someone suggest a good workaround? I was hoping to just edit the geotiff in gimp and do a color map from #000000 to #000001 or something, so that the 'almost-black' areas are no longer picked up as transparent in qgis - but, the file is too big (it's a labels layer from a USTopo) and gimp crashes in addition to complaining about some data fields on reading the file which kind of indicates the output may be dodgy anyway.

Is there a gdal tool to map colors in a geotiff?  Or, some other script or program that anyone could recommend?

Thanks in advance.


texts:
I have trouble to use function dissolve to dissolve polygons by field value in PyQGIS3.
In PyQGIS2 this line was doing the job:

QgsGeometryAnalyzer().dissolve(ki, "name.shp", onlySelectedFeatures=False,uniqueIdField=-1, p=None)

I read on this site that processing alg should be use in PyQGIS3. But I can't find here disolve. When I run processing.algorithmHelp("native:dissolve") I can't get the list of attributes, I just get some basic information. So I don't know how to use it.


texts: Can you make a map of a fictional country and cities in QGIS?

texts:
I am currently working on a project were I am making a fictional country with cities in it. I am wondering if I can do that in QGIS, and if so, how?

I have no experience with any GIS program at all. I just saw QGIS on some forum, and it seemed interesting. The fictional country is based on the Mascarene Plateau in the Indian Ocean. I  want to give the map a more professional look to it and make it look like an actual map, like in an atlas. I have a pretty clear idea on the vegetation and elevation of the landmass, and for the cities, I am still drawing the layout of them.


texts: Is there a map of map services somewhere?

texts:
I'd like to provide a map based search of map services in a web app.

Ideally, I'd like a web service that displays markers at different zoom levels: local, municipal, county, regional, state or national level.  When I click on a marker, I'd like to see a list of urls representing the map services published by that governmental entity.

For example, for a typical state I would expect to see a marker for each county that has public web mapping services.  The marker would appear on the county seat, and would contain a url to that counties' map services.

My primary focus is Texas, USA, but any examples of map service map services would be useful.


texts: Is it possible to postprocess Trimble .ssf files without Trimble Pathfinder Office?

texts:
Is it possible to postprocess Trimble .ssf files without Trimble Pathfinder Office (differential correction - our  sporadic use does not defend realtime correction subscription) ?
We only need the correction utility but the pathfinder package contains a lot more and is too expensive just for this.
I have not yet found open source software which supports ssf.


texts: PyQGIS plugin doesn't work layer list combo box

texts:
I want to create a simple QGIS plugin with a simple combo box with list from my layers

in _ _init_ _ :

    legendInterface = self.iface.legendInterface()
    listLayerName = [i.name() for i in legendInterface.layers() if i.type() == QgsMapLayer.VectorLayer]
    # add all these layer names to combo box comboLayerList
    self.ui.comboBox_2.addItems(listLayerName)

but if I start plugin in QGIS then the combo box not working is empty:

if I reload my plugin use reload using Reload plugin :

then the combo box is working. If I close QGIS and start again then combo box is not working again. If I again reload plugin then works.

Where is my mistake and my combo box not working?
how to refresh the plugin on start to work combo box?


texts: I am very new to Geographic Information Systems. I have a very basic inquery; How do I get hold of the newest ortophoto of South Africa to use as a background map for my GPS data? I'm using QGIS but I guess it doesn't matter which platform/system you use..?


texts: Rendering artifacts with linear interpolation in gdal_grid

texts:
I have problems with strange artifacts occuring when rendering GeoTIFFs using gdal_grid with the linear interpolation algorithm.

Input files

I have a samples.csv file with some geo data samples:

latitude,longitude,value
53.8035355,17.020874,19.43
53.803688,17.0284729,19.36
53.8038445,17.0360718,19.38
53.8040047,17.0436707,19.38
53.8041649,17.0513,19.45
53.8043213,17.0588989,19.47
53.8044815,17.0664978,19.55
...

For this data, I created a simple samples.vrt file:

<OGRVRTDataSource>
  <OGRVRTLayer name="samples">
    <SrcDataSource>samples.csv</SrcDataSource>
    <GeometryType>wkbPoint</GeometryType>
    <LayerSRS>WGS84</LayerSRS>
    <GeometryField encoding="PointFromColumns" x="longitude" y="latitude" z="value" />
  </OGRVRTLayer>

I also have a colormap file for the data (for use with gdaldem):

-9999.0 0:0:0:0
-28 245:245:245:255
-26 235:235:235:255
-24 225:225:225:255
-22 190:190:190:255
-20 130:130:130:255
...

I uploaded the input files here.

Procedure

First, I create a GeoTIFF with gdal_grid:

gdal_grid \
  -ot Float64 \
  -of GTiff \
  -txe 16.870876 19.691624 \
  -tye 54.980014 53.746311 \
  -outsize 682 512 \
  -a linear:radius=0.0:nodata=-9999.0 \
  samples.vrt samples.tiff

Then I colorize the tiff with gdaldem:

gdaldem \
  color-relief \
  samples.tiff colormap samples-colored.tiff

Result

If I use the linear interpolation algorithm, I get the following result:

(Notice the red bars at the top.)

If I switch to the nearest algorithm, e.g.

-a nearest:radius1=0.01:radius2=0.01:nodata=-9999.0

the output image looks fine:

It is also fine if I use the linear algorithm, but flip the image by specifying -tye in the reverse order, i.e. -tye 53.746311 54.980014 - so maybe the problem is somehow related to the negative pixel dimensions.

In my tests, I used GDAL 2.3.0.


texts: Rectify geotiff based on bounding box of another geotiff

texts:
I have 2 geotiffs:  "One.tif", and "Two.tif". I want to use gdalwarp to rectify One.tif based on the bounding box of Two.tif.

First I use gdalinfo to get the minx miny maxx maxy of Two.tif for input into the call.

I set up the call:

gdalwarp -t_srs "EPSG:32615" -multi -te 497025.024 4925030.200 497811.456 4925947.904 -ts 256 256 -overwrite One.tif Output.tif

The result, Output.tif, has the correct SRS and coordinates, resolution, etc., but there is no data in the band. It is entirely composed of 0's. When I run gdalinfo on One.tif, it gives me proper statistics for the single Float64 band within that geotiff.

Is it a no_data issue? Do I need more info for the gdalwarp call?


texts: Changing min max of Stretched Raster with ArcPy?

texts:
I have two rasters r1 & r2:

Maximum cell value of r1: 0.9

Maximum cell value of r2: 0.3

Minimum value for both is 0

I am using a color gradient: white - yellow - green - blue - red - black

When I apply the symbology of r1 on r2, the max value of r2 (0.3) is black. The symbology is set in relation to the value range of r2.

r1 = arcpy.mapping.Layer("r1.lyr")
mxd = arcpy.mapping.MapDocument("myMap.mxd")
df = arcpy.mapping.ListDataFrames(mxd)[0]
r2 = arcpy.mapping.ListLayers(mxd, "r2.tif", df)[0]
arcpy.ApplySymbologyFromLayer_management(r2, r1)

However, I want to compare the two rasters. The values of r2 should therefore be symbolized in relation to the values of r1. The max value of r2 would be yellow in comparison to the value range of r1.

Manually, you can do this in the symbology tab by changing the raster max value by hand:

However, I have a large dataset and I would have to do this a couple of hundred times.

Is there a way to automate this with ArcPy?


texts: Transforming polygons centroids in order to make hexbin map in QGIS?

texts:
I want to make a hexbin map of a particular kind. Hexbins can be used to embrace several points in the same area, like in this example :

With this method, a large city can be covered by one hexagon, and several smaller cities can be covered by the same area. But that's not what I'm trying to do.

In fact, I would do my hexbin map with a file of cities in a French district, with adjoining hexagons and one hexagon for one city, without any empty space between hexagons.

The problem is that in this case, every point must be at the same distance from the others. So my issue is that I'm wondering if it's posible to tidy up with QGis a file with polygons' centroids like this into a file with equidistant points:

How can I transform a file with points to obtain a regular pattern without empty spaces, like this example of South China Morning Post ?


texts: I would like to select all the lines within a Polygon and the lines have to be cut at the border to the Polygon. The option spatial queries takes the whole line. However I only want the part of the lines which is within the Polygon.


texts: Error adding data to Geodjango database: parse error... within geometry

texts:
I'm trying to insert some values into a Geodjango database, but getting errors. This is the table:

class Place(models.Model):
    placeid = models.IntegerField(primary_key=True)
    structidx = models.IntegerField(null=True, blank=True)
    oscode = models.PointField(null=True, blank=True)
    waste86 = models.CharField(max_length=120, null=True, blank=True)
    objects = models.GeoManager()

And this is the code I'm using to insert data:

INSERT INTO domes_place (placeid, structidx, oscode, holding)
  VALUES ('10', '1', '1', 'POINT(632000,141000)', 'N');

And this is the error:

HINT:  "POINT(632000," <-- parse error at position 13 within geometry

Please could someone tell me what I'm doing wrong?


texts: Display multiple lines as one line in QGIS

texts:
I am working with the OSM Dataset and can´t find a way to display rails in a proper way. In a big scale the rails are displaying fine but as soon the lines begin to overlap in a smaller scale I don´t like the outcome. I already filtered the lines to only have the main tracks. But especially in train stations where more rails are running next to each other I get too many overlapping lines.

Is there a way to display multiple lines as one line without changing the geometry?

I am using the QGIS version 3.4.0.


texts:
I am able to add label to circlemarker like this

L.circleMarker(points[i],{title: 'unselected'}).bindLabel('Destination').addTo(map);

This adds label which appears on mouse hover on circle marker.

But I want to add static label which will appear regardless of mouse is on that circle marker or not.

I am referring this demo http://leaflet.github.com/Leaflet.label/ for adding static label to circle marker but some how I am not able to do it.
It is working fine with markers but with circle Markers static label is not working.

Also is there any other method to add label on circle marker ?


texts:
I am trying to read in a file from https://nomads.ncep.noaa.gov
I can read all multigrid files like here
but my script fails to read files from other links that are supposed to be Grib files
such as files from here
and fails with the error below
Exception has occurred: exceptions.ValueError
no matches found
  File "/buio/buio/pygrib.pyx", line 610, in pygrib.open.select

or from here
I use this python script to read the files.
I guess the format must be different?


texts:
I would like to know if anyone knows how to incorporate images as attributes in a QGIS project, and that these are displayed even if the folder containing the project and the source files (shapes, images, etc.) change location.

I know how it is done to show the image, but when sending the project, the other person can not visualize the images in the attributes.

For reference guide with the following video:
https://www.youtube.com/watch?v=1hcZPNL9yb8


texts: Label features with letters instead of attribute field text

texts:
I am creating a map of protected areas/parks and want to display a letter label (A, B, C...) inside the polygons of all the parks in the map. The full names of the parks are too long to fit inside their respective polygons and make the map look cluttered. I know I could add a new field to the park layer's attribute table and type in the letters manually line by line, but since I have to repeat this process for several maps, I am seeking a more streamlined way to do this. Ideally the letters would appear in order on the map from left to right (W to E) and top to bottom (N to S) for easy reading, though that's not a requirement. I can't use the first letter of the parks as their label letter, since there are a lot of the names start with the same letter.

My goal is to create a map with the label letters inside the polygons, and then have a reference table with the letters, names, and other attributes on the side, where there's more room to write them out.

I'm working in ArcMap 10.6.1.


texts:
I try to understand the idea behind the different Spatial Reference Identifiers. I understand that the different SRIDs identify different projection views on the earth and that they are important for ensure comparability. This StackOverflow Post also told me:

  All of these [projections] have trade-offs between accurately measuring distance, area or direction -- you can't preserve all three

So I wonder which are the most important SRIDs and for which requirements should I use which ones. Some SRIDs I frequently stumble across are WGS84/EPSG:4326 (Lat/Lng), EPSG:3857(Web Mercator) [former EPSG:900913].


texts: Pixel Value in a Lat and Long location in Envi IDL

texts:
I am new with ENVI IDL,I want to get pixel value of each band in an image in specified lat and long location.I can do a loop among bands in image using this code

ENVI_OPEN_FILE,'E:\MASTER\Payan_Name\DataSer\Output\162\0730\1-1Georefrence_BIP\Subset_BIP.hdr', R_FID=FID, NO_REALIZE=1
ENVI_FILE_QUERY, FID, ns=ns, nl=nl, nb=nb, data_type=data_type, descrip=descrip, bnames=bnames, sensor_type=sensor_type, wavelength_units=wavelength_units, dims=dims

for b=0,nb-1,1 do begin
  indata = ENVI_GET_DATA(FID=fid,POS=b,DIMS=dims)
endfor

But I don't know how to get pixel value at lat and long,I searched and found a library which name is Rdpix is used to get pixel value at mouse location but I want it at specified lat and long location,
Can Some one please help me find the solution or give me some clues to solve my problem?


texts:
There are several discussion about implementing (or not implementing) WFS with .net from scratch.

However, if I do need to, how do I get started?

I know it's a broad question. I need help getting started. Any kind of walk through will be helpful.

I am looking for a step by step walkthrough with details assuming I'm just a GIS person with limited knowledge to a lot of programming concepts.


texts: Is it possible to rearrange the items in qgis layer context menu?

texts: I'm referring to the menu I get when I right-click on a layer in the TOC. I'd like to move Properties to the top. I can't find anything in the customization dialog that seems to apply.


texts: Storing circle in GeoJSON for OpenLayers 3?

texts:
I am loading my layers from external files - I am using geojson format for it, but it's not necessary - I can use any available format.

Is there a way, how to store a circle (ol.geom.Circle) in the geojson file?

Or is there any other file format, which can be directly used in OpenLayers 3 as layer source and which can include a circle definition?


texts: Simple georeferenced 3D data format for building models

texts:
I am looking for a simple georeferenced 3D-data format to store very simple building models. Each building should be represented as a number of poylgons, semantics do not matter. It would be ideal if this format is readable via ArcGIS, but that is not an absolute requirement.

I know there is GML/CityGML but I am looking for something simpler that ideally has either a JSON or XML serialization.


texts: Is it possible to export POINT and LINE in the same DXF from POSTGRESQL using OGR2OGR?

texts:
I have two tables in a pg database. One with annotations/points and one with linestrings.
I would like to export both of them to the same dxf file with ogr2ogr. Is it possible at all somehow?
I have managed to export the lines based on this article.
And the points as well based on this. I also checked the official gdal documentation but there are no examples for my problem.
The ogr2ogr commands:
ogr2ogr -f "DXF" dir\line_layer.dxf "PG:host=*** user=*** dbname=***" -sql "select the_geom from shcema.line_table"
ogr2ogr -f "DXF" dir\annot_layer.dxf "PG:host=*** user=*** dbname=***" -sql "select the_geom, ogr_style from shcema.annot_table"

The ogr_style column contains data like this: "LABEL(f:""Times New Roman"",s:12pt,t:Sz 5)"
So my question is: Is it somehow possible to export from these two different sources to the same dxf file?
(I am using GDAL 3.0.4.)


texts: Cannot perform buffering loop with arcpy

texts:
I have 5 buffers to perform using 5 fields of a shapefile. I try to run it with the following code:

Buffers = ['JF1','JF2','JF3','FF','FB']
env.workspace = "C:\Mydirectory\Myfile.shp"
network = "C:\Mydirectory\Myfile.shp"
out = "C:\Myoutputdirectory\Myfinaldirectory"
for buff in Buffers:
    arcpy.Buffer_analysis(network, out, buff,)

But when I do it I receive the following error:

ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000725: Output Feature Class: Dataset Myfinaldirectory.shp already exists.
Failed to execute (Buffer).

What is the problem?


texts: Cesium library: error with constructor OLCesium

texts:
I receive the error olcs.OLCesium is not a constructor

from line var ol3d = new olcs.OLCesium({map: ol2d}); taken here

my index.html looks like:

<script type="text/javascript" src="libs/v3.16.0-dist/ol-debug.js"></script>
<script type="text/javascript" src="libs/ol3-layerswitcher/src/ol3-layerswitcher.js"></script>
<script type="text/javascript" src="libs/ol3-cesium-1.16/src/ol3cesium.js"></script>

and I have no problems using ol3 library (ol-debug.js). So I think my code works fine with lib paths.

UPDATE1:
from browser inspector I can see that the file ol3cesium.js is regularly downloaded.


texts: Generating true curve elliptical polygons in file geodatabase using ArcPy?

texts:
As background, this question arose from trying to answer a question on How to generate overlapping polygons from lines output from Table To Ellipse tool?

Using the ArcMap 10.1 GUI it is very easy to digitize true curve polygons into a file geodatabase feature class using the Ellipse Construction Tool but ...

Is it possible to write true curve elliptical polygons while reading rows (arcpy.da.SearchCursor) from a table containing a centre point, a major axis, a minor axis and an azimuth for each?

I had hoped that there might be a token available with arcpy.da.InsertCursor to do this, but SHAPE@ seems to be limited by what the Geometry object supports, and that does not appear to include true curves.


texts: Projection: Argument 'crs': Invalid type when displaying ee.Image

texts:
I'm starting with GEE and I try to display in a Jupyter notebook a simple image of a small area to see if I get it right.

I followed the instruction of the documentation and write the following:

img = ee.Image('UMD/hansen/global_forest_change_2015')

ee.Image(url=img.getThumbUrl(({
    'bands':['treecover2000', 'loss', 'gain'],
    'min': 100,
    'max': 1500,
    'region': [[[-107.36, 38.68],[-107.36, 38.32],[-106.22,38.32],[-106.22,38.68]]]
})))

I was expecting a small area to e displayed but instead I get the following error :

  EEException: Projection: Argument 'crs': Invalid type. Expected: String. Actual: Type<Boolean>.

which is strange because I didn't even set the crsargument

Any Ideas ?


texts: Problem with gdal_translate when converting a GeoTiff with spherical datum to a netCDF

texts:
I have a raster file with LCC projection and spherical datum. It actually has the following proj4 parameters:

+proj=lcc +lat_1=45 +lat_2=22 +lat_0=38.7 +lon_0=20 +x_0=0 +y_0=0 +a=6370000 +b=6370000 +units=m +no_defs

I use gdal as following

gdal_translate -unscale -of netCDF "data.tif" "data.nc"

The output netcdf file has the following projection

char lambert_conformal_conic;
  :grid_mapping_name = "lambert_conformal_conic";
  :longitude_of_central_meridian = 20.0; // double
  :false_easting = 0.0; // double
  :false_northing = 0.0; // double
  :latitude_of_projection_origin = 38.7; // double
  :standard_parallel = 45.0, 22.0; // double
  :longitude_of_prime_meridian = 0.0; // double
  :semi_major_axis = 6370000.0; // double
  :inverse_flattening = 0.0; // double
  :spatial_ref = "PROJCS[\"Lambert_Conformal_Conic\",GEOGCS[\"GCS_unnamed ellipse\",DATUM[\"unknown\",SPHEROID[\"Unknown\",6370000,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]],PROJECTION[\"Lambert_Conformal_Conic_2SP\"],PARAMETER[\"standard_parallel_1\",45],PARAMETER[\"standard_parallel_2\",22],PARAMETER[\"latitude_of_origin\",38.7],PARAMETER[\"central_meridian\",20],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]]";
  :GeoTransform = "-1088999 6000 0 927000 0 -6000 ";

It seems that it gets a semi_major_axis and not an earth_radius (for the spherical datums).

Do I have to define anything else in gdal_translate command?


texts: How to connect Dotspatial with PostGIS/PosgresSQL database and perform SQL function

texts:
I am currently doing a project that needs a GIS front-end interface to visualize geographical objects (road networks), and a back-end geo-spatial database to implement network analysis (SQL based). Our development platform is based on C#. I understand QGIS is powerful and available but as we really want to build an independent platform so we do not consider it.

Now, the progresses include:

I have built a front-end Dotspatial based GIS interface, it can display the road networks and query them as we want.
I have installed the PostgresSQL and its spatial extension PostGIS.

However, the current issue is how to connect the Dotspatial GIS platform with PostGIS/PostgresSQL in order to perform the analysis we need. What I really want is to do all activies at interface side. For example, by selecting a button the software do analysis automatically, and everything is going without other operation at PostGIS database side.

How can I start? Any information or helps? Many thanks!


texts: Is there any database or API that can provide data present in Atlas books for example, list of countries, population, temperature etc.


texts: Calculating watershed runoff with landuse and rainfall data?

texts: What is the method for computing rasters of total rainfall of a given month, flow accumulation, and land use (runoff coefficient) or soil K factor? I am using ArcMap 10.


texts:
I have these LIDAR data *.xyz files (1st col x, 2nd col y, 3rd col z)

485000.16 6508997.48 1.84
485001.27 6508998.92 1.86
485001.02 6508996.68 1.85
485000.77 6508994.45 1.79
485000.52 6508992.19 1.79
485000.02 6508987.69 1.74
485000.23 6508974.93 1.89
485001.02 6508981.87 1.86
485001.54 6508986.50 1.72
485002.32 6508993.34 1.72

...and I want to import them to GRASS GIS. I've tried with one file but r.in.xyz doesn't seem to work. All the program says is that:

r.in.xyz
input=C:\Users\Siim\SkyDrive\Dokumendid\kool\magistritoo\lidarandmed\filtreeritud\508485.xyz output=508485 fs=
Reading data ...
Writing to map ...
r.in.xyz complete. 0 points found in region.
(Mon Feb 09 20:17:26 2015) Command finished (1 sec)


texts: Openlayers 3 with cordova android has 10m cap

texts:
I have been using the standard geolocation plugin with my apps for a few years now.  a few months ago I noticed that all of my devices started to get really inaccurate.  And then this week a put the two together and realized it was the plugin which was causing the issue.  I'm not sure if I accidentally updated it or what, but it now uses the browser itself instead of a native Android approach like before.  Which means that it caps the accuracy to 10m (i'm not sure why).

I tried to revert to an older version and also tried another plugin which uses the Android approach which did solve the accuracy cap, but created a new problem in that there would be rendering issues where (mostly on newer devices) the map would just go blank when moving it to the user's location.

Can anyone shed any light on my situation?  If you have an OpenLayers map rendering well with cordova without the 10m cap, could you please share what you are using for it?  By that I mean, plugin name and version, plus what ever else may be relevant (cordova cli version?).


texts: Exporting MODIS imagery in Google Earth Engine & Pixel colours

texts:
I am attempting to detect fires using MODIS thermal anomaly data (MODIS/006/MYD14A1) on Corsica island, and the code I am using below produces a suitable fire mask. I am a bit inexperienced with this, and am unsure what the different pixel colours correspond to (any information relating to pixel colours online is evading me).

Also, when I attempt to export the image to drive I am getting an error message in the Tasks tab: Error: Image.clipToBoundsAndScale, argument 'input': Invalid type. Expected type: Image. Actual type: ImageCollection.

var site = ee.Geometry.Point(9.41, 42.81);
var aoi = ee.Geometry.Rectangle(9.317758730652592, 42.86545558133746, 9.504869631531498, 42.77429487610457);
Map.centerObject(aoi, 11);

// Mid Fire period
var midStart =ee.Date('2017-08-08');
var midEnd =ee.Date('2017-08-24');

//MODIS images of the site,
var midfireimage = ee.ImageCollection('MODIS/006/MYD14A1')
                    .filterBounds(site)
                    .filterDate(midStart,midEnd);

Map.addLayer(midfireimage, {bands: ['MaxFRP', 'FireMask', 'FireMask'], min:0.0, max:6000}, 'Fire Mask');

Export.image.toDrive({
  image:midfireimage.select('MaxFRP', 'FireMask', 'FireMask'),
  region:aoi,
  description: 'MidFire',
  folder: 'Corsica_Fire',
  scale:30});


texts: How do you manage columns in an ArcSDE/Oracle feature class?

texts:
Our enterprise has a long history with Oracle and ArcSDE. And lots of tension between the Oracle people and GIS people. Oracle folk often don't like the requirement of using expensive 'nonstandard' tools (e.g. anything other than TOAD or SQLPLUS) to manipulate business tables. But I digress already!

So for example today I wanted to do something simple - a column in a 15-year old non-versioned low-precision sdebinary feature class was marked as non-nullable. Recent editing clients don't use that column any more, but for various reasons we don't want to drop it. So I wanted to change the column to nullable.

1) do it in TOAD easy peasy. Even if the feature class is locked by ArcMap sessions or ArcGIS Server map services. That worked. BUT - this does not update the geodatabase system tables. SDETABLE -o describe shows consistent with with what TOAD shows, but ArcGIS still sees the column as non-nullable. The column_registry.object_flags column (see Appendix E, Administering ArcSDE for Oracle) shows '4' for this column. Whatever that means.

2) do it in the SDE commandline tools - SDETABLE -o alter_column.  That worked too, but first required me to identify locks on the feature class, track down which ArcGIS Server services include layers connecting as those users, stop them, then run the command. Lots of literal running around. AND - this also does not update the geodatabase system tables, so ArcGIS still sees the column as non-nullable.

3) use ArcCatalog 10.1 feature class properties dialog? In my experience everything it tells you about an ArcSDE column is read-only. So I guess not.

4) a Toolbox tool? I didn't find a handy 'change column metadata' tool in 10.1.

4) anything else I missed? What will you do when the commandline tools go away? Are things any better in 10.2?


texts: Making Merge Tool run just once in ModelBuilder?

texts:
I have dozens of point shapefiles I want to merge. First I used the Feature Classes Iterator and declared the iteration output as input for the merge tool directly, the problem was that it merged the output files with themselves and the result was only the last shapefile from the folder.

What I did next, was to use the Model Builder's Collect Values tool, it did the job just fine, but took so much time because the model collects and the merges all the input feature classes once for each input file. (For example, I have 30 point shapefiles, the model collects all those 30, 30 times, and then merges 30 times).

My question is, how can I make the collect and/or the merge tool run just once? The model works fine but it'd be way better if it was more efficient. The solution may be conspicuous and simple but I can't seem to figure it out.


texts:
I have a line with several points along it. I am hoping to name these points numerically, from 1 to the last point along the line.
I've so far used Network Analyst to find my distance along the line for all points and given myself a way to organize it by location. What I am attempting to do here is name them "1, 2, 3, 4..." etc. from start to finish, but am not sure how.


texts: Reprojecting coordinates and geoJSON in leaflet

texts:
I have a Leaflet map in default projection and a WMS layer also in default (Web Mercator) projection. I would like to add polylines and a geoJSON layer to my map, but the coordinates are in EPSG4326. I have figured out how to set the entire map projection, passing L.CRS.EPSG4326 into the map object. But as I have a WMS in Web Mercator (no other projection options here), but polylines and geoJSON in EPSG4326, that does not seem to be an option.

I looked at the Proj4Leaflet plugin, but that also seems to pass a CRS into the map object.
Any ideas?


texts: I am working on labeling features. I tried the annotation labels, but I have decided that it is not what I wanted. How do I change it back or delete the annotation label?


texts: Dissolving polygons based on attributes with Python (shapely, fiona)?

texts:
I have been trying to create a function that does basically the same thing that the QGIS "dissolve" function. I thought it would be super easy but well apparently not. So from what I have gathered around, the use of fiona with shapely should be the best option here. I just began to mess around with vector files so this world is pretty new to me and python as well.

For these example, I am working with a county shapefile founded here http://tinyurl.com/odfbanu
So here are some piece of code I gathered but can't find a way to make them work together

For now my best method is the following based on : https://sgillies.net/2009/01/27/a-more-perfect-union-continued.html. It works fine and I get a list of the 52 states as Shapely geometry. Please feel free to comment if there is a more straight forward way to do this part.

from osgeo import ogr
from shapely.wkb import loads
from numpy import asarray
from shapely.ops import cascaded_union

ds = ogr.Open('counties.shp')
layer = ds.GetLayer(0)

#create a list of unique states identifier to be able
#to loop through them later
STATEFP_list = []
for i in range(0 , layer.GetFeatureCount()) :
    feature = layer.GetFeature(i)
    statefp = feature.GetField('STATEFP')
    STATEFP_list.append(statefp)

STATEFP_list = set(STATEFP_list)

#Create a list of merged polygons = states
#to be written to file
polygons = []

#do the actual dissolving based on STATEFP
#and append polygons
for i in STATEFP_list :
    county_to_merge = []
    layer.SetAttributeFilter("STATEFP = '%s'" %i )
    #I am not too sure why "while 1" but it works
    while 1:
        f = layer.GetNextFeature()
        if f is None: break
        g = f.geometry()
        county_to_merge.append(loads(g.ExportToWkb()))

    u = cascaded_union(county_to_merge)
    polygons.append(u)

#And now I am totally stuck, I have no idea how to write
#this list of shapely geometry into a shapefile using the
#same properties that my source.

So the writing is really not straight forward from what I have seen, I really just want the same shapefile only with the country dissolve into states, I don't even need much of the attribute table but I am curious to see how you can pass it on from the source to the new created shapefile.

I found many pieces of code for writing with fiona but I am never able to make it work with my data.
Example from How to write Shapely geometries to shapefiles? :

from shapely.geometry import mapping, Polygon
import fiona

# Here's an example Shapely geometry
poly = Polygon([(0, 0), (0, 1), (1, 1), (0, 0)])

# Define a polygon feature geometry with one attribute
schema = {
    'geometry': 'Polygon',
    'properties': {'id': 'int'},
}

# Write a new Shapefile
with fiona.open('my_shp2.shp', 'w', 'ESRI Shapefile', schema) as c:
    ## If there are multiple geometries, put the "for" loop here
    c.write({
        'geometry': mapping(poly),
        'properties': {'id': 123},
    })

The problem here is how to do the same with a list of geometry and how to recreate the same properties than the source.


texts: QGIS - Can't import dwg/dxf or vector layer ("Could not open layer list")

texts:
I've tried importing a dwg file into QGIS. Either I've tried importing it as DWG or DXF in all the different versions possible. I've also tried adding the file as a vector layer but it only imports "entities" without anything showing up. Is that because my polylines are not georeferenced despite being in the correct location?
Importing any version of the file gives me the "Could not open layer list" error.
I've also tried different lisp's to export the xy coordinates of the vertexes of my polylines as well, without any success.


texts: How to create a standard deviation map from three rasters in QGIS?

texts: I have 3 normalised (0-100) raster maps. Is it possible to get a map which shows the standard deviation? I'd like to see the areas with smaller and higher standard deviation.


texts: Loading NetCDF file into QGIS, file is flipped upside down

texts:
I have a georeferenced NetCDF file that I am trying to load into QGIS (3.14). I can open the file and load it but the file always renders upside down. I am not sure if this is a problem with my NetCDF file (incorrect georeference) or if this is a bug with QGIS and the GDAL NetCDF driver. I have been able to successfully plot the NetCDF file with a few other programs (ncview, python+matplotlib+basemap) and it renders correctly, so I am thinking this is an issue with QGIS. Is there a way to tell QGIS how to treat the file?
I have tried processing my file with gdal_translate and experimenting with -co WRITE_BOTTOMUP=YES and -co WRITE_BOTTOMUP=NO, as follows:
gdal_translate -of netcdf -co WRITE_BOTTOMUP=NO NETCDF:"vegetation.nc":veg_class /tmp/OUTPUT_bu_no.nc

But in both cases the resulting file still loads upside down (and interestingly in both cases ncview then renders the file upside-down as well).
I also discovered another option for the GDAL netcdf driver (here), GDAL_NETCDF_BOTTOMUP=[YES/NO], but I am unsure how to send this option to QGIS when it opens the file?
Here is the metadata for the netcdf file I am trying to load, perhaps someone can tell if the georeference info is incorrect:
ncdump -h /<path to my file>/vegetation.nc
netcdf vegetation {
dimensions:
    Y = 50 ;
    X = 135 ;
variables:
    int veg_class(Y, X) ;
        veg_class:grid_mapping = "albers_conical_equal_area" ;
    float lat(Y, X) ;
    float lon(Y, X) ;
    int y(Y) ;
        y:standard_name = "projection_y_coordinate" ;
        y:long_name = "y coordinate of projection" ;
        y:units = "m" ;
    int x(X) ;
        x:standard_name = "projection_x_coordinate" ;
        x:long_name = "x coordinate of projection" ;
        x:units = "m" ;
    char albers_conical_equal_area ;
        albers_conical_equal_area:grid_mapping_name = "albers_conical_equal_area" ;
        albers_conical_equal_area:false_easting = 0. ;
        albers_conical_equal_area:false_northing = 0. ;
        albers_conical_equal_area:latitude_of_projection_origin = 50. ;
        albers_conical_equal_area:longitude_of_central_meridian = -154. ;
        albers_conical_equal_area:standard_parallel = 55., 65. ;
        albers_conical_equal_area:longitude_of_prime_meridian = 0. ;
        albers_conical_equal_area:semi_major_axis = 6378137. ;
        albers_conical_equal_area:inverse_flattening = 298.257222101 ;
        albers_conical_equal_area:spatial_ref = "PROJCS[\"NAD83 / Alaska Albers\",GEOGCS[\"NAD83\",DATUM[\"North_American_Datum_1983\",SPHEROID[\"GRS 1980\",6378137,298.2572221010002,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"6269\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4269\"]],PROJECTION[\"Albers_Conic_Equal_Area\"],PARAMETER[\"standard_parallel_1\",55],PARAMETER[\"standard_parallel_2\",65],PARAMETER[\"latitude_of_center\",50],PARAMETER[\"longitude_of_center\",-154],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AUTHORITY[\"EPSG\",\"3338\"]]" ;
        albers_conical_equal_area:GeoTransform = "-137256.2642264267 1000.050170338725 0 2372459.286491414 0 -999.7969847651308 " ;

// global attributes:
        :Conventions = "CF-1.5" ;
        :source = "/home/UA/tcarman2/dvm-dos-tem/./scripts/create_region_input.py::fill_veg_file --xoff -137256.290295 --yoff 2321969.51702" ;
}


texts: When users add my WMS server to their ArcGIS client, the first 10 layers are enabled by default and the others are disabled by default. Is it possible to set Geoserver so the WMS has all layers disabled by default?


texts: Rasterization - choosing attribute from attribute table for cell value of raster - why?

texts:
When a polygon shapefile is converted into a raster, an attribute from
  the shapefile's attribute table must be chosen for the cell value of
  the raster.

I am pretty sure the above statement is true, but I cannot explain why.

Can anyone explain why this is the case?


texts:
It is said, that The 6° wide UTM zones, numbered 1–60.

But where is the origin? Zone 1 is from longitude [-3..+3], [0..6], [-180..-174]? Where is it said?


texts: Performing Sensitivity Analysis on Suitability Model using ArcGIS ModelBuilder?

texts:
I am working on a suitability analysis model for a landfill site selection. Basically, I am using GIS-based MCDA so I have 11 factor layers/maps that I considered for the analysis. I have already standardized these factor maps. Now, I am using the Weighted Sum tool to aggregate/combine these layers, each layers having a corresponding weight. What I want to do in sensitivity analysis is to iterate the weighted sum in my model where each iteration uses a different combination of weights for each factor layer.

How do I do this in ModelBuilder?

The combination of weights that I want to use for each iteration of the weighted sum is shown in the attachment below.Of course, I still need to do way more iterations than what I showed here.


texts: Is it necessary to add Tile layer background?

texts:
Is it necessary to add tile layer for making a GeoJSON data rendered in web browser using Leaflet or any other API ?

I have a geoJSON data (I obtained the GeoJSON data by converting .shp file to GeoJSON using ogr2ogr) of a small hypothetical network which I want to render in the browser. My network when opened in QGIS is shown

For the above network it makes no sense to add a Tile layer of OSM or any other service. Is there a way I can render the network shown above in the web using leaflet or any other API ?

I'm very new to GIS related thing.


texts: Split a GeoTIFF (BigTIFF) file into subsets with a specific overlap in Python

texts:
I have a BigTIFF file that I need to split into tiles with a set tile size and overlap. I have a script for this using PIL:

tile_height = tile_width = 1000
overlap = 80
stride = tile_height - overlap
start_num=0

def crop(infile, tile_height, tile_width, stride, img_dict, prj_name):
    im = Image.open(infile)
    img_width, img_height = im.size
    print(im.size)
    print(img_width * img_height / (tile_height - stride) / (tile_width - stride))
    count = 0
    for r in range(0, img_height-tile_height+1, stride):
        for c in range(0, img_width-tile_width+1, stride):
            #tile = im[r:r+100, c:c+100]
            box = (c, r, c+tile_width, r+tile_height)
            top_pixel = [c,r]
            img_dict[prj_name + "---" + str(count) + ".png"] = top_pixel
            count += 1
            yield im.crop(box)
img = Image
img_dict = {}

# create the dir if it doesn't already exist
if not os.path.exists(img_dir):
    os.makedirs(img_dir)

# break it up into crops
for k, piece in enumerate(crop(infile, tile_height, tile_width, stride, img_dict, prj_name), start_num):
    img=Image.new('RGB', (tile_height, tile_width), (255, 255, 255))
    print(img.size)
    print(piece.size)
    img.paste(piece)
    image_name = prj_name + "---%s.png" % k
    path=os.path.join(img_dir, image_name)
    img.save(path)

#add a json file with all image names and geospatial metadata
full_dict = {"image_name" : infile,
            "image_locations" : img_dict,
             "crs" : str(dataset.crs)
            }

with open(img_dir + '/data.json', 'w') as fp:
    json.dump(full_dict, fp)

I can't use PIL on my other rasters, as they are "BigTiff" files and not supported in PIL. I am looking for a way to translate this script into another module keeping these exact parameters. I need the parameters and naming methods to stay exactly the same as I'm using these tiles for a deep learning model that I have already created.

I have never used something like GDAL before, but I've read that this may be my best bet for Big TIFF tiling? I would really like to find a way to do this in Python.


texts:
How can external Geojson feeds be styled in leaflet? To add my external feed I'm  using the snip below:

var mygeojson = new L.geoJson();
mygeojson.addTo(map);

$.ajax({
    dataType: "json",
    url: "https://linktomyexternalgeojsonurl",
    success: function(data) {
        $(data.features).each(function(key, data) {
            mygeojson.addData(data);
        });
    }
}).error(function() {});

The Leaflet website offers examples here http://leafletjs.com/examples/geojson.html but when following the map breaks and does not display any marker. I am clearly doing something wrong and it is likely very simple and obvious to anyone but me.

Leaflet example:

var myStyle = {
    "color": "#ff7800",
    "weight": 5,
    "opacity": 0.65
};

L.geoJson(myLines, {
    style: myStyle
}).addTo(map);

My implementation of the Leaflet example:

var myStyle = {
    "color": "#ff7800",
     "weight": 5,
     "opacity": 0.65
    };

var mygeojson = new L.geoJson({
    style: myStyle})
    mygeojson.addTo(map);

    $.ajax({
        dataType: "json",
        url: "https://linktomyexternalgeojsonurl",
        success: function(data) {
            $(data.features).each(function(key, data) {
                mygeojson.addData(data);
            });
        }
    }).error(function() {});


texts:
I have downloaded a set of premade styles (*.qml). Where do I have to store this files for usage?. I am using QGIS 1.7 for windows.
thank you!
Updates
@Underdark,
thanks for answering, if found the following "solution": In the window "Layereigenschaften " (Layer-properties) click the button "Stil laden ..." (load style),now you can use this imported style. But I found no way to add these styles permanently to qgis-stylemanager
Btw: you can download this "open-street map"-styles from this site: http://www.sourcepole.ch/2010/2/28/spatialite QGIS-Styles für OpenStreetMap-Daten: Download
@Nathan, thanks for answering, yes I did mean qml


texts: In QGIS how do I edit/save data which has no geometries?

texts:
In QGIS I can import data tables (CSV in particular) which have no inherent geometry (and no lat/lon columns) into QGIS. There doesn't seem to be any obvious way to manipulate that data (edit/save).

(Obviously this data has a reference column which links to data which does have geometry).

It's documented that (at QGIS 2.6) it's not possible to edit/save CSV files:

Is there no editmode for CSV files?
How do I delete multiple data points from a .csv file through QGIS?
Why can't I edit attribute table imported through 'text as layer'

I am aware that I could edit data in another application, repeatedly returning to this, saving as CSV and importing into QGIS, but this seems very clumsy and inefficient. I'd like to remain within QGIS.

I'm guessing that there are many workarounds for this which are simple for those experienced with using databases,  but on other GIS software it's a trivial task to work with a geometry-free table.

How do I set-up a spatialite database for this purpose?


texts: Extract time series for features from large masked xarrays

texts:
I have done this before with GDAL and Python, but was hoping to be able to use xarray to scale things up.

Let's say I have three datasets:
1. A large multiband geotiff where each band is a different date of the same parameter. Spatial resolution is r0 (say). ds_A say
2. A large multiband geotiff where each band is a value that is compared to a threshold and used to mask ds_A. Times and spatial extents are the same as ds_A, but resolution here is r1. Let's call this ds_B.
3. A vector layer with polygons. Let's call it vect. Lots of features in vect (>1e3, 1e4)

I want to (i) mask ds_A with ds_B and then extract (say) the spatial average of ds_A for all times for all the features in vect.

My approach so far is this, but this takes forever, in terms of memory and speed.

def do_mask(raster, vector, attr=None):
    """Either creates a 0/1 mask (if `attr` is `None`), or a
    rasterises the label `vector` according to a field in `attr`,
    0 otherwise."""
    g = gdal.Open(raster)
    w, h = g.RasterXSize, g.RasterYSize

    out_ds = gdal.GetDriverByName("MEM").Create("", w, h, 1, gdal.GDT_Byte)
    out_ds.SetProjection(g.GetProjection())
    out_ds.SetGeoTransform(g.GetGeoTransform())
    g = None
    if attr is None:
        gdal.Rasterize(out_ds, vector, burnValues=1, allTouched=True)
    else:
        gdal.Rasterize(out_ds, vector, attribute=attr, allTouched=True)
    return out_ds.ReadAsArray()

def mask_crop_stack(raster, vector, chunks, attr=None):
    """Create either a mask where a vector is True, or a set of labels
    from all features in a vector."""
    # Open the original dataset to get the coordinates and stuff
    # for rasterisation
    ds0 =  xr.open_rasterio(raster, chunks=chunks)
    # Create the mask or labels field.
    mask = do_mask(raster, vector, attr)
    # Create a DataArray with the mask/labels, and copy the metadata from the
    # original or target array.
    mask_ds = xr.DataArray(mask, coords={'x':ds.x, 'y':ds.y}, dims=list(ds.dims)[1:])
    # Throw away unmasked stuff to make the data more manageable
    cropped_ds = ds.where(mask_ds > 0, drop=True)
    return cropped_ds

# Open first file, and crop to vector area
ds_A = mask_crop_stack("ds_A.tif",
                       "vect.geojson",
                       chunks={"band":1, "x":1024, "y":1024})
# Create labels DataArray
labels = mask_crop_stack("ds_A.tif",
                         "vect.geojson",
                         chunks={"band":1, "x":1024, "y":1024},
                         attr="ID_PARCEL")
# For some reason, the extents are slightly different
# So interpolate to match ds_A
labels = labels.interp(x=ds_A.x, y=ds_A.y, method="nearest")
# Open mask dataset, and crop spatially
ds_B = xr.open_rasterio("ds_B.tif",
                       "vect.geojson", chunks={"band":1, })
# Interpolate to match ds_A (change resolution from r1 to r0, say)
ds_B_rsmpl = ds_B.interp(x=ds_A.x, y=ds_A.y, method="nearest")
thresh = thresh_me
# Apply the mask in ds_B to ds_A
ds_A = ds_A.where(ds_B_rsmpl < thresh, other=np.nan)
# Apply a sanity check to ds_A, values need to be between 0 and 10000
ds_A = ds_A.where(ds_A > 0, other=np.nan)
ds_A = ds_A.where(ds_A < 10000, other=np.nan)

# Extract the time series and store in a dictionary
tseries = {label:ds_A.where(labels==label, other=np.nan).mean(axis=(1, 2)).to_pandas()
          for label in np.unique(labels) if label > 0}
# **Ideally**, I'd do processing pixel by pixel,
# `apply_ufunc` might not work as processing outputs large number
# of parameters per pixel.

The last line of the previous code takes forever, and often crashes the Python interpreter. I'm very new to xarray, so I'm guessing that I'm doing something wrong. My GeoTIFFs are tiled and blocky. I'm sure I'm missing something obvious here!


texts: QGIS clone QgsLayoutItemLegend or copy existing legend object

texts:
I have some legend object with my custom setting and styles:
exist_legend=QgsLayoutItemLegend()

How I can copy or clone or copy all setting and styles to new object legend, like this:
new_legend=exist_legend.clone()


texts: Are there any simple models for wildfire hazard mapping with ArcGIS ModelBuilder?

texts: I am interested in producing a simple fire hazard/severity map via ModelBuilder from commonly available datasets (e.g. DEM, NDVI, slope, aspect, land use/cover, etc).  Are there any resources or examples of simple models used to produce maps such as the example below?  Specifically, I am looking for a model that I can construct within about 10-15 minutes during a presentation to produce an interesting graphic.


texts: How to center the QGIS map on the maximum value of a raster?

texts:
I have a very large QGis raster layer showing the differences between two nearly identical rasters.

Because of this, it is very difficult to find where in this layer the few pixels that are different occur.

I have thought that one way to find them might be to have QGis zoom in and center on the maximum point of the layer. Is there a way to do this?

Other suggestions are welcome.


texts: Differential correction for certain days in the past

texts: I would like to know if there is any way to perform differential corrections on data that were collected about 1-3 months before. I know that NOAA has data here, but can I perform differential corrections in pathfinder (or by any other mean) with old data?


texts:
First things first! I am using Ubuntu Mate on 15.04. I installed Qgis and SAGA via the http://ppa.launchpad.net/ubuntugis/ubuntugis-unstable/ubuntu packages. So I use QGIS 2.12 and SAGA 2.1.4.

I wanted to use the cluster analysis for grid tool, available in the SAGA-toolbox and I couldn't find it. It should be under shapes->table but I checked and the entire table-tools are missing. Others might be missing too (Qgis is showing 235 SAGA-tools). The ones which are there are fully functionable.

In the official SAGA-programm the tool is visible though. So how can I port the missing tools into QGIS?


texts: ERROR: "Can't save to disk (check for space, write permission, or name)" while creating TIN in memory

texts:
I get the following error while creating TIN in memory
: ERROR 999999: Error executing function.
Can't save to disk (check for space, write permission, or name).
Failed to execute (CreateTin).

I don't have any issue with writing any feature class in memory, I wonder why it doesn't work for TIN.

TIN = arcpy.CreateTin_3d("in_memory/TIN", Coordinate_System, "XS_lyr P001 hardline FID", "DELAUNAY")

When I remove the "in_memory/" and just save it on local drive it works fine.


texts: Removing Linear Referencing labels after removing layers?

texts:
How does one remove Labels after adding one from Linear Referencing tab in the Find Window of ArcMap 10.4.1? Removed feature but the labels remain

The labels in the find window and under the Linear Referencing tab are displayed even though the Feature has been removed


texts: Moving a memory layer in the TOC using PyQGIS

texts:
I've got a layer (difflayer) that has been created in memory. Basically, all the layers I've got in this project are in memory and the layer creation/manipulation is all being handled via a plugin.

Now, from looking around it seems the best way to move a layer around is to clone() it, use insertChildNode to add the clone, and then delete the original. Which I'd have no problem doing, except that layers in memory don't support clone().

Does anyone have any ideas? I've also tried to duplicate the layer using the below code but get the following error:

feats = [feat for feat in difflayer.getFeatures()]
mem_layer = QgsVectorLayer("Polygon?crs=EPSG:4283&field=ID:integer", "dup layer", "memory")
mem_layer_data = mem_layer.dataProvider()
attr = difflayer.dataProvider().fields().toList()
mem_layer_data.addAttributes(attr)
mem_layer.updateFields()
mem_layer_data.addFeatures(feats)
root = QgsProject.instance().layerTreeRoot()
root.insertChildNode(2,mem_layer)

TypeError: QgsLayerTreeGroup.insertChildNode(int, QgsLayerTreeNode): argument 2 has unexpected type 'QgsVectorLayer'


texts: How to improve speed of Export to PDF by ArcGIS for Desktop?

texts:
I am working with the Fire Station Wall Map template in ArcGIS 10.1. My export to PDF is taking upwards of 30-45 minutes. Using the default PDF settings in the export menu.

My map is 60" x 36" so its not too big. My resulting file size is approx 7mb.

Any idea how I can improve the speed?


texts:
Is there a way to access the classification histogram which is used when classifying data in ArcGIS desktop?

The background is I have a series of data layers for which I need to use the same classification. What I'd like to do is to be able to show the histogram for each layer in one corner of the map - possibly even with the bars coloured in the same way as the symbols.

Has anyone encountered this or does anyone have any suggestions as to how it might be done?


texts: Google Maps V3 + multiple Fusion Table layers = inconsistent markers?

texts:
One of my colleagues is working on a map for a news story and is having some difficulty getting the markers to consistently show up. Changing the order in which the layers instantiate is causing significant variance in which markers are actually visible:

Compare: http://jsfiddle.net/aendrew/gTgHm/3/

...With: http://jsfiddle.net/aendrew/hLVzM/1/

For an obvious example, look to the right of the label for Glasgow — there are two green markers in the first version, three green in the second.

Is there some query limit I'm hitting with the Fusion Tables?

Note that each layer has less than 100 markers.


texts: I have not found shp2pgsql in postgis installation

texts:
I have been installed postgis following this guide. All is OK.

SELECT postgis_full_version();
"POSTGIS="2.1.7 r13414" GEOS="3.4.2-CAPI-1.8.2 r3921" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.10.1, released 2013/08/26" LIBXML="2.9.1" LIBJSON="UNKNOWN" TOPOLOGY RASTER"

But I have not the shp2pgsql  in /usr/lib/postgresql/9.4/bin as you can see:

foo@bar:/usr/lib/postgresql/9.4/bin$ ls
clusterdb   initdb             pg_dump         pg_standby      psql
createdb    oid2name           pg_dumpall      pg_test_fsync   reindexdb
createlang  pg_archivecleanup  pg_isready      pg_test_timing  vacuumdb
createuser  pg_basebackup      pg_receivexlog  pg_upgrade      vacuumlo
dropdb      pgbench            pg_recvlogical  pg_xlogdump
droplang    pg_controldata     pg_resetxlog    postgres
dropuser    pg_ctl             pg_restore      postmaster

Why this happens? and how can I install it?


texts: How to turn on all fields before appending?

texts: I'm using append on model builder, no test schema type. In order that field will match they need to be at the same order, and visble. How can I turn on all fields of input layer before appending inside model builder ?


texts: Processing(.org) for geospatial visualization?

texts:
Does anyone know if Processing is being used for geospatial visualization? I can't find much out there, but this is impressive: http://benfry.com/zipdecode/

Here's one I just found: http://fathom.info/projects/countyhealth.html


texts: XY Coordinate conversion for google maps

texts:
I have coordinates in XY format for System 22/17.

I wish to convert them to plot on google maps. Any help?


texts: Aggregating with attribute($currentfeature) filter in QGIS

texts:
Why doesn't the filter work (4th variable passed into the aggregate function) when I try to use in the field calculator an expression of this type:

aggregate( 'Layer', 'sum', "Field1" ,  "Field2"  = attribute( $currentfeature , 'Field2'))

It returns the total sum of Field1 rather than applying the filter.
In the aggregate function, such a filter seems to be always true, while by itself, it works.

Is it not possible to use feature-based expressions with the aggregate function, or can it work if I modify the expression?


texts: ArcPad 10: Accessing the Geography information of a feature through ArcPad Studio (VB Script)

texts:
I am trying to access the geography feature property of a feature being edited. I need to do it programatically (in order to auto populate a field using the geography) using VB Script. I cannot seem to figure out what the ArcPad object is that I want to access in the script. Any suggestions?

Thanks

EDIT

How do I access the Points Collection of an added feature in ArcPad studio?


texts: Adding coordinate constant to file in QGIS

texts:
i have a dxf file where the coordinate constants has been dropped from x coordinate. For example my real coordinate should read +2954632 but instead it reads as 54632.

How can I add the 29 000 00 in QGIS?


texts: QGIS Python plugin crashes when using processing algorithms and QThread

texts:
I'm developing a QGIS plugin which calls several processing algorithms. I call them in a pipeline where the input layer of each algo is the output of the previous one. The pipeline is called on each feature of a selected layer where the first step creates a temporary layer in the memory from a feature.

As it's a UI plugin and the algorithms can take some time, I don't want to lock the UI and try to use threads. However, using QThreads randomly crashes QGIS.

Here is a minimal code of my plugin.

def CreateFeatureLayer(input_layer, feature):
  temp_layer_name = "feature_layer_" + str(feature.id())

  feature_layer = QgsVectorLayer("Polygon?crs=EPSG:4326", temp_layer_name, "memory")
  feature_layer_data = feature_layer.dataProvider()

  attr = input_layer.dataProvider().fields().toList()
  feature_layer_data.addAttributes(attr)
  feature_layer.updateFields()

  feature_layer_data.addFeature(feature)
  return feature_layer

def BufferFeature(layer):
  result = processing.run("native:buffer",
  { 'INPUT':layer,
    'DISTANCE':0.0008,
    'SEGMENTS':5,
    'END_CAP_STYLE':0,
    'JOIN_STYLE':0,
    'MITER_LIMIT':2,
    'DISSOLVE':False,
    'OUTPUT':'TEMPORARY_OUTPUT'})

  buffer_layer = result['OUTPUT']
  return buffer_layer

def MakeDiff(input_layer, buffer_layer):
  result = processing.run("native:difference",
  { 'INPUT' : buffer_layer,
    'OUTPUT' : 'TEMPORARY_OUTPUT',
    'OVERLAY' : input_layer })

  diff_layer = result['OUTPUT']
  return diff_layer

class Worker(QObject):

  start = pyqtSignal()
  success = pyqtSignal(str)

  def __init__(self, selected_layer):
    QObject.__init__(self)
    self._selected_layer = selected_layer

  @pyqtSlot()
  def run(self):
    for f in self._selected_layer.getFeatures():
      feature_layer = CreateFeatureLayer(self._selected_layer, f)
      buffer_layer = BufferFeature(feature_layer)
      diff_layer = MakeDiff(feature_layer, buffer_layer)
    self.success.emit("Success")

class Plugin:

  def WorkerFinished(self, string):
    self._my_worker.deleteLater()
    self._my_thread.quit()
    self._my_thread.wait()
    self._my_thread.deleteLater()

  def run(self):
    selectedLayerIndex = self.dlg.comboBox.currentIndex()
    selectedLayer = self._selectedLayers[selectedLayerIndex].layer()

    self._my_thread = QThread()

    self._my_worker = Worker(selectedLayer)
    self._my_worker.moveToThread(self._my_thread)
    self._my_worker.success.connect(self.WorkerFinished)

    self._my_thread.started.connect(self._my_worker.run)
    self._my_thread.start()

Some runs successfully finishes, but there are times where it crashes QGIS randomly.

Calling the pipeline in the main thread (in Plugin.run()) doesn't crash QGIS.

Does anyone know what can cause the random crashes?


texts: Copying a feature from a shapefile to another while updating a field

texts:
I have two shapefiles.
After iterating on all features of one of them, the code reads a particular field value:

 from osgeo import ogr
 driver = ogr.GetDriverByName('ESRI Shapefile')
 SHP_F = myshaoefilefile1.shp
 SHP_CL = myshapefile2.shp

 SHP_F_DataSource = driver.Open(SHP_F, 1)
 SHP_CL_DataSource = driver.Open(SHP_CL, 1)
 # get the shapefiles layer
 SHP_F_layer = SHP_F_DataSource.GetLayer()
 SHP_CL_layer = SHP_CL_DataSource.GetLayer()
 # estraction of features from  layers
 SHP_F_featureCount = SHP_F_layer.GetFeatureCount()
 SHP_CL_featureCount = SHP_CL_layer.GetFeatureCount()
 for i in range (0, SHP_F_featureCount ):
    for j in range (0, SHP_CL_featureCount ):
          if SHP_CL_layer.GetFeature(j).GetField(4) in ['string1', 'string2', 'string3']
             # copy the feature from Layer2 to Layer1 ??
          if SHP_CL_layer.GetFeature(j).GetField(4) in ['string4', 'string5', 'string6']
            # update some field of the inLayer1.GetFeature(i) ??

If the field value is equal to a pre-defined string, I need to copy that feature to the first shapefile.
Also I need to update some fields, both integer and string type.
Is there a simple way to copy features from a shapefile to another? The two shapefiles have the same geometry type and the same projection and reference system (WGS84 UTM 33N).
The code is part of a batch python script (python 2.7 version).


texts: I have converted to ESPG 4326, but I am still getting an error that my coordinate system is wrong?

texts:
I created a shapefile in ArcGIS and I converted my CRS in QGIS to EPSG: 4326, WGS 84. I have done this method before and it allows me to plot my shapefile over a basemap in python. However, I am continuously getting an error that is meant when I have not converted my CRS to EPSG 4326, even though I have???

ValueError: shapefile must have lat/lon vertices  - it looks like this one has vertices
in map projection coordinates. You can convert the shapefile to geographic
coordinates using the shpproj utility from the shapelib tools
(http://shapelib.maptools.org/shapelib-tools.html)

Usually, when I change the CRS, the problem is solved but for some reason it is not fixing the issue.

I have followed this procedure many times, and it always works except for this time! http://fascicul.us/convert-display-shapefile-matplotlib.html

I have also converted my shapefile to EPSG: 4326, WGS 84 with ogr2ogr and I am still getting the same error when I attempt to plot it.

I can plot the shapefile with geopandas, but not with readshapefile (which is where I get the error). I would prefer to plot with readshapefile.


texts: how to rename the band-names of a layer stack?

texts:
I created a layer stack of 20 images (all GTiff-files) using the gdal_merge tool in QGIS, available at:

raster > miscellaneous > merge

The result is a GTiff layer stack of the 20 images, but the band names are band1, band2, band3...band20, looking at the properties.

Is there a way to rename the band-names, either afterwards, or during creation of the stack?

They are all from different sensors/dates and will be more meaningful to work with other names than band1,band2,band3.

Thanks!


texts: Converting points to raster and back again in QGIS?

texts:
To convert raster grids to vector points In QGIS there we have different possibilities (grass tool "to.vect"; saga tools "grid values to points").

My problem is a little bit complicated. In the beginning I have a point layer. For further data processing I need to convert the vector points to a raster file (every point is a raster now). In the end I want to convert this raster (raster points) backwards to a point layer for the kriging tool.

The grass tool (to.vect) is very slow and the saga tools (grid values to points) gives me more than just the points I need to have (I get a very high density point shape). Is there another solution just to convert the few rasters into a point layer?


texts: Converting images to different file formats using GeoTools

texts:
After searching through GeoTools documentation and following their tutorials, I have failed to find a class that can convert images to different format (for example translate a .tif to .jpg).

Does GeoTools have a Translate function, like GDAL, in some other git repository?


texts: How to create feature type from xml attribute value?

texts:
Here is a sample from an xml child element:

<element>
<aspects>

</aspects>
<aspects>

</aspects>
</element>

I want to extract the value "Grundlose Kuhle" from the aspects.a.id="1".v element.

How can I get FME to put all "v"-values from ID="1" elements into a list so that I can write them to a featuretype?


texts:
I assume I am making a very rookie mistake here, but I can't find an sDNA plugin in QGIS. I installed and authorised the software according to the instructions. I also checked "show also experimental plugins" box. Does sDNA needs to be installed in a specific folder in order for QGIS being able to find it?

On the other hand, ArcGIS is able to find the sDNA, but shows an error saying "click to repair project item" which dosn't seem to fix the problem.

I am running Windows 10, GIS version 3.4 Madeira, ArcGIS Pro and sDNA v4_0_3

Any help would be greatly appreciated.

Thanks!


texts:
I have two datasets representing the same area in Romania, but in different CRS, A is local, B is WGS 84.

dataset A with:

Stereo_70
Authority: Custom

Projection: Double_Stereographic
False_Easting: 500000,0
False_Northing: 500000,0
Longitude_Of_Origin: 25,0
Scale_Factor: 0,99975
Latitude_Of_Origin: 46,0
Linear Unit: Meter (1,0)

Geographic Coordinate System: GCS_Dealul_Piscului_1970
Angular Unit: Degree (0,0174532925199433)
Prime Meridian: Greenwich (0,0)
Datum: D_Dealul_Piscului_1970
  Spheroid: Krasovsky_1940
    Semimajor Axis: 6378245,0
    Semiminor Axis: 6356863,018773047
    Inverse Flattening: 298,3

and dataset B with:

WGS 84
WKID: 4326 Authority: EPSG

Angular Unit: degree (0,0174532925199433)
Prime Meridian: Greenwich (0,0)
Datum: WGS_1984
  Spheroid: WGS 84
    Semimajor Axis: 6378137,0
    Semiminor Axis: 6356752,314245179
    Inverse Flattening: 298,257223563

In order to convert B to A, I imported A's CRS to B=> Datasets did not overlap, they are miles away from each other.

Could anyone explain what exactly is the problem?


texts: Model built in ArcGIS 10 won't open in 10.1

texts: I've built a model to batch clip rasters.  It uses iteration and works fine on machines running arcGIS 10.  However, colleagues that have upgraded to 10.1 cannot open it.  the red x appears with a "unable to execute selected tool" pop up.  Below is a picture of the model (I can't export the script for some reason):


texts:
I have seen similar questions and read the responses but did not understand any of the answers even after reading them several times closely, including Adding XY Points to ArcMap?

I am adding XY data to an existing ArcMap document (ArcGIS Desktop 10.4). The existing feature class that I want to add the points on top of has a coordinate system of NAD_1983_StatePlane_Florida_East_FIPS_0901_Feet and a projection of Transverse_Mercator.

I added XY data, setting X Field = Longitude and Y Field = Latitude. I allowed the Coordinate System of the Input Coordinates to remain the default, which was Projected Coordinate System NAD_1983_StatePlane_Florida_East_FIPS_0901_Feet and Geographic Coordinate System  GCS_North_American_1983.

The points created by ArcMap are not in the correct location. I know this is a projection issue, but I am not sure how to fix it.

If you direct me to change or define projections, specify exactly which feature classes you mean - my existing feature class or the newly added XY points. I read one answer that said to change the projection of the data frame, but did not say how this is to be accomplished.

OLON         OLAT
-80.19112    26.11602
-80.15004    26.200484
-80.14676553 26.17532654
-80.19112    26.11602
-80.15004    26.200484
-80.20728    26.24123


texts: Random Forest out of the bag and confusion matrix

texts:
I have implemented Random Forest classifier to classify remote sensing data in R. The original code comes from here:

How to perform Random Forest land classification?

Everything works fine, but what I need is to obtain confusion matrix and out of the bag error for fast classification accuracy assessment. I have created junk dataset for validating functionality of the whole script itself (only 2 classes). I added parallel computing to speed up computations.

I am wondering if my implementation is correct and If I will be able to obtain confusion matrix and out of the bag error when I will add much more features describing classes later. Here is part of my code:

sdata <- readOGR(dsn = vyber_shp, layer = "trenink")

# extraction from raster to be classified

rdata <- data.frame(extract(r, sdata))

# create random forest model in parallel

cl <- makeCluster(detectCores())

registerDoParallel(cl)

rf.mdl <- foreach(ntree=rep(125, 4), .combine = combine, .packages = 'randomForest')
%dopar% {randomForest(x=rdata, y=sdata$Class, ntree=500,
proximity=TRUE, importance=TRUE, confusion=TRUE,
do.trace=TRUE, err.rate=TRUE)
    }

stopCluster(cl)

# classify raster in parallel

beginCluster()

predikce <- clusterR(r, raster::predict, args=list(model=rf.mdl))

endCluster()

# write and save result of classification as raster image

klasifikace <- writeRaster(predikce, filename='Klasifikace_RandomForest',
format='HFA', options='INTERLEAVE=BSQ', datatype='INT2S', overwrite=TRUE)

varImpPlot(rf.mdl)

rf.mdl$confusion


texts: Features with style categorized by date disappear

texts:
A shapefile of walking routes has attributes of name, date and length. If I choose Categorized style of display based on field "name" I can specify a different colour for each route and they appear as expected. If I try to categorize the features on the date field I can choose different colours for each date but the features disappear on the map.

What am I doing wrong?

Update: while creating Categories when you press the button for Classify it always creates one extra category for items which don't fit into any of the other categories. I normally remove this category by clicking on the minus sign icon. However, if I leave this category in all my lines appear on the map. If I remove the extra category all the lines disappear from the map.

Can anyone explain what is causing this?


texts: Using local tangent plane projection in QGIS?

texts:
Is local tangent plane projection is a kind of gnomonic projection?
is it possible to define it as a custom projection in QGis?


texts: Displaying SQL Server View (not table) on QGIS 3.18

texts:
I have an SQL Server database connected to QGIS 3.18. I am able to add tables including geography data type to the map without a problem. However, I want to be able to add a 'View' (Orn_View_spatial) with geography data type from SQL Server. It loads into the Browser panel successfully, but when I drag and drop it onto the map, QGIS returns an error message:
Layer is not valid: The layer dbname='MBA-LBA_Ornaments'
  host=DESKTOP-DLA7UV8 estimatedmetadata=true
  srid=4326 type=Point disableInvalidGeometryHandling='1'
  primaryKeyInGeometryColumns='0' table="dbo".
"Orn_View_spatial" (Geo4) is not a valid layer and can not be added to the map.
Reason: No primary key could be found on table Orn_View_spatial

Apparently it is not possible to create a primary key on a View in SQL Server--is there any way around this issue?


texts: Why is Format of output data grayed out on Save As dialog for Rectify?

texts:
I started georeferencing and I could not save my new "rectify" file, I don't know if is because of the raster or I did something wrong.

Then I find out that my control layer is not accurate, it doesn't match to other shapefiles in the same geodatabase (same projected coordinate reference), neither with a basemap.

How do I activate the Format of output data option on the Save As dialog for Rectify to make sure that I have an accurate control layer?


texts: Adding V-LATE extension to ArcMap 10.3 running on Windows 10?

texts:
I failed to add the V-late extension to ArcMap 10.3 running on Windows 10.
After trying to add the V_LATE_2beta_argis10.tlb in the customise mode with administration rights, I get the Message box "no new objects added".

Is there any known bug?


texts: Converting many kml files into one kmz and visualize them in GoogleEarth?

texts:
Is it possible one kmz file to have many kml files?
If yes, how can I convert the kml files to a kmz in order to be able to show them in GoogleEarth?
In case of one kml file, I read that its enough to place the kml in a folder, zip it and change the extension. I tried that and it works. But in the case of many kml files, it only shows the first one taken from the list:

When Google Earth opens a KMZ file, it scans the file, looking for the first .kml file in this list. It ignores all subsequent .kml files, if any, in the archive. If the archive contains multiple .kml files, you cannot be sure which one will be found first, so you need to include only one (link)


texts: Can you use variables in URI field in cartodb CSS?

texts:
I am wanting to use a different marker-file based on a value returned via SQL in cartodb.
I am wanting todo something like

marker-file:  url('[icon]');

Where the value in [icon] is a full URL to a PNG image.
When I do the above there is no error displayed, but there is no marker displayed on the map.
If I use the below then the marker is displayed.

marker-file: url('http://test.wildwalks.com/wildwalks_custom/icons/map/access/marker-green-seat-136.png');

The value returned by [icon] is identical to this full URL.
Any suggustions on how to use SQL values to define marker-file urls?


texts: Understanding background to Single Hash (#) Vs. Double (##) in QGIS Python Scripting?

texts:
First of all I understand the double hash (##) is used to define the Inputs / Outputs in a script's GUI.  If there's more to it than that please let me know.

In all other Python coding I've done double hashes are just another variation on the single hash (meant for comments).

Can anyone provide background information into how QGIS manages to not see a ## line as an unexecuted comment?

I assume it's simply built into the QGIS Python compiler to recognize ## as defining GUI inputs & outputs.

Again I understand it now, and bottom line is it works, but it struck me as an odd way to set the GUI.  Just curious.