texts: Creating loop in raster algebra which compares each cell to values of other rasters in time series?

texts:
I have a list of rasters representing NDVI values (and plan to do similar thing with some other index values as well) and in that raster I also set all values representing clouds to NoData.

Now, for each raster I want to calculate difference in values between that raster and previous one BUT if the previous raster has NoData value (so it's a cloud) in a cell, I want to skip that raster and get values from next previous raster (and if that raster also has NoData then skip to the next and so on) and then calculate delta values.

So with ArcPy it would look something like this:

def DeltaNDVI(pre, post):
    DNDVI = Raster(post) - Raster (pre)
    return DNDVI
for i in range (1, len(NDVIList)):
    preprepre = NDVIList[i-3]
    prepre = NDVIList[i-2]
    pre = NDVIList[i-1]
    post = NDVIList[i]
    DeltaRaster = Con(IsNull(pre), Con(IsNull(prepre), DeltaNDVI(preprepre, post),DeltaNDVI(prepre, post)), DeltaNDVI(pre, post))

This obviously isn't an optimal solution since the list I have contains a dozen rasters and there are multiple lists with different number of rasters in them. There should be some kind of a loop there but I don't think that I can use loop with ArcGIS map algebra.

The list is sorted based on a date and I have date of acquisition for each raster in date(day, month, year) format so calculating difference in days is not a problem(and I also know that the usual interval for Landsat 8 scenes is 16 days, but my list doesn't have all rasters from that year so the difference between 2 time adjacent rasters in a list can be for example, 32 instead of 16 days).

So with that in mind I'd like to get another band as an output where I would se the difference in days for each delta pixel.

As you can see in the image, the target raster is "input4" (or in the code above "post"). If "input3"(or "pre") has actual value in its cell then output cell equals input4-input3, but if input3 has NoData (ND) in a cell then output cell equals input4-input2 (unless input2 also has NoData) and so on.

if the target raster in question was "input3" then it would only be compared to input2 and input1 and first raster in a list doesn't get compared to anything obviously. If there would be cloud on the same spot in every raster(quite unlikely) then final output would be NoData. And if the target raster itself has a cloud then output is also NoData.

Band 2 equals to difference in days for the calculated delta values (and in this case the time cycle between two Landsat 8 scenes is 16 days).

Can this be done with ArcPy script?


texts: Importing gdal in Python with anaconda changes environment

texts:
I have installed gdal through anaconda. When I run a subprocess in Python before importing gdal I can get that gdal tools are located in my conda environment. For example subprocess.run(['which', 'gdal_rasterize']) would show /path/to/env/bin//gdal_rasterize

This is good because whenever I run gdal_rasterize within a subprocess it runs within the environment, however as soon as I import gdal into my Python interpreter or script (which I use gdal for other uses) whenever I run the same subprocess.run(['which' , 'gdal_rasterize']) I now get /usr/bin/gdal_rasterize

In fact importing gdal also changes the path to which version of Python I'm using in the environment. subprocess.run(['which', 'python']) changes to /usr/bin/python

This all wouldn't be an issue but I'm working on a server with tools that are less up to date than the conda environment and certain gdal functions have different parameters. Just trying to figure out why the environment would be changed for some reason.

Does anyone have any ideas?

I am on Linux (CentOS 6.6)


texts: I want to open an Esri .gdb with QGIS. I have adjusted GDAL and it works fine with an example .gdb provided on another thread on the topic. I still cannot open the .gdb of interest to me and suspect it is because it may not be a 10.something gdb, but probably an older .gdb file version. I have currently no access to Esri software. Anybody out there knowing a solution to the problem?


texts: Is it possible to add text and numbers into the same column/field of attributes of a shapefile in ArcGIS 10.3?

texts:
I have to add a field to the attributes of a shapefile and the data reads as:

2
BDL
22
18
11
BDL
BDL
25
3

Could these data be added to the same field/column of the attributes of the shapefile in arcmap 10.3?


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

texts:
My issue is closely related to:

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

and somewhat similar to:

Why isn't QGIS detecting CRS from .prj file?

but it seems that the root of the problem is different.

I am downloaded a shapefile (and all its necessary additional files) from a government repository. When I load the shapefile into QGIS, the CRS is not automatically detected even though a .prj file is in the folder.

How can I tell QGIS : "please look here and use this .prj file."

DETAILS:

Despite indicating that they use a standard projection (Canada_Albers_Equal_Area_Conic, ESPG: 102001), it actually seems to be a custom projection with a slight tweak to that standard so I can't just choose a CRS from the menu. [If I select "Canada Albers Equal Area Conic", then the polygons are displaced relative to where they should be based on other files for which I know the projection is correct]

The .prj file looks like this:

PROJCS["Canada_Albers_Equal_Area_Conic",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",
SPHEROID["GRS_1980",6378137.0,298.257222101]],
    PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],
    PROJECTION["Albers"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-95.0],
    PARAMETER["Standard_Parallel_1",49.0],PARAMETER["Standard_Parallel_2",77.0],PARAMETER["Latitude_Of_Origin",49.0],UNIT["Meter",1.0]]

When I load the file in R, the projection is automatically detected correctly:

> crs(test)
CRS arguments:
 +proj=aea +lat_1=49 +lat_2=77 +lat_0=49 +lon_0=-95 +x_0=0 +y_0=0
+datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0

To compare with standard Albers:

CRS arguments:
 +proj=aea +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0
+datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0

If I compare the output provided by the http://epsg.io/102001 website (suggested in an answer), I get the following for the OGC WKT (standard way of writing projections?) for the standard Canada Albers projection:

PROJCS["Canada_Albers_Equal_Area_Conic",GEOGCS["GCS_North_American_1983",DATUM["North_American_Datum_1983",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["False_Easting",0],PARAMETER["False_Northing",0],PARAMETER["longitude_of_center",-96],PARAMETER["Standard_Parallel_1",50],PARAMETER["Standard_Parallel_2",70],PARAMETER["latitude_of_center",40],UNIT["Meter",1],AUTHORITY["EPSG","102001"]]

And the following as the ESRI WKT way of writing the projection:

PROJCS["Canada_Albers_Equal_Area_Conic",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Albers"],PARAMETER["False_Easting",0],PARAMETER["False_Northing",0],PARAMETER["central_meridian",-96],PARAMETER["Standard_Parallel_1",50],PARAMETER["Standard_Parallel_2",70],PARAMETER["latitude_of_origin",40],UNIT["Meter",1]]

I am noticing slight differences in the language/terms used to set the projection.


texts:
I have column that contains dates that are created from a numpy array to fc, since numpy doesn't handle date types I have had to create another field to update the rows.  When using this script in the ArcPy window, no errors are returned, however my values do not change for UPDATE_DATE_DATE.  What could cause this?

UPDATE_DATE_DATE is a DATE type field while OPEN_TIME is a string type field.  And I would like to set the string values as DATE type in UPDATE_DATE_DATE.

Here is an example of the date format for OPEN_TIME, 9/25/2015 12:29:04 AM.

cursor = arcpy.UpdateCursor(fc)
for row in cursor:
    row.setValue("UPDATE_DATE_DATE", row.getValue("OPEN_TIME"))
cursor.updateRow(row)


texts:
I am using the Esri proxy to avoid CORS issues as I am testing against clients data held on another domain.

The wierd issue I am having is that WFS data only seems to work via one route, whereas mapping data seems to need another.

I have hosted the proxy locally using IIS. The following works for the WFS layer (ie using esriConfig.defaults.io.proxyUrl)

esriConfig.defaults.io.proxyUrl = "/ESRI-Proxy/proxy.ashx";

    // *************************************************************
    // Add a WFS Layer
    // *************************************************************
    var url = "https://clientsServer/arcgis/services/TEST/MyMapService/MapServer/WFSServer";

    var opts = {
        "url": url,
        "version": "1.0.0",
        "name": "GGS_Existing_Cables",
        "wkid": 4326,
        "maxFeatures": 1000
      };

    var layer = new WFSLayer();
    layer.fromJson(opts);

    var map = new Map(getMapDivId(), {
        basemap: "gray",
        sliderorientation : "horizontal"
    });

    map.addLayer(layer);

However when I then bring in the mapping layers I get CORS warnings (including for the WFS layer).

The mapping code is as follows:

// *****************************************************************
    // * Create two ArcGISDynamicMapServiceLayer layers.
    // *
    // * One shows the OpenStreetMap layer as a basemap. The other loads
    // *   pipeline map data off the west coast of Africa from iTechs
    // *   ArcGIS Server Instance.
    // *****************************************************************
    var imageParameters = new ImageParameters();
    imageParameters.format = "jpeg"; //set the image type to PNG24, note default is PNG8.

    //Takes a URL to a non cached map service.
    var westAfricaPipelineLayer = new ArcGISDynamicMapServiceLayer(geturlToWestAfricaPipelines(), {
     "opacity" : 0.5,
     "imageParameters" : imageParameters
    });

    var streetMapBaseLayer = new ArcGISDynamicMapServiceLayer("http://server.arcgisonline.com/arcgis/rest/services/world_street_map/mapserver", {
     "opacity" : 0.5,
     "imageparameters" : imageParameters
    });

    map.addLayer(westAfricaPipelineLayer);
    map.addLayer(streetMapBaseLayer);

     // *****************************************************************
     // Add the Nile Delta WMS Layer
     // *****************************************************************
     var layer1 = new WMSLayerInfo({
         name: '1',
         title: 'Nile Delta'
     });

    var resourceInfo = {
        extent: new Extent(-126.40869140625, 31.025390625, -109.66552734375, 41.5283203125, {
        wkid: 4326
        }),
        layerInfos: [layer1]
    };

    var wmsLayer = new WMSLayer(geturlToNileDeltaWmsMap(), {
        resourceInfo: resourceInfo,
        visibleLayers: ['1']
    });

    map.addLayers([wmsLayer]);

However if I add the following, the mapping layers work but the WFS layer still stops working:

urlUtils.addProxyRule({
      urlPrefix: getiTechRestfulBaseUrl(),
      proxyUrl: getUrlToProxy()
    });

    urlUtils.addProxyRule({
      urlPrefix: getiTechWMSBaseUrl(),
      proxyUrl: getUrlToProxy()
    });

getUrlToProxy is as follows

const URL_TO_PROXY = "http://localhost/ESRI-Proxy/proxy.ashx";

I would have thought that they amount to the same thing, i.e. esriConfig.defaults.io.proxyUrl and urlUtils.addProxyRule?

I have either (hopefully not) found a bug, i'm doing something very stupid, or I am using the two proxy routes in an unsafe manner (or a combinataion of...).


texts: Save Raster Layer As... not working in QGIS 3.4

texts:
I am trying to save a VRT layer as a TIFF file using the layers panel in QGIS 3.4.

I successfully created and saved the VRT from multiple ASC LiDAR tiles. My method is to work from the layers panel, right-click the VRT layer and select Export/Save As. The Save Raster Layer As... window opens, but whatever I enter -- raw data or rendered image, or any other options -- nothing happens and no file is created.

There are no error messages. The Cancel button does not work: the window stays open. I can close the window using the Close Window 'x' in the bar.

I can convert the VRT to a TIF if I use Raster/Conversion/Translate but that is not so intuitive, and the option appears to be there in the layers window.

QGIS 3.4.9 LTR, Linux Xubuntu 18.04 LTR, 64-bit


texts:
I have a custom function that pulls the markers from a json file called markers.json :

var markers = {
  "resources": [
    {
      "name": "AITokarServer",
      "coords": [
        {
          "x": -1210.0,
          "y": -1770.0
        },
     },
    {
      "name": "AIShigiServer",
      "coords": [
        {
          "x": -1310.0,
          "y": -1570.0
        },
     }
}

This is the function:

  function showResourcesByName(resName) {
    for (var i = 0; i < markers.resources.length; i++) {
        for (var j = 0; j < markers.resources[i].coords.length; j++) {
            var name = markers.resources[i].name;

            if (name == resName) {
                var x = markers.resources[i].coords[j].x;
                var y = markers.resources[i].coords[j].y;

                L.marker([y, x], {icon: bor}).addTo(map).bindPopup(name);
            }
        }
    }
  }

Using the Function to show markers on map:

  showResourcesByName("AIBorServer");

How can I change the icons of each "name" field after I call them on the function?

Something like:

showResourcesByName("AIBorServer"),(BorIconHere);
showResourcesByName("AIShigiServer"),(ShigiIconHere);


texts: Get a gps navigatable file from a Google map available on a website

texts:
I am trying to download this - http://googlemapkml.publishpath.com/Websites/googlemapkml/files/Content/5880032/GEDriving66.htm - map on to my phone. My end goal is to be able to drive on this route, so it would be great if this map could be opened on my phone (Android) and I could get GPS-assisted navigation.

I am totally new to these things.


texts: Loading big GeoJSON polygon data from GeoServer using OpenLayer bbox loading strategy

texts:
This is my project as beginner to learn how to use OpenLayers and GeoServer.
I'm using GeoServer version 2.19.0 to share my shapefile with so many polygon in it. Using EPSG:4326 as the projection, I can show it all using OpenLayers 6.5 without much hassle. But, showing all the polygon at once to user is not good, so I try to make filter using bbox loading strategy, hoping it would render only visible layer which inside of current map view.
But my code doesn't work like the way I needed, the GeoJSON vector doesn't rendered based on the view extent.
At first I think if the projection was the culprit, so I transform the coordinates of 'extent' to EPSG:4326, since my shapefile used it as the projection.
    var extentProjected = ol.proj.transformExtent(extent, 'EPSG:3857', 'EPSG:4326')

And then using bbox feature in GeoServer WFS service, I send the query.
        'http://localhost:8080/geoserver/ows?service=WFS&version=1.1.0&' +
        'request=getFeature&typeName=katingan:Indo_Desa_region&maxFeatures=3000&' +
        'outputFormat=application/json&' +
        'bbox=' + extentProjected.join(',')

It doesn't work, GeoServer sent blank features of GeoJSON.
After that, I change the query and set explicitly the SRS in it to this:
var url = 'http://localhost:8080/geoserver/ows?service=WFS&version=1.1.0&' +
        'request=getFeature&typeName=katingan:Indo_Desa_region&maxFeatures=3000&' +
        'outputFormat=application/json&srsName=EPSG:4326&' +
        'bbox=' + extent.join(',') + ',EPSG:3857'

It work, the layer shown but can't be filtered when I zoom to specific area of map.
Here's my shapefile: download shapefile
and Here's my code:
var vectorSource = new ol.source.Vector({
    format: new ol.format.GeoJSON(),
    url: function(extent) {
        // var extentProjected = ol.proj.transformExtent(extent, 'EPSG:3857', 'EPSG:4326')

        var url = 'http://localhost:8080/geoserver/ows?service=WFS&version=1.1.0&' +
            'request=getFeature&typeName=katingan:Indo_Desa_region&maxFeatures=3000&' +
            'outputFormat=application/json&srsName=EPSG:4326&' +
            'bbox=' + extent.join(',') + ',EPSG:3857'
            console.log(extent, url)
        return (
            url
        )
    },
    strategy: ol.loadingstrategy.bbox
})
var vector = new ol.layer.Vector ({
    source: vectorSource,
    style: new ol.style.Style({
        stroke: new ol.style.Stroke({
            color: 'rgba(0, 0, 255, 1.0)',
            width: 2,
        })
    })
})
var map = new ol.Map({
    target: 'map',
    layers: [
      new ol.layer.Tile({
        source: new ol.source.OSM()
      }),
      vector
    ],
    view: new ol.View({
      center: ol.proj.fromLonLat([113.91794724978148, -2.217485671684518]),
      zoom: 4
    })
});

Do you have any clues for me, so the map rendered faster, and it doesn't always show all the GeoJSON layer at once?


texts: r.thin (GRASS) not giving correct result

texts:
I have raster image of road feature. I want to extract line feature from it. So I run r.thin tool on it of parameter -
{ 'GRASS_RASTER_FORMAT_META' : '', 'GRASS_RASTER_FORMAT_OPT' : '', 'GRASS_REGION_CELLSIZE_PARAMETER' : 0, 'GRASS_REGION_PARAMETER' : None, 'input' : 'C:/Users/User/Music/test111.tif', 'iterations' : 200, 'output' : 'TEMPORARY_OUTPUT' }

and then r.to_vect. But the result of r.thin is giving one or two yellow pixels.

What is wrong I am doing?


texts: Determining units of SHAPE_Length and SHAPE_Area fields and reducing/rounding those values?

texts: When I run the INVEST SDR I got an error as "Possibly due to too larger number with respect to field width" that causes by these large integers of the shapefile area data that is shown as below in attribute table, does anyone knows a solution for this and how to know the units of these data


texts: What is Kappa coefficient in HSI classification?

texts: What is the Kappa coefficient and how it is calculated in the HSI classification process?


texts: Does anyone know why layers would continually be lost from a project and refuse to reload even though the project is regularly saved?


texts: Open Source solution to batch process shortest distance between points along road networks?

texts:
I have the locations geocoded (a fairly small dataset, ~200), but I haven't found a straightforward way to load these data into QGIS and perform network analysis between nodes.

Any suggestions?

I'm aware of the Road graph plugin, but I don't know how to load a batch of data.

I've also tried to figure out how to do this with TOPOJSON and d3.js without much success.


texts:
I need to parse a KML file which contains the line already shown in the title. I know it's an integer value and I'm currently storing it in code as private final int mForrestIdGisSu; but I'd like to know what this value is if at all possible. FID I think is "Forest I.D." and GIS might be realated to the GIS software, but I'm not sure about SU. Has anyone ever seen this before?

P.S. The KML itself is used to describe the production boundaries of different parts of a forest.


texts: Is there a setting or settings in QGIS 3 that allows users to modify application tool-tip behavior?

texts:
I can't read the longer tool-tip text in QGIS because the tooltop text is longer than the width and it doesn't wrap past 2 lines:

Is there a setting or configuration change to either widen or heighten the tool-tip?

In QGIS 2.18, the same tool tip doesn't wrap and is wide enough to read:

High probability that this is a bug, but I wanted to see if this is resolvable via a setting.


texts: Append (only) new features to a layer in QGIS

texts:
I've got a vector layer Target. I want to periodically add features that arrive (with the same attributes) in new.csv, which was created in the field.

I'd like to append them to Target if and only if they're not already there. I don't want to duplicate existing features, and I want to keep any tweaks done in Target, i.e. I don't want to replace existing features from new.csv either. (The features have a unique and consistent ID in a Name field.)

What's the simplest way to do this, in a way that can easily be (semi)automated to run weekly? I can think of the following, but seeking advice before I reinvent the wheel:

Use the merge processing algorithm and then postprocess to delete duplicates, making new.csv have a flag in a each feature's comment field to know which duplicate to nuke.
Dust off my rusty SQL and run some sort of APPEND query in DB Manager, with an EXCEPT clause to skip what's already in Target?
Write a PyQGIS script to traverse the features in new.csv, throw out those where Name is already in Target, and then use layer.dataProvider().addFeatures(featureList) to add what's left to Target.

Suggestions? I don't think it matters, but the feature geometry are points in EPSG:4326. Target is one of several layers in a gpkg, and new.csv is of course a separate comma separated file, which can be imported as a layer with point geometry.


texts: I'm using a text box in the print composer (QGIS 2.16).  If I write more then one line, how can I set the line spacing between the lines?


texts: Check if circle is within polygon/rectangle/circle in leaflet and turfjs

texts:
I am using leaflet and turfjs. I am able to check whether rectangle/polygon/point is within another polygon.
But I am not able to check the same for circle. I want to restrict user to draw circle inside other shape (polygon/circle/shape) etc.

Could anybody please help me to know how to do this?

This is how I check whether new shape is within polygon or not. It is working fine when new shape is rectangle or polygon. But for circle, it throws errors at  for "layer.getLatLngs()".

// create parent feature
var parentContainerBoundary = {
    "type": "FeatureCollection",
    "features": [
        {
            "type": "Feature",
            "properties": {},
            "geometry": {
                "type": "Polygon",
                "coordinates": [parentCoordinates]
            }
        }
    ]
};

        var pointWithinContainer = turf.within(newShapePoints, parentContainerBoundary);
        return pointWithinContainer.features.length != newShapePoints.features.length;


texts: Unable to integrate CartoDB into Leaflet Map

texts:
This is my first day with CartoDB, trying to integrate our visualization into our existing Leaflet implementation, and having issues. In short: no visualization.

Below is a sample of the problem followed by the code involved. I have a createVis test that works fine, followed by a createLayer test that's not working. The layer seems to add correctly but the tiles don't display.

I've tried:

With and without lower layers
Adding the layer in the chain and
explicitly with "layer.addTo(map)"
Explicitly making it visible with
"layer.show()"

No luck so far, and running out of ideas.

http://www.zonability.com/carto.html

        var url = 'http://zonability.cartodb.com/api/v2/viz/e711de5e-8f77-11e4-a96a-0e0c41326911/viz.json';

        // CreateVis
        cartodb.createVis('map1',url);

        // Integrate into existing Leaflet map
        var map = new L.Map('map2',{
            center: [30.242713340145,-97.770476601062],
            zoom: 9
        });
        map.on('layeradd',function(layer){ console.log('Layer added'); });
        L.tileLayer('http://{s}.tiles.mapbox.com/v3/zonability.glj45k32/{z}/{x}/{y}.png').addTo(map);
        L.tileLayer('http://{s}.tiles.mapbox.com/v3/zonability.gljpogdp/{z}/{x}/{y}.png').addTo(map);
        cartodb.createLayer(map,url).on('done',function(layer){
            layer.addTo(map);
            console.log('Done');
        }).on('error',function(err){
            console.log('Error: '+err);
        });


texts: Downloaded shapefiles become 1KB when I try to upload to QGIS

texts:
This problem is happening when downloading multiple GIS files from different sites.
I downloaded 2012.zip from https://www.sos.wa.gov/elections/research/precinct-shapefiles.aspx to see the precinct maps for 2012 in Washington state.

Everything looks normal, so I extract the data from the zipfile so I can upload it into QGIS with no problems. I go to add a shapefile to my new project, and select the shp file.

I leave everything in default, and go to add the shapefile.

Once I've done so, I see that a layer name "Statewide_Prec_2012" has been added but nothing shows up. The area where the map should be is completely blank, despite the layer being selected. I "zoom to layer" and the screen stays the same. I open the attributes table and it's blank. Now, if I go and look at the extracted data from earlier, this is what I find.

Does anyone have any clue what is happening? I never ran into this issue using ArcGIS but switched to using QGIS so I can work on projects on both of my computers.


texts: When I add a feature to a vector layer, all the features of the layer disappear Openlayers

texts:
I'm working in a Geocoder tool that can put new features in a vector layer directly just pressing a button. When press the button all the features of
the layer disappear in the map. The features still exists because I can get
the length increased by one using after that:

 alert(vector.getSource().getFeatures().length)

I did the following function as button click handler

function addPoint(lng,lat)
{
  var pp = new ol.geom.Point( [
    ol.proj.transform([lng,lat], 'EPSG:4326', 'EPSG:3857'),
]);
var featurething = new ol.Feature({
   geometry: pp
});
vector.getSource().addFeature(featurething);
}


texts:
In QGIS I have two  vector layer shapefiles which should overlay approximately one on the other.  They have both been downloaded from the same government website (data.vic.gov.au) with the same projection (GDA94 Zone 55).  When loaded into QGIS they show the same CRS (EPSG:28355) which is also the project CRS, "on the fly" projection is on,  the content of their .prj and .qpj files is identical.  However when loaded they appear on the screen many kilometres apart.

What is going on?


texts: Automating map creation for 50,000 unique routes using ArcGIS Desktop?

texts:
I have a table with 50k entries. There is a to field and from field with location in either lat/lng or street address. All these location are in the USA

I would now like to automate the creation of unique jpg maps for each entry showing the route to get from point a to b. I have a ArcGIS 10 license along with spatial analyst (not sure if these will help).

Any suggested workflows to automate the creation of these maps?


texts: Why can't I perform any operation over a 700MB MIF file containing 5 multipolygons? It takes forever

texts:
I'm trying to work with a MIF file containing 5 really big multipolygons, but every operation I try takes hours, or just crash because of memory issues. Mainly using ogr2ogr, to transform MIF file to other format, like SHP or PostGIS. Unsuccessful, so far.

The file weights 700MB. I could load it into PostGIS using ogr2ogr. I just left the command running all night and it finished silently. The size of the table in disk is about 160MB. So, I suspect something went wrong (looks too small)

UPDATE

This is the header of MIF file

Version 300
Charset "Neutral"
Delimiter ","
CoordSys Earth Projection 8, 28, "m", -3, 0, 0.9996, 500000, 0
Columns 2
  Float Float
  String Char(32)
Data

After that, 5 regions with:

124222 polygons
462270 polygons
798815 polygons
1124915 polygons
1706436 polygons

About the PostGIS table, I'm trying to check it. It contains 5 rows, like the original MIF file, but every single operation I try to execute over this table using the geometry column also takes hours. I tried to add a new integer field, and after 4 hours is still running. With ps, the process says "ALTER TABLE waiting".

This is the schema of the table

    Column    |       Type       |                         Modifiers
--------------+------------------+----------------------------------------------
 ogc_fid      | integer          | not null default nextval('my_table_ogc_fid_seq'::regclass)
 wkb_geometry | geometry         |
 float        | double precision |
 string       | character(32)    |
Indexes:
    "my_table_pk" PRIMARY KEY, btree (ogc_fid)
    "my_table_geom_idx" gist (wkb_geometry)
Check constraints:
    "enforce_dims_wkb_geometry" CHECK (st_ndims(wkb_geometry) = 2)
    "enforce_srid_wkb_geometry" CHECK (st_srid(wkb_geometry) = 23030)

So, how could I handle this big MIF file? I can't do anything with it, neither with the result PostGIS table.


texts:
I'm having an issue with QGIS version 3.2 Bonn. For my class, a question asks, "What city corresponds to row number 1713 for the Colleges and Universities shapefile?". I'm having issues answering this, because each time I load the shapefile into QGIS, the attribute table order is randomized, giving me a different city each time.

How can I open the "default" version of the attribute table for this layer, so that the order of rows is not randomized each time?


texts: How to calculate the average value of a raster within polygons?

texts:
I am using the GRASS plugin in QGIS.  I am trying to calculate the average value of a raster dataset within a set of polygons.

I have created a new mapset and have imported both the vector and raster layers into the GRASS mapset.

I am now trying to calculate the set of statistics using the v.rast.stats tool.

It starts to run and then comes up with the error an error occurred while converting vector to raster and then 'raster map not found'.

v.rast.stats raster=secant_grass@zonalstats
vector=plots_grass@zonalstats
colprefix=g/Applications/GRASS-6.4.app/Contents/MacOS/scripts/v.rast.stats:
line 190: 632 Segmentation fault v.to.rast in="$VECTORFULL"
out="${VECTOR}_$TMPNAME" layer="$GIS_OPT_LAYER" use=cat type=area
--quiet ￼An error occurred while converting vector to raster ￼Raster map not found ￼ nothing removed Finished with error

Any help would be greatly appreciated.

Regards


texts: GeoServer CQL Filter doesnt work with Unicode

texts:
I'm using CQL filter to extract certain data from my GeoServer layer.
var source = new ol.source.Vector({
        loader: function (extent, resolution, projection) {
            var url = "http://localhost/geoserver/visualization/ows?service=WFS"
                + "&version=1.1.0&request=GetFeature"
                + "&typeName=visualization:".concat(inputTable)
                + '&outputFormat=text/javascript'
                + "&format_options=callback:loadFeatures"
                + '&srsname=EPSG:3857'
                + filter;
            $.ajax({ url: url, dataType: 'jsonp', jsonp: false });
            console.log(filter);
            console.log(url);
        }
    });

Notice the filter string as the final parameter at url. If I normally input the feature name using alphabet text, it will return the result. But if I use other than alphabet (UTF-8, EUC-KR) it won't return anything.
&cql_filter=branchCode='N01821'  -> works
&cql_filter=방향='S'             -> korean word, doesn't work

If I print the URL, the Korean strings give me this error
<ows:ExceptionText>Could not parse CQL filter list. Lexical error at line 1, column 1. Encountered: "\ubc29" (48169), after : "" Parsing : 방향='S'.</ows:ExceptionText>

Anyone knows the reason why and the solution?


texts:
I am working on a project based on QGIS API c++ 2.8.2 and qt4.8.

I have two methods in a class, the first create a layer and add a point
Here is the code for the creating:

LibCartoError CMyQgsS57::createPointLayer(QString pLayerName, QMap<QString,
        QVariant> pFeatureName)
{
    QgsVectorLayer * lVectorLayer = new QgsVectorLayer("Point?crs=epsg:4326", pLayerName, "memory");

    lVectorLayer->startEditing();
    QgsVectorDataProvider * lDataProvider = lVectorLayer->dataProvider();

    QMapIterator<QString, QVariant> lAttrib( pFeatureName );
    QList<QgsField> lMyAttribField;

    while (lAttrib.hasNext()) {
        lAttrib.next();
        // add attribute
        lMyAttribField << QgsField(lAttrib.key(), lAttrib.value().type());
    }
    lDataProvider->addAttributes( lMyAttribField );

    // save attributes
    lVectorLayer->commitChanges();

    // add to the map
    QgsMapLayerRegistry::instance()->addMapLayer( lVectorLayer );

    mMapLayer.insert(pLayerName, lVectorLayer);

    ((QList<QgsMapCanvasLayer>*)myListForPrint)->insert( 0, QgsMapCanvasLayer(lVectorLayer, TRUE) );
    // return
    lError.error_code = NoError;
    lError.error_detail = "La creation de la couche Point " + pLayerName + " est effectuée";
    return lError;
}

I have a timer that delete this point if the result of an operation is good, else the point is kept in the map (or canvas).

Here is the code for deleting  the layer:

LibCartoError CMyQgsS57::deleteLayer(QString pLayerName)
{
    LibCartoError lError;

    int lIndexOfLayer;
    bool lLayerFound = true;
    for (int k = 0; (k < ((QList<QgsMapCanvasLayer>*)myListForPrint)->size() && lLayerFound); k++)
    {
         if(((QList<QgsMapCanvasLayer>*)myListForPrint)->at(k).layer()->name() == pLayerName) {
            lIndexOfLayer = k;
            lLayerFound    = false;
        }
    }

    if( !lLayerFound )  {
        // delete the layer
        QgsVectorLayer* lLayertoDelete = NULL;
        lLayertoDelete = (QgsVectorLayer*) mMapLayer.value( pLayerName );

        QgsMapLayerRegistry::instance()->removeMapLayer( lLayertoDelete->id() );
        mMapLayer.remove( pLayerName );

        ((QList<QgsMapCanvasLayer>*)myListForPrint)->removeAt(lIndexOfLayer);
    }

    lError.error_code = NoError;
    lError.error_detail = "Layer :: " + pLayerName + " is deleted";
    return lError;
}

If I launch the system monitor, my application takes 4.4 MB each time the delete and create is invoked. If only the createpoint is invoked the use of memory is stable, but if I create / delete  a lot of time I can see the use memory increasing.

Is there something missing in my code?


texts: PostGIS Tiger Geocoder Not returning, stno, street, styp for geocoding

texts:
I recently loaded TIGER 2018 data for California. After fixing all the errors I have loaded all the data into the database. I checked each of the tables and all of them contain data.

From this example
How to make a PostGIS TIGER Geocoder in Less than 5 Days
I ran this query to check.

SELECT g.rating,
    ST_AsText(ST_SnapToGrid(g.geomout,0.00001)) As wktlonlat,
    (addy).address As stno,
    (addy).streetname As street,
    (addy).streettypeabbrev As styp,
    (addy).location As city,
    (addy).stateabbrev As st,(addy).zip
FROM geocode('424 3rd St, Davis, CA 95616',1) As g;

I get this result.

 rating |         wktlonlat          | stno | street | styp |    city     | st |  zip

 100    | POINT(-120.91795 38.29745) |      |        |      | Buena Vista | CA | 94502

As you can see the stno, street and styp are empty. I am not sure what the problem is as all the data has been loaded.


texts: Update map tiles in Leaflet only when map panning ends?

texts:
Currently when the map is dragged a few hundred API requests are sent due to many layers.

I would like the requests to only be sent when the map comes to a stop.

How can I achieve this?


texts:
I am trying to create a certain type of label common in the UK forestry sector. It consists of a circle with a horizontal dividing line, with one attribute in the upper portion and a second in the lower. I have concatenated the two parts of the label, with a newline character, and giving the label a circular border, but I can't work out how to add the horizontal dividing line.

I'm using the most recent edition of QGIS, Essen, in 64-bit windows.


texts: In which format am I working and how can I plot it correctly on GIS?

texts:
I've been working recently replicating this paper: "Constructing Street Networks from GPS trajectories". They got a repository with the implemented code. Anyway, the issue is that their data is in the following format:

391821.292854,5821122.049075

Reading carefully the code and the paper I have figured out that these coordinates are in meters after Mercator projection and are probably in ESRI format.

The problem that I have is that I cannot figure it out how to plot them correctly on QGIS. The file that I have is on txt. I'm using the Delimited Text plug in to import the data and playing with the CRS configuration without results. The data appears nowhere near Berlin (where the data is from)


texts: How to count the number of polygons within a grid?

texts:
I have two vector layers: 1) with polygons representing landcover types and 2) a grid, overlapped with the first one.

I want to obtain, for each square within the grid, the number of land cover types. For example:

In the red square it would have 2 landcover types (2 and 7), and the left square to that one only 1 land cover type (2).

How do I obtain this values, using QGIS?


texts: Can you print cross-hatch area fills like they appear on screen?

texts:
I'm trying to use a cross-hatch symbol on polygons in QGIS.

When displayed on screen, they look just fine - both in the map and in the composer.

The problem is when they are rendered to an image and/or printed, the hatching is much smaller in size (between the hatches) than on screen. The result is a much darker looking polygon that simply looks black, rather than a cross-hatch pattern that can be seen through to the colours below.

There is a print resolution setting (DPI) that when changed to 96dpi, it will render as it does on screen. But how to fix for 300dpi printed maps?

Thanks!


texts:
In OpenLayers 5, I have some LineString that represent the streets. Now I want that when the map is zoomed also the LineString should be zooomed.

For example in the following image, I have two linestrings at zoom level 1. At zoom level 2 near the blue lane, this should be zoomed and the red lane should be horizontally scrollable. How can I do this?


texts: I am trying to build a model that reads through about 30 different point feature classes and appends only the selected features to a new feature class. Some of the feature classes don't have any selected records, so it is appending every record from that feature class. Is there a way to exclude these features from being appended? In the past we used a VBA macro to read through each feature class and uses selection sets to determine if any of the records are selected, but I can't figure out how to get that same functionality using Python.


texts: The world map that I am using from NaturalEarth.com has France and French Guyana as one feature.  I want to split these into two separate shapes. I'm sure there is a quick fix for this, and the others I'm sure to soon encounter. I am using ArcMap 10.2 and need to use a feature that will make it so that when I click on Guyana, France isn't also selected.


texts: Understanding projection, extent and scale in QGIS?

texts:
I am creating a map using QGIS that has shapefiles in different projections. I was able to reproject them to the correct ESPG code that I need, however they do not line up.  My base layer "illinois state plane zones" has the extents;
North: 42.508347916666146
West: -91.51351833333352
South: 36.969971916666374
East: -87.49521391666703

Scale: 119:1

The layer that is not lining up correctly has extents in different units, so I can understand why they will not overlay each other, but how can I rectify this?

I have re-projected and that doesn't seem to help.

Layer 02062019 Extents are:
North: 1801408.0041401046
West: 1196576.2196446825
South: 1799584.9969252362
East: 1198603.4151709464

Scale:1:2543

Also the scales don't match.

Can I georeference a shapefile to another shapefile to correct this issue, at least as far as the scale issue is concerned?

Re-projecting didn't seem to help at all and I need to figure this out.

What could be causing the extent problem?


texts: How to add a new database connection to QGIS with Python?

texts:
I have a plugin for QGIS that helps to populate databases.  The user has the option to create a new database or to use an existing connection.  If the user chooses to make a new database, I would also like to add this database as a new connection in the browser panel :

I can find lots of help on connecting to the database, which I have been able to do, but nothing for how to add a new connection.  Is this possible?


texts:
I'm following this tutorial about Mapbox GL and flask which has been going great until the step to add markers to the map. I've tried to distill it to make it as simple as possible but it's still not working.

I'm a JavaScript newbie and at a loss to the best way to debug a JavaScript map.

Does anyone have ideas on best way to debug JavaScript, or can anone see any obvious errors going on here?

        <meta charset='utf-8' />
        <title></title>
        <meta name='viewport' content='initial-scale=1,maximum-scale=1,user-scalable=no' />

        <link href='https://api.tiles.mapbox.com/mapbox-gl-js/v0.24.0/mapbox-gl.css' rel='stylesheet' />


texts: Easier calculation of a surface using QGIS

texts:
I'm looking for an easier way to calculate the surface on QGIS:

My project CRS is in Lambert93.

Currently, when I'm calculating a surface using QGIS, when the surface seems wrong, I always have to:

"Save as..." my current layer under a new CRS (Lambert93 in my case, (see Capture 1).
Then the new layer appears, and the $area calculation seems correct. (Capture 2)

My question is: this process doesn't seems very efficient as it creates new files on my computer; is there any other way to do this? (ie. without Save as... new layer, neither creating new files, ...)

Capture 1 : Save as the layer under a new CRS (Lambert93)

Capture 2 : my $area is now correct !


texts:
I am trying to add a legend in Leaflet based on a GeoJSON layer like this:

This is my code:
L.geoJson(Viales_Transitables).addTo(map);
var Viales_TransitablesStyle = {
    'weight': 1.3, // grosor de línea
    'color': "#ffffff",// color de linea
    'stroke':0.2,
    'opacity':10,
};
L.geoJson(Viales_Transitables, {
    style: Viales_TransitablesStyle
}).addTo(map);

function getColor(d) {
    return d === 1 ? '#fa0000' :
    d === 2 ? '#fa9b00' :
    d === 3 ? '#f6fa00' :
    d === 4 ? '#44fa00' :
    '#FFEDA0';
}

function style(feature) {
    return {
      color: getColor(feature.properties.Ancho_Tran),                                           weight: 2,
      opacity: 1,
      fillOpacity: 0.7
    };
  }

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

I would like to add this legend in Leaflet but I don't know how.


texts:
I'm trying to load a WFS layer in PyQGIS. My WFS uri is for county boundaries here in England. It's open data and pretty light, if you were wanting to recreate the issue.

Here's my code:

uri = 'http://www.geostore.com/OGC/OGCInterface?SERVICE=WFS&UID=UDATAGOV2011&PASSWORD=datagov2011&INTERFACE=ENVIRONMENTWFS&VERSION=2.0.0&LC=0'
layer = QgsVectorLayer(uri, "WFS_Layer", "WFS")
if not layer.isValid():
    print "Layer failed to load!"

Sure enough, the layer is not valid.

I'm following (to the best of my limited ability) the QGIS instructions here: http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/loadlayer.html

I'm losing the thread at this part:

Where in the declaration of vlayer is the uri being passed in? I feel I'm missing something fundamental.

Can anyone post working code? It'd be awesome to be able to get PyQGIS and WFS working nicely!


texts: Combine multiple binary rasters to a single raster

texts:
I have 5 different species richness maps which I am trying to join in order to show one cohesive map. Each species raster layer has values of 1 and 0; 1 if the species lives there and 0 if the species does not.

In other words, I want something like this:

if species 1==1,1
if species 2==1,2
if species 3==1,3

And so on until I reach 5. I want everything else to be classified as a 0.

I am very new to using raster calculator and am unsure how to set this up.

I just realized that this will not work since the species will overlap. I will need to figure something else out.


texts: I would like to use TomTom multiNet data with pgRouting to create drive time/distance polygons. I have gone through the documentation regarding network topology creation, but am still confused to proceed further. Could anybody help me out with the steps I need to take to make multiNet routable with pgRouting. Any pointers are welcome.


texts: Attempting thermal difference analysis in Google Earth Engine?

texts:
I am attempting to do a means of NDVIS and a means of Thermal thermal for Summer (July &August) 2016.

When I run the annotated code below the error I get in the side console panel is: Use print (...) to write to this console. And in the red box right below: "image" is not defined in this scope.

What's wrong with my code?

//specify a portion of Landsat archive to work with as image collection
var collection = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
  .filterDate('2016-07-01','2016-08-31')
  .filterBounds(ee.Geometry.Point(-95.32287567601213,29.79395929049292));

//add a function for NDVI calculation & apply it to the collection using collection.map():
var addNDVI = function(image) {
  return image.addBands(image.normalizedDifference(['B4', 'B3']));
};
var ndviCollection = collection.map(addNDVI);

//add a function for NDVI calculation & apply it to the collection using collection.map():
var addIR = function(image) {
  return image.addBands(image.normalizedDifference(['B10', 'B11']));
};
var ndviCollection = collection.map(addIR);

//normalize difference
var addNDVI = image.normalizedDifference(['B4', 'B3']);
var addNDVI = image.normalizedDifference(['B10', 'B11']);

//calculate mean value of each band from the collection as a new multi-band image
var mean = ndviCollection.reduce(ee.Reducer.mean());

//select the mean NDVI band and Thermal band from this reduced image to visualize and display:
var vizParams = {bands: ['nd_mean'], min: -0.5, max: 1, palette: ['FF0000', '00FF00']};
Map.setCenter(-95.32287567601213,29.79395929049292);
Map.addLayer(mean, vizParams, 'NDVI');

Map.addLayer(addIR,{min:0,max:700}, 'IR');


texts: Devising arbitrary coordinate system for site-plan that cannot be georeferenced?

texts:
In ArcGIS 10.1 (advanced lisence) I have imported the plan of a site. The plan has no coordinates, and I couldn't georeference the plan since it has undergone changes during the last 30 years. So, the map has no spatial reference. I used polygons to vectorialize some features I am interested in. Since my ultimate goal is to perform some spatial statistic analysis (e.g., calculating Moran's I for a numeric attribute stored for each feature), how can I devise an arbitrary coordinate system for the site in order to be able to use distance-based statistics?

As for earlier similar question: a somewhat similar question here (projection and coordinate system for non georeferenced archaeological site) got no replies. Further, another often-quoted thread of this forum (Using customized Coordinate System in ArcGIS Desktop?) does not address my issue. At the best of my understanding, what is suggested there entails that "...You have to know the geographic coordinates lat_0 and lon_0 of the origin of your CRS". This does not seem to apply to my issue since I do not know the coordinates of the origin of the coordinate system.


texts: Set symbol on lines according to line features length

texts:
I need to display arrows alongside to line features. Lines represent water pipes. The issue I'm facing is that lines are a concatenation of numerous line features of various lengths. Therefore, symbols are attached to every line, whatever its length. This makes symbols either overlaying or touching eachother (as in example below).

Inside my mapfile, I define the display rules of my symbols in the STYLE element:

    STYLE
        SYMBOL "wp_arrow"
        GAP -100
        COLOR 255 0 0
        OFFSET -5 -99
        SIZE 6
    END

I would like these rules to apply only to line features being more than 10 meters long (or x map units).


texts: How to batch create GeoServer stores for each user?

texts:
I have  unix system with a number of regular users, each with his/her own database in PostGIS. I need to create a PostGIS store for each user.

I know how to do this in the web interface.

Is it possible to create the stores in a batch file/script instead of doing it one by one?


texts: QGIS - tips for visualizing three overlapping layers?

texts:
Does anyone have a suggestion on how best to visualize three overlapping years as seen below?

The yellow/brown part is the base background color of the map.

The grey on the top of the screen represents adjacent countries that are greyed out

Then there are three added layers:

1) green on bottom with no transparency,

2) red on top with 85% transparency, and

3) blue on top with 85% transparency

I think it is difficult for people to visualize the three different layers. Any suggestions on how I could clean this up and make it easier to read? Perhaps on one layer use full transparency with only an outline? Or change the colors so it is easier to see what is overlapping with what? Or use polka dots or cross-hatches? Or do I simply have too many layers and I need to separate them into different maps?


texts: Using gdal to convert SInt16 DEM to UInt16 DEM, while setting nodata value to 0xFFFF

texts:
Can this be done on the command line?
I'm working with the Mars MOLA DEM GeoTIFF. Unfortunately, my framework can't deal with signed values, so it gets really confused when it looks at the data. I want to translate it to unsigned, and I want to adjust all the nodata values to 0xFFFFU.
The data set has values ranging from -8,500 to 21,226, and the nodata value is -32,768 (0x8000). I'd like to add 8500 to all the pixel values, change any nodata(0x8000) pixels to the new nodata value of (0xFFFF), and set the nodata value tag to 0xFFFF (65535).
I can't figure out how to get gdal to change any nodata values in the input dataset.


texts: How to calculate runoff in Rainfall Analysis?

texts:
I have 30 years  daily rainfall data. How can i calculate the mean monthly runoff, minimum monthly runoff and maximum monthly runoff with the available rainfall data using GIS?

I would like to use SCS curve number method for run off computation.


texts: Authentication problem: can't access proprietary ArcGIS portal

texts:
I am currently developping a MFC application that loads html and javascript code, using the ArcGIS JavaScript API (4.17), to display a map, where I am adding feature and graphic layers. I must as well add that I am not using the ArcGIS Web AppBuilder. All good until now and I want now to upload these data/feature layers on a specific proprietary ArcGIS portal. The authentication pattern I am trying to develop is the named user login workflow.
I tried thus to follow this tutorial and use the code provided by it: https://developers.arcgis.com/javascript/latest/sample-code/identity-oauth-basic/index.html
My application is registered on the proprietary ArcGIS account and I have its client/application ID.
As many, I have this "redirect_uri" error displayed when I am launching my application and try to sign in, and even when I am naively double click on the unmodified identity-oauth-basic.html file provided by the aforementioned tutorial. Following that, I got a bit confused with what I can find on the net to solve this problem.
From what I gathered from the net, I understand that I need to specify on the ArcGIS portal the URL of the page users of my app are redirected to after they successfully log in.
But in my case, everything is local and I am basically embedding an HTML file through my application. Thus I don't know what to specify as a "redirect url" in the "Authentication" pane of the registered application on the ArcGIS account.
I tried to add "http://localhost/" and "https://localhost/" as some tutorials and forum posts may specify it, but the problem is still here. I wouldn't actually know what I should type in a web browser to point to my HTML file, except from its absolute path.
I added as well as a "redirect url" the address I am specifying in my OAuthInfo instance. Here as well I am getting confused about the different specified URLs. The proprietary portal URL I am trying to work with has the following format: https://.maps.arcgis.com. But some URLs are specified with the '/arcgis' (redirects to a non-found page in my case) at the end and the aforementioned tutorial is even adding the '/sharing' extension (redirects to a page where I can be redirected to the ArcGIS API reference). What are the conventions here and what should I use in my OAuthInfo instance and on the application authentication settings?
Can someone clarify these different points so I can move forward with this authentication problem?


texts:
I have added a shapefile in QGIS using Python plugin. I want to add this shapefile when QGIS is loading so I have to add code in the initgui() method.

I have used following code in the initgui() method:

self.iface.addVectorLayer("C:/abc.shp","Global","ogr")

But when I open QGIS it asks me to save the file but not add/show the shapefile whether I saved  it or not.


texts: Using Polygon Neighbors for boundary polygons to find boundary lengths in ArcMap?

texts:
I am working in ArcMap. I have a map of districts in South Carolina, and I'm trying to find the boundary length that every district has with every other district, and the boundary of South Carolina itself. I understand that the GeoProcessing tool Polygon Neighbors can find the boundary lengths between polygons, but I also need the length that a district has with the edge of the state.

I am trying to verify that whatever values I get here correspond accurately with the perimeter of each district, so I don't simply want to use the perimeter tool and subtract off the sum of other neighbor lengths.

I have considered creating a buffer polygon around the state to act as an "outside" district, but I'm not sure how to do that.


texts: Using if-then with ArcMap raster calculator

texts:
I have several buffers to be added together (Cold, Hot, Raining, etc.), but they only need to be added if they exist.

So I currently have create buffer if Cold else no buffer. (I did convert polys to raster.)  I want to set up a calculator that runs all of the buffers that exist.  But I won't know which exist (large project with lots of weather specifications).

So far I am able to add them fine if everything is there, but if one specified buffer doesn't exist then it errors out.


texts: Filter function based on attribute. Leaflet

texts:
Im trying to create a webmap based on Leaflet - google spreedsheets example

I have a a "group" layer with three diferent group names.

I have added some easy buttons to click the layers on and off, but I have only been able to make it work on the entire "group" layer.

My buttons:

var buttons = {
  Agriculture: L.easyButton({
    position: 'topright',
    states: [{
      stateName: 'add-markers',
      icon: 'Agriculturee',
      title: 'add random markers',
      onClick: function(control) {
        map.addLayer(group);
        control.state('remove-markers');
      }
    }, {
      icon: 'Agriculturee',
      stateName: 'remove-markers',
      onClick: function(control) {
        map.removeLayer(group);
        control.state('add-markers');
      },
      title: 'remove markers'
    }]
  }),

so im using map.addlayer, but it should more be something like map.addlayer(group.Transportationhub)

what function should a add to my code or how should a define each of the names in the group as individual layers?

"Group" is defined as:

var group = L.featureGroup(markerArray);

The icons are styled different based on the values in the spreedsheet:

  function determineLayers(points) {
    var layerNamesFromSpreadsheet = [];
    var layers = {};
    for (var i in points) {
      var pointLayerNameFromSpreadsheet = points[i].Group;
      if (layerNamesFromSpreadsheet.indexOf(pointLayerNameFromSpreadsheet) === -1) {
        markerColors.push(
          points[i]['Marker Icon'].indexOf('.') > 0
          ? points[i]['Marker Icon']
          : points[i]['Marker Color']
        );
        layerNamesFromSpreadsheet.push(pointLayerNameFromSpreadsheet);
      }
    }

Log.console(layerNamesFromSpreadsheet) :

I have successfully used the easy buttons before with different filters based on the GeoJson, but this is different when the data is from spreedsheets.

Any simple funcktion a can use to seperate all the different names in the group?


texts: In Windows, multiple CPUs can be triggered in parallel, but in Linux, only one CPU can be used

texts:
When using QGIS, my environment is in docker (opensuse15.2)
Other answers are as follows:
Use this:
from PyQt5.QtCore import QSettings, QThread

# parallel rendering

QSettings().setValue("/qgis/parallel_ rendering", True)

for use max cores

from qgis.core import QgsApplication

threadcount = QThread.idealThreadCount ()

QgsApplication.setMaxThreads (threadcount)

It doesn't work in Linux


texts: Calculating Length of Linestrings in WGS84 in Miles

texts:
I have a few thousand LINESTRING features in a MySQL database. I'm trying to figure out the cumulative length of all the lines -- in miles. Using the GLength function I can retrieve the length in (what I assume to be) degrees, but I'm not sure how to get this to miles.

I also have QGIS and ArcMap, which both return the same degree values for the lines.

From what I've read the data needs to be projected, but I'm unsure what to project it to.

I'd like to perform the calculation within MySQL, if possible.


texts: Setup of development environment of QGIS on windows with qt creator

texts:
Can anyone help me in setting up environment for working on QGIS on windows using qt-creator?
I am using qt creator version 1.3.1,qt version 4.6 and QGIS version 1.8.

Thanks.


texts: What simple solutions are available for sharing big data with end users?

texts: I have approximately 1700 raster tiles showing classified land cover with a total size of approximately 300GB.  I am looking for a simple solution to privately serve these data to non-technical end users.  The end user needs to visually assess these tiled raster products.  Very little analysis would be required, only the ability to make the raster tiles transparent with a color map and overlay them on aerial imagery.  My initial thought was to create a gigantic mosaic of these tiles and serve the file via FTP.  Is there a better method to serve these data to the end-user, perhaps one that could elegantly handle tiled data?  I have very little web mapping experience, which should weigh into the solution. An 8TB, dual processor server and the ArcGIS family of products are available to me.


texts:
I am running QGIS 1.8 (SpatiaLite 3.0.1) and QSpatiaLite 6.0.1.

I have a spatial table EventPrecise (with geometry column p in EPSG:4326). I want to join this table with another, non-spatial table Event and produce a new spatial table out of it with EPSG:21781, while at the same time filtering on a non-spatial attribute t. This is my query:

SELECT e.id, e.oid, e.type, e.subtype, e.t, ST_Transform(ep.p, 21781) AS p
FROM event e LEFT OUTER JOIN eventprecise ep ON e.id = ep.eid
WHERE datetime(e.t) > datetime('2013-02-05 18:20:00') AND datetime(e.t) < datetime('2013-02-19 08:31:00');

When running this query in QSpatiaLite with options "Create Spatial Table & Load in QGIS", and geometry column = "p", I get the error:

'%s' doesn't correspond to a geometry field, try with another field name

Creating a View works though.. is this a bug?

EDIT: Turns out that it works when doing an INNER JOIN, so p is never null - still, is this normal behaviour?


texts: Arcpy - Create an Historical Marker on an archiving enabled dataset in Oracle SDE

texts:
I am looking for a way to create an Historical Marker programmatically with arcpy.

The archiving is enabled on a dataset in Oracle SDE.

I think, it's related with arcpy Version Toolset, but I can't find the documentation about this.

#this is an non-working attemp because CreateHistoricalMarker_management() method does not exists
import arcpy, datetime

my_dataset = "c:\\gis\\connection_file.sde\\an_archived_dataset\\a_feature_class"
marker_name = "A new dummy marker"
marker_datetime = datetime.datetime.utcnow()

arcpy.CreateHistoricalMarker_management(my_dataset, marker_name, marker_datetime)

>>>AttributeError: 'module' object has no attribute 'CreateHistoricalMarker_management'


texts:
When I enable an enterprise geodatabase for a database in postgresql, the st_geometry type get installed with it. So Is it necessary to use "Create Spatial Type" tool to create st_geometry type in postgresql again?

One benefit of this tool is to upgrade existing st_geometry type. But for a fresh install, is it necessary?


texts: Calculating the area of multiple polygons using GeoPandas

texts:
I am writing a Python script to calculate the area of polygons, as I have a small polygon, I am getting the area of polygon 0.000. A shapefile contains multiple shapefiles. I am thinking due to the area unit, it is printing 0.000. Shapefile projection is CRS EPSG:4326 - WGS 84 - Geographic.  I have to change the unit of area to Hectare in the script to calculate the area. I need to change the unit of area hectare(ha) in the below script.
Below is my working code, please let me know, how to change the area unit the below code
# Import necessary modules
import geopandas as gpd

# Set filepath (fix path relative to yours)
fp = r"D:\Data\Merged_Shapefile.shp"

# Read file using gpd.read_file()
data = gpd.read_file(fp)
selection = data[0:200]
for index, row in selection.iterrows():
    poly_area = row['geometry'].area
    print("Polygon area at index {0} is: {1:.3f}".format(index, poly_area))


texts: Programmatically accessing declarative ToggleButton?

texts:
I've just started with Dojo/Javascript and even though I kind of understand all the AMD way of coding, I keep getting stuck on these (probably) easy things.

How can I programmatically access a widget created in a declarative way?

I have a toggle button created like this:

<button id="toggleLeg" data-dojo-type="dijit/form/ToggleButton" data-dojo-props="iconClass:'dijitCheckBoxIcon', checked: true, onChange: function(value){legendVisibility(value);}">Legend/Layers</button>

And I would like to call the legendVisibility() function when clicking on it using the true/false returning values from the button.

The function simply toggle an accordion container on and off using its CSS property (the accordion is included into a content pane with id="legPanel").

    function legendVisibility(value) {
        checked = value;
        var rightPane = dom.byId('legPanel');
        if (checked == false) {
            //I want the right pane out
            domStyle.set(rightPane , 'display', 'none');
        } else {
            //I want the right pane in (preserving its content)
            domStyle.set(rightPane , 'display', 'block');
        }
    }

If I put it outside the main require function it clearly throws a ReferenceError saying that "dom is not defined" when declaring the var rightPane.

How can I correctly include it into the require function and programmatically use the widget?

Should I get rid of the onChange: function(value){legendVisibility(value);} when creating the widget?

The modules I need are "dojo/dom" to get the reference of the node and "dojo/dom-style" to change the CSS properties.

Do I need the "dojo/on" too to use the event propagated by the click?
(I also tried with the Toggler function with no luck!)

        registry.byId("toggleLeg").on("change", function() {
            checked = this.value;
            var rightPane = dom.byId('legPanel');
            if (checked == false) {
                //I want the right pane out
                domStyle.set(rightPane, 'display', 'none');

                registry.byId("mapDiv").resize();
                registry.byId("mapDiv").reposition();
            } else {
                //I want the right pane in (preserving its content)
                domStyle.set(rightPane, 'display', 'block');

                registry.byId("mapDiv").resize();
                registry.byId("mapDiv").reposition();
            }
        });


texts: How to attach one spatialite database to another in db manager

texts:
I'm trying to attach one database to another so that data from it can be used in queries.

I'm using

ATTACH DATABASE '{full url}' AS {database_name};

No errors actually occur, but when I run

PRAGMA database_list

Only the main database shows.


texts: Combine two grids with a small area of overlapping pixels to a new raster in QGIS/SAGA

texts:
The first grid is 1 byte and the second is 2 byte.

The new raster should have a blended transition.

How can i do this?


texts: Using OpenLayers 2, GeoServer WMS, and OSM CloudMade?

texts:
While I have some familiarity with GIS, this is my first attempt using OpenLayers.

I am trying to display the following layers in OpenLayers:

CloudMade Base Map
WMS Data overlay from my GeoServer  (http://54.244.249.54/geoserver/PalmBeach/wms?service=WMS&version=1.1.0&request=GetMap&layers=Florida2&styles=&bbox=-82.97059340790229,24.320616400434016,-80.1393763936513,26.907996410185536&width=512&height=467&srs=EPSG:4326&format=application/openlayers)

I have used the CloudMade sample as a starting point (http://developers.cloudmade.com/wiki/21/CloudMade_Tiles), and have successfully changed the parameters to center and zoom the map, but I cannot get the WMS to display as a layer, despite it being listed as a layer in the layer selector.

I'm not certain if its a projection issue (I changed the projection in GeoServer to match the Open Layers projection EPSG:4326) or just my unfamiliarity with OpenLayers, but would appreciate any guidance.

My code is below:

<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN"  "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">

    <title>Using CloudMade Tiles with OpenLayers</title>
    <link rel="stylesheet" href="openlayers/theme/default/style.css" type="text/css" />


texts: With OpenLayers is there anyway to keep existing features whilst adding more on each refresh?

texts:
Currently I use a Cluster strategy and a Refresh Strategy and get KML data from a url.

The KML data plots points on a vector layer and the cluster strategy groups close features together.

At each refresh interval the features are removed and new data is read from the server.

This is currently working.

What I want to do is the same as the above but when the refresh occurs I don't want the existing features to be removed. Instead I want the new features to be added to the existing features and the clusters to be recalculated.

I've spent a couple of days trying to find out how to do this but am failing miserably at the moment.

Any thoughts?


texts: How can I avoid a specific area of raster?

texts:
I have a raster polygon. I need to avoid certain area (another raster) from that raster. I need to avoid a watershed sub basin area that has its area on other basin also. How can I remove that basin?

For example:

I have the raster or shapefile of a country and I also have the shapefile or raster of a state in that country. I need to get the shapefile or raster of the country without the state I have.

I need to mask, but if I extract by mask, it gives the output as the boundary of the layer that we specified, but I need the rest of the area from the input.


texts: To project, or not to project? Extracting raster values with R

texts:
I have a raster from which I need to extract values at point locations using raster::extract(raster_stack, points_sf).

The raster is at 30 arc-second resolution and came from WorldClim with proj4string: "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"

i.e. it's coordinates are geographic on an ellipsoid. My point location data has the same lon-lat WGS84 CRS.

I've noticed that about 3/4 of the values I extract differ depending on whether or not I've transformed both the raster and the points to a planar projection with proj4string:

"+proj=laea +lat_0=43.15268 +lon_0=-70.30744 +x_0=12019341.4
+y_0=1730136  +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs".

My question is, should I expect points to move between raster cells this much as a result of projection? And which method is more accurate? In case spatial scale matters, the bounding box for the points is:

     xmin      ymin      xmax      ymax
-73.68631  41.01861 -67.89930  45.23919

Please let me know in the comments if this question needs a reproducible example.


texts: Right way to compute monthly MODIS Albedo Daily 500m (MCD43A3.006)

texts:
I am trying to compute monthly composite albedo from daily MCD43A3.006 product.
From the description, it is said that each 500m/pixel daily image is generated using 16 days of data, centered on the given day. I am not quite sure what is implied.
So my current program just takes all observations within a month and compute average value.
I wonder if it is the right way to do it or not.
Link to earth engine code
// compute monthly composite albedo
var monthlyComposite = function(modcol, id, folder, scale, maxPixels, region, nimg) {
    nimg = nimg || 500;
    scale = scale || 500;
    maxPixels = maxPixels || 1e11;

    var composite = modcol.mean();

    var palette =['000000', 'FFFFFF'];
    Map.addLayer(composite.select('Albedo_WSA_vis'), {min:0, max: 500, palette: palette}, 'Monthly albedo' + id);

    Export.image.toDrive({
        image: composite.toFloat(),
        description: id,
        folder: folder,
        fileNamePrefix: id+"_albedo_",
        region: region,
        crs: 'EPSG:4326',
        scale: scale,
        maxPixels: maxPixels})
};

var getImages = function(date1, date2, id, shape) {
    var imgs = ee.ImageCollection('MODIS/006/MCD43A3')
                          .filter(ee.Filter.date(date1, date2));
    imgs = imgs.select(['Albedo_WSA_vis']);
    monthlyComposite(imgs, id, "Global_albedo_count", 500, 1e12);
};

var year = 2012
var year_str = "2012_"
var next_year = 2013

// month 1
var Jan_1 = ee.Date.fromYMD(year,1,1);
var Jan_2 = ee.Date.fromYMD(year,2,1);
var Jan_id = year_str.concat("01");
getImages(Jan_1, Jan_2, Jan_id);


texts: Time slider for two things change in map

texts:
I want to use the time slider in GEE to change the year, while showing the crop distribution and groundwater for the year.
But my code can't control the display of both at the same time.
Is there a way to solve this?

// Run this function on a change of the dateSlider.
var showcrop = function(range) {
  // Asynchronously compute the name of the composite.  Display it.
  range.start().get('year').evaluate(function(name) {
    var data = cropLandcover.filter(ee.Filter.calendarRange(name,name,'year'));
    var layer = ui.Map.Layer(data, cropLandCoverVis , name + ' Cropdata');
    Map.layers().set(0, layer);
  });
};

var showwater = function(range) {
  // Asynchronously compute the name of the composite.  Display it.
  range.start().get('year').evaluate(function(name) {
    var data1 = table.filter(ee.Filter.calendarRange(name,name,'year'));
    var layer1 = ui.Map.Layer(data1, {color: 'FF0000'}, name + ' Waterdata');
    Map.layers().set(0, layer1);
  });
};

// Asynchronously compute the date range and show the slider.
var dateRange = ee.DateRange(start, end).evaluate(function(range) {
  var dateSlider = ui.DateSlider({
    start: range['dates'][0],
    end: range['dates'][1],
    value: null,
    period: 365,
    onChange: showwater
  });

https://code.earthengine.google.com/80915b70d76a37d646f865d40ab7dc2d


texts: Calculate field tool vs Field Calculator

texts:
I have successfully populated a field in Arc using the Field Calculator, however I would like to re-create this in python as part of a tool I am making. I was going to do this buy using 'Calculate Field' and copying the python snippet. However I can not get this to work.

This is the successful expression that works in the Field Calculator:

However, when I try the same in Calculate field I get this error:

This is how I am running it:


texts: Description for pixel values in different resampling methods

texts:
so I am looking for the description of pixel values for the following methods:

Nearest Neighbor
Bilinear Interpolation
Cubic Convolution
Bicubic Spline
LaGrange

If someone has got some knowledge of it, I would be very thankful


texts: How to georeference a LANDSAT scene in QGIS 1.8?

texts:
I tried to georeference a LANDSAT scene in QGIS 1.8 but I failed.
My first question is
Can I use geografical coordinates as x,y coordinates of control points? or I need projected coordinates (from Campo Inchauspe/argentina 4 for example)?

I used the geographical coordinates (lat; long in dd mm ss) and I obtained a black small (less than 1 meter) image.


texts: How to remove repetitive text in attributes?

texts:
I would like to clean up fields that contain repeated "Street Names". As seen below, I have fields "Name_From" and "Name_To". These fields sometimes repeat the same name (highlighted). I would like to remove these repeated street names.

As seen below in the Name_From field, the highlighted record says "Preller & Preller". I would like to end with only "Preller". ("Preller & Preller" becomes only "Preller".)

This repetition is also present in the Name_To field.

.


texts: How to prevent Google My Maps from truncating lat long values

texts:
My Lat and Long coordinates get truncated from 37.747574 and -25.675669 to
37.74757, -25.67566 when I import them into Google My Maps.

Is there a way to prevent this from happening?


texts:
I have one ASCII text file with eight columns. The first four provide the x-coordinates (e.g. longitude) of the four corners of a grid cell, while the last four provide the y-coordinates (e.g. latitude) of the four corners of the same cell in the same order. Furthermore, the fourth x-coordinate of line k is equal to the first x-coordinate of line k + 1, and the third x-coordinate of line k is equal to the second x-coordinate of line k + 1; similarly for the y-coordinates. This means that the cells described in lines k and k + 1 are adjacent having two corners in common. In this context, the number of lines in the file is the number of grid cells.

Is there a way to draw in QGIS the grid cells described in this way?


texts: python handling case in address (lowercase,uppercase)

texts:
I'm working on a python QGIS plugin that takes address info from specific parcels and lists it. The data is in various cases and I'm having a hard time figuring out how to handle it.

My data looks like this:

num addr                city        zip
45  NE   CRESCENT DR    Issaquah    98027

The script merges the attributes into a string:

45 NE CRESCENT DR, Issaquah, WA 98027

I've been using .title() which fixes the case in CRESCENT to Crescent, but messes up the NE and DR which should stay caps. Additionally for numerical streets like 234th Ave I need the "th" to be lowercase.

I would like the output to be:

45 NE Crescent DR, Issaquah, WA 98027

Any Advice?


texts: Copy and paste value from a feature sharing the same id

texts:
I would like to copy and paste the value from another feature sharing the same id in QGIS.
The feature to take the value from is always the only one with a value other than NULL per "id".
Here what I am trying to achieve:


texts:
I have a NetCDF file of ice velocities of the Greenland ice sheet but they are in polar stereographic coordinates (y: -655750 to -3355750 and x: -639750 to 855750 meters).

I am wondering if there exists a command in NCO or CDO that can convert from NSIDC polar stereographic north (which is in meters) to WGS84 degrees. I checked this ncremap module, but what is this 1x1.nc file in the NCO documentation? Is it just some random data in the WGS84 lat/lon coordinates? Does it have to be the same extent as the original data? Will this module automatically recognize the NSIDC Sea Ice Polar Stereographic North projection in the original file?

Below is some info about the NetCDF file:

Format:
           netcdf4
Global Attributes:
           Conventions                = 'CF-1.6'
           product_version            = 'ESP 2.1'
           date_created               = '2018-06-25'
           title                      = 'greenland_iv_500m_s1_20171228_20180228_v1_0'
           institution                = 'ENVEO IT GmbH'
           summary                    = 'Greenland IV from Sentinel-1'
           project                    = 'ESA Greenland Icesheet CCI'
           tracking_id                = '4d99b8ce-9418-4a84-aa47-92d0dd90a746'
           time_coverage_start        = '2017-12-28'
           time_coverage_end          = '2018-02-28'
           sensors_used               = 'Sentinel-1'
           methods_used               = 'Offset-tracking'
           velocity_projection_method = 'Horizontal Flow'
Dimensions:
           y = 5401
           x = 2992
Variables:
    crs
           Size:       1x1
           Dimensions:
           Datatype:   int32
           Attributes:
                       grid_mapping_name                     = 'polar_stereographic'
                       standard_parallel                     = 70
                       straight_vertical_longitude_from_pole = -45
                       false_easting                         = 0
                       false_northing                        = 0
                       unit                                  = 'meter'
                       latitude_of_projection_origin         = 90
                       spatial_ref                           = 'PROJCS["WGS 84 / NSIDC Sea Ice Polar Stereographic North",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",70],PARAMETER["central_meridian",-45],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY["EPSG","3413"]]'
    y
           Size:       5401x1
           Dimensions: y
           Datatype:   double
           Attributes:
                       units         = 'm'
                       axis          = 'Y'
                       long_name     = 'y coordinate of projection'
                       standard_name = 'projection_y_coordinate'
    x
           Size:       2992x1
           Dimensions: x
           Datatype:   double
           Attributes:
                       units         = 'm'
                       axis          = 'X'
                       long_name     = 'x coordinate of projection'
                       standard_name = 'projection_x_coordinate'
    land_ice_surface_velocity_magnitude
           Size:       2992x5401
           Dimensions: x,y
           Datatype:   single
           Attributes:
                       _FillValue   = 3.402823466385289e+38
                       units        = 'meter day-1'
                       description  = 'magnitude of horizontal ice velocity'
                       grid_mapping = 'crs'
                       coordinates  = 'y x'


texts:
I have loaded my bike tours as .gpx files into QGIS.

It already looks quite nice but I was wondering if it might be possible to create something similar to the strava heatmap with tools of QGIS.

I have tried to generate points along the lines and use kernel density but it just does not look nice as a raster layer.

Any suggestions?


texts:
I have a matrix like this:

   10 0.646519601   0   0.085929349 0.255477717 0.012073334
   40 0.138403594   0   0.065979943 0.776693742 0.018922721

I want to keep the first column, and get rest numbers sorted by row. The result should be like this:

   10 5 1 2 4 3  #0.64 is the largest; 0 is the smallest in the row#
   40 4 1 3 5 2


texts:
i've got some problems with QGIS and I do need your help..

I've got 2 point layers. Layer 1: 1700 points, Layer 2: 1704 Points.. Both are ESRI-Shapefiles

1700 Points are identical. only 4 points are different (different coordinates).
Now I want to compare both Layers and create a new layer from the result, so the final layer should only have these 4 points. My problem is the following:

If i choose Layer 1 first, the result is 1700 points. If I choose Layer 2 first, the result is 1704 points.. I guess the problem is, that both layers are completely different (not a single attribute is identical).

how I can I get the correct result, based on the geo-coordinates? That should give me only 4 points output. I have to do this for around 46.000 points, so a hint for my problem would be awesome!

I hope you understand what I'm talking about, bc I'm new to QGIS and my english is not the best..

king regards,

Patrick

EDIT:
Field names Layer 1:

OBJECTID, FID_Hausko, nasoid, lebenszeit, lebensze_1, bundesland, ortschlues, ortschlu_1, ort, ortsteil, strassensc, strassen_1, hausnumm_1, unverschlu, infotext, xwert, ywert, aktuell, FID_Gemein, OBJECTID_1, gemeindezu, nasoid_1, id_1, lebensze_2, lebensze_3, anlass1, anlass2, schluessel, bezeichnun, gemeindeke, herkunft, herkunft_e, aktuell_1, hatgeometr

Field names Layer 2:

id, lat, lon, ags, kid, plz, ort, ortsteil, str, hnr, adz, kvz_id, apl_bez, nb, adsl, vdsl, vec, versorgt_d, knb, hh_me, hh_ibv, vzk_ist, db_ist, kabelnetz, mobilfunk_, mobilfun_1, mobilfun_2, mobilfun_3, fttb, ist_versor, hh_d_g, hh_u_g, ausbaug, exist, kein_Ausba


texts: How to Increase the length of column type geometry in postgres(pgAdminIII)

texts:
I am just wondering is there any way to increase the default length of column type geometry in postgres.

Since we are trying to save large number of vertices polygon in to the geometry column geom in postgres database.It is saving nothing.


texts:
I am having issues with my SQL query in the arcpy.da.UpdateCursor(). The code seems to work fine if I remove the SQL query  so I am fairly certain that it is the issue. The syntax that I would use in ArcMap is last_edited_date >= GETDATE() AND last_edited_user NOT LIKE 'ESRI'. How do I write this in Python?

import arcpy

arcpy.env.overwriteOutput = True

#Define Variables
fc = r'C:\Users\...\FSAArchiveTest.gdb\FSA_Archive_Test'

try:

    edited_archive = "last_edited_archive" #This is a STRING field
    edited_user = "last_edited_user"  #This is a STRING field
    edited_date = "last_edited_date" #This is a DATE field

    #THIS IS THE SQL QUERY THAT IS NOT WORKING
    sql = "last_edited_date >= GETDATE() AND last_edited_user NOT LIKE {0}".format('ESRI')

    with arcpy.da.UpdateCursor(fc,[edited_archive, edited_user, edited_date], sql) as cursor:
        for row in cursor:

            row[0] = (row[1] + " " + str(row[2]))
            cursor.updateRow(row)

        print("It worked")

except:

    print("did not work")


texts: Editor Tracking setting on a FGDB preventing data display in QGIS

texts:
I have a File Geodatabase containing a polyline feature class. I can drag this into QGIS (3.0.0) and view the data in a read-only format. This is fine, for the purposes of the majority of users of this data they simply need to be able to reference the data.

This data, however, will be subject to some editing by multiple users in an Arc environment. As such, I want to have Editor Tracking enabled to keep a record of changes that are occurring. However, once Editor Tracking has been enabled, this data no longer displays in QGIS. I can still drag the layer in, but the attribute table is empty, and nothing shows in the project view screen.

What is the exact nature of the change that occurs to a feature class with Editor Tracking enabled that might be the cause of this behaviour?

What can I do with this feature class so that I am still able to keep track of edits whilst simultaneously providing read-only access to users in a QGIS environment?


texts: Exporting to pdf in QGIS 2.8 composer does not respect line styles

texts: I'm using QGIS 2.8. When I export a map using "export to pdf" in Composer, the style of lines of my vector data in the resulting pdf file is not what is seen on the screen. When I export to jpeg it does not do this and lines are correctly displayed. Here are two screenshots: the first is how my map should look like (exporting to jpg); the second is how the maps comes out when I export it to pdf. The style I defined for the contour lines is not respected when I export to pdf (I get thicker lines and all the same).


texts: Extracting Climate classifications from raster using multiple polygon shapefiles

texts:
I have performed a cluster analysis on geographic locations (national and sub-national level) in Matlab.  I have also written code to write the shapefiles for each of the locations in a cluster to it's own folder.

Using the Köppen classification found here, I have added this raster to ArcGIS 10.2.2, such that each pixel is a value between zero and thirty. I then add all of my shapefiles from the specified folder.

I want to produce a table that gives me the name of the country/region (which is in the polygon shapefiles) along with the Köppen classifications that were found within that region.

I started by merging all the shapefiles into one and then tried a number of different routes:

I did a raster clip using my new shpefile to clip the raster and then created an attribute table for the raster.
I have tried converting the raster to a polygon after clipping it and creating a spatial join.
I also tried an 'extract by mask' using the 'feature outline masks' tool to convert the raster clip to a mask.

Each time I try something I end up with either:

(1) An attribute table that has no region names, but does have the Koppen number/code along with a count of the number of times the Koppen code occurs. This is, I assume, a count of the number of gridcells with that classication.

or

(2) A table with the shape lengh and area of polygons, along with the classification code.
I havent been able to produce an attribute table showing the name of the region, along with the Koppen classification number(s) for that region (one region can have more than one classification). I don't want counts, just to know what regions are found within each country/region.

Any suggestions?


texts: Strange coordinates for San Francisco: 5979762, 2085921 from Shapefile

texts:
I have a Shapefile (.shp, .dbf, .shx, .prj) that's supposed to represent street data for San Francisco, which I'm writing a parser for. The problem is, when I parse the bounding box according to the ESRI format for Shapefile here, I get back around (5979762, 2085921) for the min coordinate of the bounding box (note that I actually get floats, but I rounded them down here for clarity). These are not UTM coordinates for SF, as I originally thought, and they're definitely not lat/lons.

I'm not sure what units these are or how to interpret them. The content of the .prj file is as follows:

PROJCS[
    "NAD_1983_StatePlane_California_III_FIPS_0403_Feet",
    GEOGCS[
        "GCS_North_American_1983",
        DATUM[
             "D_North_American_1983",
             SPHEROID["GRS_1980",6378137.0,298.257 222101]
        ],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.0174532925199433]
    ],
    PROJECTION["Lambert_Conformal_Conic"],
    PARAMETER["False_Easting",6561666.666666666],
    PARAMETER["False_Northing",1640416.666666667],
    PARAMETER["Central_Meridian",-120.5],
    PARAMETER["Standard_Parallel_1",37.06666666666667],
    PARAMETER["Standard_Parallel_2",38.43333333333333],
    PARAMETER["Latitude_Of_Origin",36.5],
    UNIT["Foot_US",0.3048006096012192]
]

My problem is that I really don't have any experience with .prj files, so I don't know how to interpret this format. Can anyone explain what this format is and how I should interpret it, or where I can find any more resources on how to do so?


texts: PostGIS st_difference similar to ArcGIS Erase?

texts:
I am attempting to replicate the ArcGIS erase tool in PostGIS.

Here is the problem I am working with

Two Layers

land use land cover (150,000) multipolygon table
multipolygon soils table (30,000)

both have spatial indexes

here is my query

CREATE TABLE agriculture.important_ag_soils_undeveloped AS
--cte selecting farmland data
with farm_to_clip AS (
  SELECT *
  FROM base_layers.ssurgo_soils_highlands_updated_stateplane
  WHERE farmlndcl IN ('Farmland of statewide importance','Farmland of unique importance','All areas are prime farmland')
),
--cte selecting land use land cover data to erase from farm_to_clip data above
lulc_to_erase AS (
  SELECT shape FROM base_layers.lulc_huc_hl_clip
  WHERE lu12 IN (1110,1120,1130,1140,1150,1200,1300,1500,1600,7300,1211,1400,1410,1411,1419,1420,1440,1810,1462,1463,1499,1710,1800,1804))
--now I am running this sort of messy query where I perform a st_difference for each farm polygon record and run a subquery as the second input unioning the landuse land cover where it intersects the farm polygon
-- which yeilds the difference, and I wrap the st_difference function in the coalesce function so if it is st_difference(f.shape,NULL) then f.shape is returned for that row.
SELECT
  f.*,
  coalesce(
    st_difference(f.shape,
      (SELECT st_union(l.shape)
       FROM lulc_to_erase AS l
       WHERE st_intersects(f.shape,l.shape)
      )
    ),
    f.shape) AS geom
FROM farm_to_clip f;

this query takes about 5 minutes to run and I know there must a more efficient way of writing this query.

from @tilt suggestion I removed a the subselect and coalesce function and just used a left join with the st_difference function here is my updated query:

create table agriculture.important_ag_soils_undeveloped2 as

WITH farm_to_clip AS (
    SELECT *
    FROM base_layers.ssurgo_soils_highlands_updated_stateplane
    WHERE farmlndcl IN ('Farmland of statewide importance','Farmland of unique importance','All areas are prime farmland')
    ),

lulc_to_erase AS (
    SELECT shape from base_layers.lulc_huc_hl_clip
    WHERE lu12 in (1110,1120,1130,1140,1150,1200,1300,1500,1600,7300,1211,1400,1410,1411,1419,1420,1440,1810,1462,1463,1499,1710,1800,1804))

SELECT f.*,
    ST_Difference(f.shape,l.shape) geom
FROM farm_to_clip f
LEFT JOIN lulc_to_erase l ON ST_Intersects(f.shape,l.shape);

Problem: the new method of using the left join returns a significant more amount of polygons in the new layer than the first query I used.

My question becomes which query gives me the farm layer erased by the lulc like the way ArcGIS performs its erase?

UPDATE

I used the arcgis erase tool on the farm_to_clip and lulc_to_erase and the layer that was created is the same layer that I created in my first query. My second query did not seem to work correctly at all, but I cannot figure out why..?

I am going to try unioning my lulc_to_erase geometry before I send it through the select st_difference function

farm = light transparent green
lulc = hatched orange

QUERY 1 RESULTS

QUERY 2 RESULTS


texts: Converting SLD to CSS style in GeoServer?

texts:
I am trying to convert SLD to CSS styles in latest geoserver version 2.8.1.

I have added css plugins(from geoserver-2.8.1-css-plugin.zip) but it is not working properly.

Is there any other plugin for converting from SLD to CSS style?

What is the solutions for this?


texts: Many-to-many relationship and QGIS Forms

texts:
I have a PostgreSQL/PostGIS database with a base table (Pois_points) and two subtables which constitute a many-to-many relationship with the base table, as you can see on the picture below (click on Many to many relationships).

I want to build a QGIS form that depicts this relationship between the tables when adding and editing items. I have seen examples of a one-to-many relationship in QGIS, which is relatively easy. But many-to-many is not documented or maybe not possible at all?


texts: Install Georeferencer Plug-in QGis 1.7.0 OS X

texts:
I'm trying to install the Georeferencer plug-in in QGis 1.7.0 on Mac OS X (10.5.8).
However the plug-in doesn't show up in the list. I have used "enable 3rd parties directory" and I have checked "manage plug-ins", it's not installed yet.

I believe that I can't install a newer version of QGis on OS X 10.5.8, if I'm wrong please correct me...
Any suggestions? I'd be very thankful. If there's no way to install the plug-in - are there other ways to georeference raster imagery?


texts: Using Multiple Versions of ArcGIS API for JS in Same Application

texts:
We've a large web based GIS application running which uses an older version of ArcGIS API (v3.7) for JS. The ArcGIS JSAPI v3.7 is lacking some functionality which are available in v3.15. Now I also need to include ArcGIS JS API v3.15 for a specific function.

Changing API version for whole application is not an option.

So, Is it possible to include some of the modules from v3.15 for a specific function and destroy when that function is executed?

Scenario:

//code uses v3.7
// ....
// ...
function executeLatestFeature() {
  // need to use some functions from v3.15
  //...
}
//code uses v3.7


texts:
I found this answer on my question.
But when I try to do it, I get this error:

The command:

pdal translate D:\Information\bigrotate_frontview.pcd D:\Information\output.las

I get this error:

Failed to find match for field 'intensity'.
PDAL: writers.las: Unable to convert scaled value (NaN) to int32 for
 dimension 'X' when writing LAS/LAZ file D:\output.las.

So I searched for that error and I found this.

So I tried the option --writers.pcd.xyz=true

My command looks like this:

pdal translate --writers.pcd.xyz=true D:\file.pcd D:\output.las

I got the following error:

PDAL: Argument references invalid/unused stage: 'writers.pcd'.


texts: Aligning polygons in ArcGIS for Desktop?

texts:
I am using ArcGIS 10.3 for Desktop.

I am working on a map that I have divided into 5 zones. I have drawn 5 polygons to represent each zone. I need to align the polygons where the edges meet.  I am trying to use the "Align Edge" tool in the Topology toolbar. I select the tool, click an area of one polygon where there is a sliver I need to remove.  I see a dotted line and then click the other polygon where I want the edges to align.

I get an error message that reads:

  The edge cannot be modified. The operation was attempted on an empty
  geometry.

I don't know why this is happening or how to rectify.


texts: Locating plugins in the QGIS gui programatically from PyQGIS

texts:
I've let my love of shiny new QGIS plugins get the better of me, and I'm finding it difficult to remember which menu / submenu / panel that a given plugin is in. I'm using well over 100 now, and going on a purge usually means I find another use for the plugin I just removed...
So I thought I'd try to write a plugin to get that info so I could search for it. Things that the Plugin manager knows, plus (most importantly) how to locate it in the GUI.
Processing is great for finding a lot of the built-in QGIS functionality, but it doesn't find plugins which weren't written to make use of its framework.
Looking at the answer to listing installed plugins via pyqgis I was able to get a list of plugin names.
There seems to be a class, QgsPluginRegistry, which looks promising in the C++ code, but I can't find a SIP definition file for that, so I don't think this class is exposed to Python.
Using a bit of tracing through the QGIS source, I was able to write a short snippet to sniff the configuration files for the plugins
for x in qgis.utils.findPlugins("/home/steven/.qgis2/python/plugins"):
    # gets a (name,ConfigParser) tuple
    print "_"*80
    print x[0] # plugin name
    for key in x[1].options('general'):
        print "\t%s = %s" % (key, x[1].get('general',key))

This gives me access to the same info as the QGIS Plugin Manager... changelogs, authorship, category and so on.
qgsAffine
    name = Affine Transformations
    qgisminimumversion = 2.0
    description = Apply affine transformations to selected geometries.
    version = 1.0.2
    author = Mauricio de Paulo and Erik Timmers
    email = [redacted]
    changelog = [snipped]
    experimental = False
    deprecated = False
    tags =
    homepage = [redacted]
    repository = [redacted]
    tracker = [redacted]
    icon = icon.svg
    category = Vector

What it doesn't tell me is where, exactly, is it in the menu structure? Can this be got from PyQGIS in a clean way?
Or to find this, do I need to start sniffing round the Qt GUI classes or writing python parsers (doable, but seems a bit over the top?) The category value (vector, raster etc) helps narrow the search, but there are plenty of plugins (e.g. MMQGIS) which create their own menus.


texts:
I want to display "map tips" in QGIS with python for my custom map canvas QgsMapCanvas().

How can I do this? I can switch manually with View/Map Tips but it won't work with my custom canvas in my plugin, it works only with iface.mapCanvas().


texts: Getting random coordinates based on country

texts: Are there any webserveices/APIs that can take a country name and return a random coordinates (not including the sea/ocean sections) inside that country?


texts:
I have an OpenLayers map which uses an ImageLayer as the base layer and I want to create a custom projection so that (0, 0) is at the top left of the image rather than the bottom left.
I tried setting the projection's axisOrientation to 'esu' but that did not change the axis in any way. I tried creating my own proj4 projection, but the origin is still appearing at the bottom left of the image.
  proj4.defs(
            'newproj',
            '+proj=somerc +axis=esu +units=m +no_defs'
  );
  register(proj4);
  this.proj = get('newproj');
  this.proj.setExtent(this.getImageExtent());

  var view_options = {
            projection: this.proj,
            center: getCenter(this.getImageExtent()),
            zoom: 2,
            extent: this.getImageExtent(),
  };
  this.view = new View(view_options);

  this.map = new Map({
            target: 'map',
            view: this.view,
  });