texts: Convert GDAL command line commands to Python script or code

texts:
I'm looking to convert my GDAL commands to Python code using either GDAL or Rasterio, as I've heard this is better to use.

I've been scratching my head at this and I could be missing some simple stuff, but I've tried to check as much Python API documentation as possible, but I've been unable to find any for the potential complexity of my commands. They are;

gdal_translate -co COMPRESS=JPEG -co JPEG_QUALITY=90 -co PHOTOMETRIC=YCBCR -co TILED=YES -co BLOCKXSIZE=128 -co BLOCKYSIZE=128 uncompressed.tif compressed.tif

gdaladdo --config COMPRESS_OVERVIEW JPEG --config PHOTOMETRIC_OVERVIEW YCBCR --config INTERLEAVE_OVERVIEW PIXEL -r gauss compressed.tif 2 4 8 16 32 64 128 256 512 1024

With the aim being I'd like to be able to run scripts on many folders of TIFFs which will not be georeferenced.

It doesn't have to be using GDAL or Rasterio these are just the only Python modules I know that should be able to do the job.


texts: Mosaic to New Raster doesn't actually output anything

texts:
I am attempting to Mosaic some Raster files together using ArcMaps/Catalog's Mosaic to New Raster tool. It goes through the motions of adding each .tif file we are using, then after a bit it says that is has Succeeded.

If I check the reports and look at where it says it outputs there it has the GDB I specified.

However, when I navigate to the GDB there is nothing there. I have run this a number of times and each time is the same. Nothing appears in any of the GDBs I attempt to send the new Raster Mosaic to.

Anyone run into this?


texts: Writing QGIS field calculator expression?

texts:
I have an attribute table with a real data type column named "Volume". Some of these cells have no value, and I want to replace the NULL with a dummy value like -999.
In the Field Calculaor I select "Update Existing Field", select the "Volume" field and I attempt the following:

1) I code CASE WHEN NULL THEN -999 END - Result: the entire column becomes NULL.

2) Then I try CASE WHEN  "Volume" IS NULL THEN -999 END - Result, the NULL valued cells now show -999 BUT all the other formerly filled cells are now NULL.
After many hours of trying other expression variations I can't figure this out. My guess is I am missing a basic concept.

Oh, and how can I "undo" a bad field calculator update besides quitting QGIS and  discarding the updates.


texts:
In the following how would I customize the popup options?  I want to add in the option for the keepInView  to be true.

var myLayer = L.geoJson(polygon, {
  onEachFeature: yourOnEachFeatureFunction
}).addTo(map);

function yourOnEachFeatureFunction(feature, layer){
  if (feature.properties.name) {
      layer.bindPopup(feature.properties.name);
  }
}


texts: Modelbuilder - iterate consecutive workspaces and select specific rasters from each

texts:
This is probably easiest explained with a picture,

I'm trying to iteratively perform a calculation using specific rasters from each two consecutive folders in a folder. For example from above picture, Iterate through the Date Folders in the Year Folder, and select the B4.TIF from date 1 and date 2, then process using these two tifs.

I can parse and slice filenames for the output files and folders, however just cannot figure out how to selectively iterate?

Have tried to parse filename and slice/calculate value for an inline variable based on each date to use as a wildcard but I'm just getting rasters from the last folder in the list.

I know using python would be easier, but I don't have the skills hence the post!


texts:
I work in Qgis 2.18.17. I try to install "Processing R Provide" as it shows in video
https://www.youtube.com/watch?v=nEJJE4DRiPQ&t=193s
but as I understood it is usefull only for windows when I use Ubuntu. I founded this plugin in repository
https://plugins.qgis.org/plugins/processing_r/
But I could not find this plugin inside QGIs.
Also I tried to install Qgis library in R but seemed useles
How should integrate R into Qgis?


texts: How to build and deploy GeoExplorer on jboss 5.1.x?

texts:
I am trying to deploy opengeo GeoExplorer on jboss 5.1.x by using this command:

ant deploy -Dcargo.host=localhost -Dcargo.port=8080 -Dcargo.container=jboss51x -Dcargo.username=abc -Dcargo.password=abc123

but GeoExplorer fails to build successfully.

Here is the log containing some of the exceptions:
http://paste.debian.net/189517

I have deployed the geoexplorer.war file in jboss 5.1.x but deployement failed. I get the following error when I run geoexplorer:

  ringo/global not found or not readable

What do I need to do in order to deploy GeoExplorer in jboss 5.1.x?


texts:
I'd like to enable jsonp on GeoServer 2.3.0. The docs say to set ENABLE_JSONP to true.

A response from the GeoServer users list says to edit web.xml, (set System variable ENABLE_JSONP=true and outFormat=text/javascript )but I have no idea which xml tags to use and where to put it.

Here are installed versions:
GeoServer 2.3.0 (Web Archive)has been installed on a remote host with tomcat 7.0.39


texts: Performing multiple if statements with ArcPy UpdateCursor?

texts:
Beginner with Python trying to perform multiple calculations on a queried set of rows.  I'm am then looping thru a set of queries.  But I need to perform all calcuations on returned rows before starting next query in the array.  The first set of if else statements runs correctly.  The subsequent statements don't update every record.  Here's what I have so far...

for fl in flArray:
calcFields = ['FL', fl, 'EValue', 'Diff', 'OOC', 'E_GRTR', 'gCAT', 'eCAT']
queryFL = flArray[fl]
cursor = arcpy.da.UpdateCursor(out_layer_file, calcFields, queryFL)

 # Calculate difference
 for row in cursor:
  row[3] = row[2] - row[1]
  if row[3] > 0:
        row[5] = 1
  else:
        row[5] = 0
  cursor.updateRow(row)

 ### Calculate g Category

 for row in cursor:
  if row[1] < .15:
     row[6] = 0
  elif .15 <= row[1] <= .22:
     row[6] = 1
  elif .23 <= row[1] <= .3:
     row[6] = 2
  elif .31 <= row[1] <= .37:
     row[6] = 3
  elif .38 <= row[1] <= .45:
     row[6] = 4
  elif row[1] > .45:
     row[6] = 5
  cursor.updateRow(row)

 ### Calculate e Category

 for row in cursor:
  if row[2] < .15:
     row[7] = 0
  elif .15 <= row[2] <= .22:
     row[7] = 1
  elif .23 <= row[2] <= .3:
     row[7] = 2
  elif .31 <= row[2] <= .37:
     row[7] = 3
  elif .38 <= row[2] <= .45:
     row[7] = 4
  elif row[2] > .45:
     row[7] = 5
  cursor.updateRow(row)

 ### Calculate OOC

 for row in cursor:
  row[4] = abs(row[6] - row[7])
  cursor.updateRow(row)
  cursor.reset


texts: Is the QGIS cache directory automatically created?

texts:
This question is directly a result of a question posted earlier about the units of cache size:

QGIS - what is cache size measured in

I went to the Options menu in QGIS, to the Network tab.  Here is a screenshot of said tab.

As you can see, my cache directory path is set to C:/Users/Russell/.qgis/cache/.

I went to look at the directory to see how large the directory was or the total size of the files contained within.  To my surprise, the cache directory does not exist.  I have my folder options set to show any hidden files and folders, so it is not simply hidden from my view.

I've been using QGIS for a number of years now with no problems, and it does not seem slow.  This leads me to believe it is being cached in another location.

The computer is running Windows 7.  I installed QGIS using the OSGeo4W installer.  Is there a setting somewhere that overrides the QGIS setting?

It seems worth putting this out there to see if there are other cases of this happening, or if it is isolated to my machine.  As I said, I am not seeing an impact on performance, but that could simply be because the computer is speedy enough.  It would be nice to know why the files for my installation don't seem to match what is shown in the program itself.

Any insight would be appreciated.

Thanks!


texts: Problem with Y offset shift of LAS file (las2las)

texts:
I'm trying to set the origin of a georeferenced .las point cloud to be coordinates 0,0,0, essentially stripping the georeferencing. The only way I could think to do this was subtract the minimum X, Y and Z value of each point using lastools (las2las). I was able to do this with the X and Z coordinates, but when I subtracted the Y coordinates something strange happened: starting with a minimum Y value 2,887,191 I subtracted 2,887,191 and got 4,294,967. I haven't been able to get it to 0.

Anyone have any idea what's happening here?

EDIT:

To be clear: I'm interested in setting the minimum x,y and z values of the .las file to 0, 0 and 0, ignoring the header and stripping georeferencing entirely but preserving the spatial relationship of the points. The best mathematical way I can think to do this is to subtract the minimum values of each axis from every point. Does anyone know of tools that can do that?


texts: Preparing two images with different dimensions to use in Agisoft PhotoScan

texts:
I have two aerial photos that are flowed over the national forest. The photos were scanned by a contractor and they came in various dimensions for each photos.

The one I have had smaller dimension and the other had bigger dimension.

How do I make the smaller image bigger dimension so I can use it to perform overlay each other with Agisoft PhotoScan?

Both image I have do not have no coordinate and projection.

I have ERDAS Imagine 2018 and ArcGIS Desktop 10.5.1, and ArcGIS Pro 2.4.2 here on my computer.


texts: Latitude and longitude values in Google Takeout location history data sometimes out of range

texts:
I've exported my geolocation history data as a JSON from Google Takeout and the latitude/longitude values are out of normal range sometimes. It seems to happen for some locations only, while others are normal and correct.

Here's a small extract of three dates with the invalid lat/long values and where I actually was at the time, is there some common transformation I need to use here that fixes this:

Date        lat         long        Location
13/01/2014  423.3738877 106.6510714 Jakarta, Indonesia (-6.23,106.76)
28/10/2015  397.6754516 115.7765074 Perth, Western Australia (-31.91,115.88)
8/08/2017   409.237868  372.8135318 Miranda, Brazil (-20.13, -56.88)

In the raw JSON file, the above three data points look like this (in the same order as listed above):

{
  "locations" : [ {
    "timestampMs" : "1389633689564",
    "latitudeE7" : 4233738877,
    "longitudeE7" : 1066510714,
    "accuracy" : 20
  }, {
    "timestampMs" : "1446040136103",
    "latitudeE7" : 3974349226,
    "longitudeE7" : 1157765074,
    "accuracy" : 1212
  }, {
    "timestampMs" : "1502203327757",
    "latitudeE7" : 4094251339,
    "longitudeE7" : 3728135318,
    "accuracy" : 16
  } ]
}


texts: How to create a table from a point in Google Earth and display in ArcGIS (shapefile format)

texts: Is it possible to create attribute fields for points in Google Earth (KML format) and display them as an attribute table in shapefile format in ArcGIS?


texts: Using gdal2tiles XML output with Proj4Leaflet

texts:
I use gdal2tiles to create tiles from a referenced image (with profile parameter set to "raster"). Then I use Proj4Leaflet to display the results. In order to do that I need to specify parameters of CRS such as resoultions, bounds and origin. After processing i get tilemapresource.xml that seems to contain all the information needed:

<?xml version="1.0" encoding="utf-8"?>
<TileMap version="1.0.0" tilemapservice="http://tms.osgeo.org/1.0.0">
  <Title>OUTPUT.tif</Title>
  <Abstract></Abstract>
  <SRS>+proj=aeqd +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs</SRS>
  <BoundingBox minx="4081867.72157528297976" miny="-2321600.04517470393330" maxx="6014219.16503371112049" maxy="477065.79241934354650"/>
  <Origin x="4081867.72157528297976" y="-2321600.04517470393330"/>
  <TileFormat width="256" height="256" mime-type="image/png" extension="png"/>
  <TileSets profile="raster">
    <TileSet href="0" units-per-pixel="14828.59620879369504" order="0"/>
    <TileSet href="1" units-per-pixel="7414.29810439684752" order="1"/>
    <TileSet href="2" units-per-pixel="3707.14905219842376" order="2"/>
    <TileSet href="3" units-per-pixel="1853.57452609921188" order="3"/>
    <TileSet href="4" units-per-pixel="926.78726304960594" order="4"/>
    <TileSet href="5" units-per-pixel="463.39363152480297" order="5"/>
    <TileSet href="6" units-per-pixel="231.69681576240149" order="6"/>
  </TileSets>
</TileMap>

Assuming

origin = XML.Origin
resolutions = XML.untis-per-pixel array
bounds = XML.BoundingBox

I get correct Leaflet map - image sits where it should, coordinates returned by events are correct as well.

But bounds parameter actually stands for "Bounds of the CRS, in projected coordinates" according to the documentation (https://kartena.github.io/Proj4Leaflet/api/#l-proj-crs). And BoundingBox attribute in the .xml shows bounds of the tile layer, but not the whole CRS.

Docs and many tutorials also set resolutions parameter as array [... 2^(n+3), 2^(n+2), 2^(n+1), 2^n].

resolutions: [8192, 4096, 2048] // from the docs

However resolutions in the .xml file are never 2^n. Again, using them is fine as the result is correct.

I would be really grateful for some explanation of such behaviour:

Why resolutions of zoom levels are usually set to 2^n while in fact they are equal to some random numbers?
Should I really use extent of the image as bounds (extent of the CRS)?

Thanks in advance for your answer!


texts: How can I convert files from a county government to shapefiles?

texts:
I work for a company that creates software for multiple listing services and real estate agents. We're expanding our coverage area into a new market, and I've been given some files for two different counties which they are calling their "shapefiles". One file is an access database, and the other is a fixed width column text file. Is there some way for me to convert these files to actual ESRI shapefiles?

Here are some examples -

Here's what I'm given in the fixed width file for the "vector" column (from Palm Beach county, Florida):

A0CL10D10L12D40R24U40L2U10
A4L30D2L2D6R2D14L2D6R2D7D11L7D14L14D20R26R8U6L8D6U6R8D11R17U14R5U6L5D6U6R5D6R8D6R14U38L14D12L8CL5D7R5U7
A1L24D12L5D10R5D28R12U4CR22U20L10D10L12D10

I have no idea what the above could be. If you want to see the file here's a small sample: sample.txt

The Access file (from Broward County, Florida):

The access file just says "Long binary data" for the geometry columns. It has a lot of tables named GDB_xxxx (GDB_SpatialRefs, GDB_GeomColumns, etc). I'm guessing this has something to do with Microsoft's spatial format.

Any help you guys can give would be greatly appreciated!


texts: Is there an easy way to use Postgis GeoJSON in Openlayers 3

texts:
I was using google-vector-layers and leaflet-vector-layers from Jason Sanford (https://github.com/JasonSanford) to easily display, style and add customized popups for data from Postgis-databases.  This works in combination with a modified version of the PHP-Database-GeoJSON from Bryan McBride.

Is there nothing comparable to use with Openlayers 3 out there? I must admit, that I don't have the programming skills to write such a library. Might be that someone knows about a comparable code. After doing a lot of Googling around, I found a lot of answers to particular problems and managed to implement things on how to load GeoJSON using Ajax and boundingBox strategy, styling the vector layers and to add popups, but I'm still missing an easy way on how to put all these things together.

I'm wondering if there is maybe an existing solution I didn't find till today, since I think PostGis > GeoJSON > Openlayers 3 (including styling and displaying features via popup) should be a quite standardized way?


texts: Accessing specific layers in ArcGIS Map Service to toggle them on/off?

texts:
I have access to a number of ESRI Mapservices and would like to be able to access specific layers within the mapservices in my online map application. I am working on creating this with the ESRI Javascript API.

I have looked, at and hoped that I could use, this sample map:

http://help.arcgis.com/en/webapi/javascript/arcgis/demos/map/map_legendvisible.html

and substitute some of the data I have access to.

When I do substitute mapservices that we have here, the sample map will toggle on and off the entire mapservice and not the specific layer I was hoping it would.

Here is the mapservice I am trying to access with this code:

http://gis1.co.frederick.va.us/ArcGIS/rest/services/Planning/testplanning/MapServer

I would love to hear from someone about a solution. Nothing seems as straight forward as it looks...


texts: Assigning function to a key in QGIS Plugin

texts:
I have trouble assigning function to a hotkey in my plugin. I have a class that does things with a mouse click. I'd like to assign another function to a keyboard key or sequence. I've tried two different approaches:

https://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/snippets.html#how-to-call-a-method-by-a-key-shortcut
This assigns given key sequence as a shortcut, but does nothing.
def scale():
    iface.mapCanvas.zoomScale( 2500 )

shortcut = QShortcut(QKeySequence(Qt.ControlModifier + Qt.Key_Q), iface.mainWindow())
shortcut.setContext(Qt.ApplicationShortcut)
shortcut.activated.connect(scale)

Both don't return any errors, but just doesn't work. I made sure that key sequences aren't in use. What am I missing?


texts: Creating one dataset from several overlapping point datasets using ModelBuilder?

texts: I am trying to build a model in which I create 9 point data sets that are subsets of one large point data sets. The 9 small point data sets are suitability analyses for tree species. So one point can be suitable for several different species. I'd like merge, join, spatial join these 9 data sets into one data set.  Ideally, the resultant point data set would then have fields (species 1, species 2 species 3, etc.) that represent the tree species suitable for each point location.  I feel like this should be easy to do, but just can't figure out how to do it.


texts: Using Google Maps or Bing Maps provided as WMS in MapInfo Professional?

texts: Is there any way that I could use google maps or bing map or other free map provider as WMS in MapInfo Professional 10.0?


texts: Writing a python processing script with QGIS 3.0

texts:
Following the update to QGIS 3.0, it has become very difficult to find any information concerning the writing of processing scripts in QGIS 3.0.

@Underdark (see here ) has provided a basis for the skeleton. This code also seems to have been added in QGIS, when writing a new script from template (QGIS 3.0.2).

However, I couldn't find any way to help Python newbies like me to understand how to change that code, especially for the input and output layers.

My goal is to write a script taking 2 raster layers and a double as input, outputting two layers.

What would be the changes required to the example code to allow that?

For QGIS 2.x I would have used the following syntax :

##Layer1=raster
##Layer2=raster
##myDouble=Double
##OutLayer1=output raster
##OutLayer2=output raster

From what I understand, the changes have to be made in the following procedure, but I'm not sure what to put in place.

def initAlgorithm(self, config=None):
    self.addParameter(QgsProcessingParameterFeatureSource(
        self.INPUT,
        self.tr("Input layer"),
        [QgsProcessing.TypeVectorAnyGeometry]))
    self.addParameter(QgsProcessingParameterFeatureSink(
        self.OUTPUT,
        self.tr("Output layer"),
        QgsProcessing.TypeVectorAnyGeometry))

On May 16th, the QGIS Python API documentation was released. However it is still unclear for me how to use it here. (Which might very well be a lack of Python knowledge)


texts: How to get the closest line to a selected polyline

texts:
polylines A, B, and C meet at point P1. Whats a quick way that I could find out which line is closest to polyline A using arcobjects ?

At first I used IProximityOperator by casting PolylineA and then using ReturnDistance method on the other polylines, but I kept getting a 0 in return. I thought it was because they were connected.

var proxyOpt=(IProximityOperator)polylineA.Geometry;
foreach(var polyline in polylines){
  // I keep getting Zero
  var distance = proxyOpt.ReturnDistance(polyline.Geometry)
 }

Should I try connecting the polylines, and using the angles between the 2 lines instead ?


texts: Using different symbols and colors for each row in a shapefile attribute table

texts:
I have a shapefile where I have a column of individual locations (ID field below), and then a column of yes/no data for each location. Some locations were visited twice and got different results.

I need to display this data so that each location is its own color, and where Yes locations are all triangles (regardless of color) and No locations are all circles.

Is this possible?

I tried it out using the Symbology - Multiple Attributes tool, it'll let me manually set colors for yes or no, but not unique symbols for different locations.

I'm using ArcGIS 10.2


texts: Getting grid points for DEM purpose that is bounded by surveyor points?

texts:
I have a list of surveyor points, and I will need to create a DEM model out of it ( in geotiff format if possible for space consideration, if not, then in ASCII grid format).

I need the DEM for watershed and catchment delineation purpose.

Instead of using existing software like ArcGIS, I want to write my code to do it ( in C++ or C#), so that I can later integrate them into my own GIS software. is there anything in GDAL or other open source tools that do this well? Or do I have to write from scratch?


texts: GeoDjango and limit with bounding box query?

texts:
Just some context first, I have a postgres database with PostGIS extension, running under a GeoDjango wrapper. I have a weird bug that hopefully someone can shed some light on.

I'm running a query like this

geom = Polygon.from_bbox(bbox)
query = Venue.objects.filter(location__within=geom)

This works fine for most cases and is doing exactly what I want, except for in one specific situation where I get a result that looks like this

(for reference, the bbox param above is always defined by the SW and NE corners of the map which in this case are

sw_lat:50.98777649091917
sw_lng:-0.9302881336211613
ne_lat:52.15534947527341
ne_lng:0.06946772575383875

)

I've been trying to narrow down the problem, and I'm starting to think it might be something to do with crossing from negative to positive longitude, but I'm not sure. Or maybe it is something to do with indexing on the database level?

Here is the SQL that is generated:

SELECT * FROM "web_venue" WHERE ST_Within("web_venue"."location", ST_GeomFromEWKB('\001\003\000\000 \346\020\000\000\001\000\000\000\005\000\000\000(\366\304G\026|I@\000(\\\317C\266\362\277(\366\304G\026|I@\000@\341z\036\272\305\277y\006C\216\230\021J@\000@\341z\036\272\305\277y\006C\216\230\021J@\000(\\\317C\266\362\277(\366\304G\026|I@\000(\\\317C\266\362\277'::bytea)) LIMIT 50

(I took out all the explicit select statements and just put '*' for brevity.)


texts: Control view of Multi polygon in Google map

texts:
We have multipolygon WKT. We have no issue to draw all of them in Google map.
Below are the codes. The only issue is that we want only the first polygon to be visible but the rest to transparent. How to control on that? Any idea please?

function drawPoly(data, pl) {
    var pointsData = data.split(",");
    var cords = [];

    // iterate over each points data and create a latlong
    // & add it to the cords array
    var len = pointsData.length;
    for (var i = 0; i < len; i++) {
        var xy = pointsData[i].split(" ");

        var pt = new google.maps.LatLng(xy[1], xy[0]);
        cords.push(pt);
    }

    // now add it to the map
    // Construct the polygon
    alert("I=" + pl);
    var polygonFirst = "";
    if (pl == 0) {
        var polygon1 = new google.maps.Polygon({
            paths : cords,
            strokeColor : '#1E90FF',
            strokeOpacity : 0.8,
            strokeWeight : 2,
            fillColor : '#1E90FF',
            fillOpacity : 0.35,
            zIndex : 1
        });
        // polygon1.setMap(map);
    } else {
        var polygon2 = new google.maps.Polygon({
            paths : cords,
            strokeColor : '#1122AA',
            strokeOpacity : 1,
            strokeWeight : 2,
            zIndex : 1
        });
        // polygon2.setZIndex(1);
        polygon2.setMap(map);
    }

}


texts: Lebanese army projection + datum to WGS84

texts:
Most public offices in Lebanon give their data out as Esri shapefiles, with a strange non-standard projection, making it very difficult to use with any data from other sources, including any data in WGS84.

The .prj file in each of the maps looks like this (broken into multiple lines for comprehension):

PROJCS[
  "SGL_Stereo",
  GEOGCS[
    "GCS_SGL_GEO",
    DATUM["D_SGL",SPHEROID["Clarke_1880_IGN",6378249.2,293.4660212936265]],
    PRIMEM["Greenwich",0.0],
    UNIT["Degree",0.0174532925199433]
  ],
  PROJECTION["Double_Stereographic"],
  PARAMETER["False_Easting",0.0],
  PARAMETER["False_Northing",0.0],
  PARAMETER["Central_Meridian",39.15],
  PARAMETER["Scale_Factor",0.9995341],
  PARAMETER["Latitude_Of_Origin",34.2],
  UNIT["Meter",1.0]
]

When I open such files in QGIS, the following CRS is extracted:

 * Generated CRS (+proj=sterea +lat_0=34.2 +lon_0=39.15 +k=0.9995341 +x_0=0 +y_0=0 +a=6378249.2 +b=6356515 +units=m +no_defs)

The closest to a solution that I found was to reproject to EPSG:22770 (Deir ez Zor / Syria Lambert), and manually adjust the content by a few meters. This is very unsatisfying and ultimately inaccurate.

I want the maps given to me by the municipality in WGS84. How do I do that?


texts:
I want to get path and row number for my area of analysis to procure LISS IV raster images. Where can i find the same.

I understand similar question has been asked before but the solutions given in those post are not working for me.


texts:
I am trying to calculate the habitat fragmentation of a forested region in Africa in order to protect the areas less fragmented. I have used the tool "Landscape Fragmentation Tool V2 requiring Spatial Analyst" but I have an error during its execution. Could you recommended me another tool to calculate fragmentation of habitats using arcgis 10.1?

The input that I used was a land cover map with 1 and 2 data, and the edge width was 50.

This is the error:

    Executing: Landscape_Fragmentation_Tool Bool_fragm 50 C:\CW2\project_2 fragmentation
Start Time: Mon May 12 22:51:18 2014
Running script Landscape_Fragmentation_Tool...

OUTPUT DATASET NAME IS:

    C:\CW2\project_2\fragmentation.img

DATA EXTRACTED FROM LAND COVER
EDGE FRAGMENTATION IDENTIFIED
INTERIOR, PATCH, AND PERFORATED FRAGMENTATION IDENTIFIED
CATEGORIZING CORE FOREST PATCHES...
  File "C:\arcgis_extras\Landscape_Fragmentation_SA\LFT_SA.py", line 254, in <module>
    gp.RegionGroup_sa("core1_", core_grp, "EIGHT", "WITHIN")

<class 'arcgisscripting.ExecuteError'>: ERROR 010213: Error in reading raster c:\fragdata\core1_.
ERROR 010067: Error in executing grid expression.
Failed to execute (RegionGroup).

Completed script Landscape_Fragmentation_Tool...
Failed to execute (Landscape_Fragmentation_Tool).
Failed at Mon May 12 23:33:04 2014 (Elapsed Time: 41 minutes 46 seconds)

Thanks!

UPDATE:

Here you can find more lines in case you have some idea:

 #-------------------------------------------------------------------
# CLASSIFY CORE FOREST BASED ON PATCH SIZE

gp.workspace = TempWS

gp.AddMessage("CATEGORIZING CORE FOREST PATCHES...")
print "CATEGORIZING CORE FOREST PATCHES..."

gp.Reclassify_sa (fragMap2_f, "Value", "4 1", "core1_", "NODATA")

# region group core class...
core_grp = TempWS + "\\core_grp.img"
gp.RegionGroup_sa("core1_", core_grp, "EIGHT", "WITHIN") <-- The error is here

# convert core groups to txt file...
coreGrpFile = "%s\\coreGrpFile.txt" % TempWS
gp.RasterToASCII_conversion(core_grp, coreGrpFile)

# convert 4 class frag map to txt file...
fragMap_4c = "%s\\fragMap_4c.txt" % TempWS
gp.RasterToASCII_conversion(fragMap2_f, fragMap_4c)

cur = gp.SearchCursor(core_grp)
row = cur.next()

remap_dct = {}


texts: ArcGIS Online:Unable to launch application: Unable to access application's configuration

texts:
I have a custom configuration for an ArcGIS Online Web Application, and have followed all instructions.  I have changed the Web App settings and receive this message.
Unable to launch application: Unable to access application's configuration
My JSON file is configured correctly.   What could be the cause of this?


texts: How can I export geotagged photos from Google Earth Pro?

texts: I have a project and I am assessing social media based landscape photographs. How can I download photos as point shapefile and export them from Google Earth Pro to ArcMap?


texts: How would I go from a set of points representing routes to a set of map tiles for online consumption?

texts:
Possible Duplicate:
  Web-based user editable routing application

I have a set of routes as a text file with coordinates.  For each coordinate I also have whether or not it is an intersection with another route.  I'm completely new to GIS and I'd like to know how I can go from these paths to an online map showing the routes.

I'd like to use all free software and would like the resulting tiles to be usable in something like openLayers.  I would also like to perform routing algorithms on the paths later, so I think using a spatial database would be appropriate.  I'm just not sure what steps would be involved and what software should be used.  Also, since I'm completely new, if you can recommend any books, tutorials, or videos that would be great too.  Thanks!


texts: How to copy attributes from many columns to one, QGIS3

texts:
I have many columns, about 30 and each of those has some attributes. I want to merge those columns to one and collect all attributes in a row of columns. I tried with this post but it is not for I am looking for.
Combining Columns QGIS Attribute Table
I have an example in ms excel what I want that without commas and spaces. Attributes have to be separated with one comma.

With this script, if some fields are empty, I have got null fields. Sorry for my English


texts: Filtering a collection the map is not showing the filtered results

texts:
I filter an image collection by a country which seems to works well but when I show the results on the map they appeared not filtered:

var countries = ee.FeatureCollection("USDOS/LSIB_SIMPLE/2017");
var vn_shape = ee.Feature(
  countries
    .filter(ee.Filter.eq('country_na', 'Argentina'))
    .first()
    );
print (vn_shape);
Map.addLayer(vn_shape, {color: 'black'}, 'Country');

var firesAus = ee.ImageCollection('FIRMS')
             .filterBounds(vn_shape.geometry())
             .filterDate('2019-01-01', '2020-01-31');
print('Size Country Fires', firesAus.size())

Map.addLayer(firesAus, {color: '00FFFF'}, 'Fires');
Map.setCenter(-64.03, -31.54, 3)

I wanted the fires appeared only inside the country, what I am doing wrong?


texts: ArcObjects - find nearest feature to a point - has to be FAST

texts:
I've seen lots of discussion on this, especially using IIndexQuery.NearestFeature() on a FeatureIndex instance. Something like this:

private IdDistancePair FindNearest(IGeometry g, IFeatureClass fc, string where = "")
{
IFeatureIndex fi = new FeatureIndexClass();
fi.FeatureClass = fc;

//you can optionally set a FeatureCursor on the IFeatureIndex
//to reduce the number of features indexed.
//also, limiting subfields may speed up the indexing operation
IQueryFilter qfilter = null;
if (!string.IsNullOrEmpty(where))
{
  qfilter = new QueryFilterClass();
  qfilter.WhereClause = where;
  IFeatureCursor fcursor = fc.Search(qfilter, true);
  fi.FeatureCursor = fcursor;
}

//create the index
fi.Index(null, ((IGeoDataset)fc).Extent);

int fid;
double distance;
((IIndexQuery)fi).NearestFeature(g, out fid, out distance);

return new IdDistancePair(fid, distance);
}

However, it appears that the FeatureIndex is an in-memory thing unrelated to any existing spatial index on the query feature class. It takes a long time (dozens of seconds) to build. (and the above doesn't work for me anyway, always returning a FID of -1 after a 45s wait).

I only need to do one query per layer, repeated for ~5 layers. The docs on the FeatureIndex class suggest that the above approach is expensive in such a circumstance.

I need this to complete in say 10 seconds or less. I CAN limit the search to a radius, say a few thousand meters. But there are still likely to be thousands of query features inside that radius, so I'd prefer to avoid any brute-force loop/measure the distance solution in favor of one that uses a spatial index, if possible.

Any suggestions?


texts: How to load NASA JPL WMS Global Mosaic into QGIS

texts:
I was playing around with uDIG and realized that its WMS capabilities are awesome! The Global Mosaic perfectly fits a map I created! However, I was planning to use the WMS layer as an orientation for cleaning up my polygons, yet geoprocessing (merge etc.) is not very much straightforward in uDIG (is it?).

Therefore I wanted to load the WMS into QGIS (1.7.2). After I set up the connection as
  http://wms.jpl.nasa.gov/wms.cgi?

the connection works but when loading any service I get the following response:

Antwort: <?xml version='1.0' encoding="UTF-8" standalone="no" ?>
<!DOCTYPE ServiceExceptionReport SYSTEM "http://wms.jpl.nasa.gov/exception_1_1_1.dtd ">
<ServiceExceptionReport version="1.1.1">
<ServiceException>
This server no longer provides full WMS services!
</ServiceException>
</ServiceExceptionReport>

can anyone suggest a solution for correctly loading this WMS into QGIS?
thank you!

Jens


texts: Is there a way to automatically find errors made by the Watershed tool?

texts:
I've delineated the watershed of near 7000 reservoirs with the Watershed tool in ArcGIS 10.1 and converted them in polygons. As you can see in the image bellow, some mistakes were made because the pour point didn't snap properly.

Is there a way to automatically find these errors in my Shapefile so I can redo them?


texts: Reading page size from MapDocument object?

texts:
I have a set of 30 MXD's with different page size. I would like to export maps with their name file which include the mxd's page size. i have this code:

import arcpy,os,sys
import arcpy.mapping
from arcpy import env

env.workspace = r"C:\Project"
for mxd in arcpy.ListFiles("*.mxd"):
    mapdoc = arcpy.mapping.MapDocument(r"C:\Project\\" + mxd)
    df = arcpy.mapping.ListDataFrames(mapdoc, "Layers")[0]
    arcpy.mapping.ExportToJPEG(mapdoc,r"C:\Project\out\\"+ os.path.basename(mapdoc.filePath).replace ('.mxd','') +".jpg")
    print 'ExportToJPEG'
    mapdoc.save()
    del mxd

the row that start with "arcpy.mapping " is Connected with the next row.

How can I do this?


texts: ArcSDE installation; "sdesetup: command not found"

texts:
I am trying to install ArcSDE component for postgresql database from A DVD file.

I have copied files of following folders on the server :

command_64
pg_64

in both of theses files there are install files. for both I used :

./install -load
yes
mount point: DVD location or wherever you are installing from
install directory: enter to accept /home/sde
Package numbers to load: enter to accept all
yes
no

now I have sdeexe101 file which contains bin, etc, liband locale folders in it.

I also  As the root user, copied the st_geometry.so, libsg.so, and libpe.so files from $SDEHOME to the lib directory in the postgres location.

[root@machine pgsql]# cp /home/sde/sdeexe101/lib/st_geometry.so /usr/lib/pgsql
[root@machine pgsql]# cp /home/sde/sdeexe101/lib/libsg.so /usr/lib/pgsql
[root@machine pgsql]# cp /home/sde/sdeexe101/lib/libpe.so /usr/lib/pgsql

when I tried to execute sdesetup to create geodatabase and authorize it using:

sdesetup -o install -d POSTGRESQL -u sde -p sde...

I got this:
sdesetup: command not found

what am I missing here? Are pg_64 and commands_64 files that I needed?

I appreciate any help!


texts: Error in QGIS 2.2 Raster Merge Tool in OS X

texts:
I'm having this strange error in QGIS 2.2 for OS X while trying to run the raster merge tool (in the Raster - Miscellaneous - Merge). This happens just by clicking on the "Select" button to choose the files that will be merged.

Does anyone has an idea of what can be the problem? Is it due to the Python version I have installed

Thanks!

An error has occured while executing Python code:

Traceback (most recent call last):
  File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/GdalTools/tools/doMerge.py", line 92, in fillInputFilesEdit
    files = Utils.FileDialog.getOpenFileNames(self, self.tr( "Select the files to Merge" ), Utils.FileFilter.allRastersFilter(), lastUsedFilter )
  File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/GdalTools/tools/GdalTools_utils.py", line 394, in getOpenFileNames
    return self.getDialog(parent, caption, QFileDialog.AcceptOpen, QFileDialog.ExistingFiles, filter, selectedFilter, useEncoding)
  File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/GdalTools/tools/GdalTools_utils.py", line 351, in getDialog
    dialog.selectNameFilter(selectedFilter[0])
TypeError: QFileDialog.selectNameFilter(QString): argument 1 has unexpected type 'QPyNullVariant'

Python version:
2.7.5 (default, Aug 25 2013, 00:04:04)
[GCC 4.2.1 Compatible Apple LLVM 5.0 (clang-500.0.68)]

QGIS version:
2.2.0-Valmiera Valmiera,

Python path: ['/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/processing', '/Users/gds506/.qgis2/python/plugins/LecoS', '/Applications/QGIS.app/Contents/MacOS/../Resources/python', u'/Users/gds506/.qgis2/python', u'/Users/gds506/.qgis2/python/plugins', '/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins', '/Library/Frameworks/SQLite3.framework/Versions/B/Python/2.7', '/Library/Frameworks/GEOS.framework/Versions/3/Python/2.7', '/Library/Python/2.7/site-packages/numpy-override', '/Library/Python/2.7/site-packages/matplotlib-override', '/Library/Frameworks/GDAL.framework/Versions/1.10/Python/2.7/site-packages', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python27.zip', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-darwin', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-mac', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-mac/lib-scriptpackages', '/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-tk', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-old', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-dynload', '/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/PyObjC', '/Library/Python/2.7/site-packages', '/Library/Python/2.7/site-packages/PIL', '/Users/gds506/.qgis2/python/plugins/mmqgis/forms', '/Applications/QGIS.app/Contents/Resources/python/plugins/fTools/tools']


texts:
I have a method that creates a new Shapefile layer and returns the layer.  As soon as I try to operate on that layer, I get a Segmentation fault.  I verified the code works by pulling all the code out of the method and sticking it in the calling method.  Thoughts on why this happens?  This is a different problem than I posted that ended up being a PIP/Conda conflict.

def CreateMergedLayer(merged_layer_file_name):

  """
  Initializes an ESRI Shapefile with the fields of two
  layers that will be merged - could be improved by creating
  a set of fields from both (assuming they have the same
  definition if sharing names), but just made purpose-built
  for this function - if i have to do again, i'll refactor

  Parameters
  ----------
  merged_layer_file_name : string
    The full path and file name of the new merged layer

  Returns
  -------
    The pointer to the layer
  """
  print ("About to create shapefile for {}".format(merged_layer_file_name))

  driver = ogr.GetDriverByName("ESRI Shapefile")
  data_source = driver.CreateDataSource(merged_layer_file_name)

  # create the spatial reference, WGS84
  srs = osr.SpatialReference()
  srs.ImportFromEPSG(4326)

  new_layer = data_source.CreateLayer("la_streets_with_block_centroid_connectors", srs, ogr.wkbLineString)

  new_field = ogr.FieldDefn("GeoID", ogr.OFTString)
  new_field.SetWidth(16)
  new_layer.CreateField(new_field)

  new_field = ogr.FieldDefn("FULLNAME", ogr.OFTString)
  new_field.SetWidth(100)
  new_layer.CreateField(new_field)

  return new_layer

Calling this method here:

merged_layer = CreateMergedLayer("/tmp/la_streets_with_block_centroid_connectors.shp")

print ("Created Layer, getting definition: layer: {}".format(merged_layer))
merged_layer_defn = merged_layer.GetLayerDefn()

This outputs the following without getting past the GetLayerDefn() request:

About to create shapefile for /tmp/la_streets_with_block_centroid_connectors.shp
Created Layer, getting definition: layer: <osgeo.ogr.Layer; proxy of <Swig Object of type 'OGRLayerShadow *' at 0x114ad4630> >
Segmentation fault: 11


texts: Create a flow vector NetCDF from X and Y velocity components

texts:
I want to create a flow vector NetCDF data from two .tif raster, one representing X-velocity and another representing Y-velocity component. Both raster groups have identical spatial extent. I have converted each raster to NetCDF in ArcGIS and then merged it using NCO tool (ncecat). I had hoped I can symbolize the resulting NetCDF data as a vector, using arrows in QGIS but it is not recognized as a vector layer.
So, my question is how can I create a flow vector NetCDF data from two velocity components rasters (or from velocity magnitude and flow direction rasters) using NCO tools or Python? Is there any other way to do it beside these two?


texts: Iterating over shapefile to mask and extract raster data in Python

texts:
I have a shapefile that is a grid with multiple polygons that have a specific id in them ranging from 1 to 50. Each polygon overlays an agriculture research plot that I would like to calculate vegetation indices. I then want to append this information to an Excel file.
But I am having trouble with getting the loop to work, where only the first polygon appends data to the Excel sheet, it seems like it doesn't even loop through. Any ideas or improvements?

import rasterio
import fiona
import numpy as np
import shapefile as shp
from openpyxl import load_workbook

with fiona.open("path/to/shapefile.shp", "r") as shapefile:
    shpfl = shp.Reader("path/to/shapefile.shp")
    recs = shpfl.shapeRecords()
    count = 0
    wb = load_workbook(wb_output)
    ws = wb.active
    for feature in shapefile:
        shape =[feature["geometry"]]
        plotNum = recs[count].record[2]
        sheetRef = plotNum + 1

        #open and mask the multispectral image to the plot
        with rasterio.open(multispec) as src:
            out_image, out_transform = rasterio.mask.mask(src, shape, crop = True)
            out_meta = src.meta

            #extract color bands
            red = out_image[0,:,:].astype('float64')
            nir = out_image[3,:,:].astype('float64')
            red_edge = out_image[4,:,:].astype("float64")

            #ndvi
            num = np.subtract(nir, red)
            denom = np.add(nir, red)
            np.seterr('ignore', 'ignore')
            ndvi = np.divide(num, denom)

            flatndvi = ndvi.flatten()
            avg_ndvi = np.nanmean(flatndvi)

            cell = ws['c' + str(sheetRef)]
            cell.value = avg_ndvi
            wb.save(wb_output)


texts:
I am a new MODIS Level 1B user. I have been using the MODIS Swath tool to process my images.

However, I have over 1000 images to process. I was wondering if there a script available in IDL or R (very little coding experience) where I can batch process all the images?


texts: Load information about feature for Geojson

texts: Is it possible to display whole "info" for vector layer? I have a lot of geojsons with different feature in one map. I would like to show information.WMS can display info in this case (image). How I should do this for geojson?


texts: Points freezed in OpenLayers when adding data to different layer

texts:
I have a map using OpenLayers with 2 points vector layers ( layer A  & Layer B ) , layer A already contains hundreds of points and navigating well and layer B without features, whenever I add new features to layer B, features will be added without any problem but when I try to zoom in & zoom out the map, I have found that points in layer A freezed and not in correct places and I got the below script errors in OpenLayers.js file.

Uncaught NotFoundError: Failed to execute 'insertBefore' on 'Node': The node before which the new node is to be inserted is not a child of this node

OpenLayers version 2.12


texts: Filter diagrams with atlas coverage layer in QGIS

texts: Please advise how I can hide diagrams outside of the coverage layer when using QGIS atlas function. I can do this for polygons,rasters etc without problem but whatever I do I cannot hide diagrams from other layers. Please see example below. Basically I do not want anything outside of the polygons..


texts:
I know how to produce a map with the Mollweide projection, both using QGIS and R (I'm most interested in doing all of this in R). The proj4string that I've used in R is +proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs

But I am interested in producing a map using what I've seen called the Mollweide oblique projection, or the Atlantis projection, or alternatively the "Hammer map". See here for examples.

I'm not quite sure where to start, though.

Does this involve adjusting some of the parameters in the above Mollweide proj4string?

For starters, here is how I get a map in Mollweide projection:

library(maptools)
data(wrld_simpl)

mollCRS <- CRS('+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs')

mollMap <- spTransform(wrld_simpl, mollCRS)


texts: I'm new to GIS and I'm using QGIS 3.4.5. I've been downloading individual SRTM tiles containing DEM data and merging them together in order to see the elevation of every point in a given area (such as the entire USA). It's slow and tedious to download all these files individually. Is there some XYZ tile server or other service I can connect to that will serve the DEM data for wherever my GIS canvas is located? Just like how I can add a Google Maps layer or a hillshade layer.


texts: Python script terminates after calling QgsProject.write()

texts:
Python script using QGIS package terminates suddenly when calling QgsProject.write.
Idea is to create QGIS Project with select list of files and save it to disk, in automated manner.

Script:

from qgis.core import QgsProject, QgsRasterLayer

def foo():
    filePath = "C:/Users/username/Documents/projects/GeoTiffSample.tiff"
    savePAth = "C:/Users/username/Documents/projects/test2.qgs"
    proj = QgsProject.instance()
    proj.setFileName(savePAth)
    testlayer = QgsRasterLayer(filePath, baseName = "testLayer")
    proj.addMapLayer(testlayer)
    proj.setFileName(savePAth)
    proj.setDirty(True)
    try:
        print("Prewrite")
        ret = proj.write()
        print("Postwrite")
        pass
    except Exception as exc:
        print(exc)
    print(ret)

print("START")
try:
    foo()
except Exception as exc:
    print(exc)
print("END")

Output:

START
Prewrite

Documentation is kind of lacking on any info. I suspect unhandeled exception in QGIS library.


texts: leaflet layerGroup is rendered behind baselayer

texts:
I have several baselayers and overlays. Everything is fine except the lauerGroup is rendered behind the baselayer.
First I define the baselayers and add them to a 'baselayer' variable.

// base layers
naturalEarth = L.tileLayer('http://{s}.tiles.mapbox.com/v3/mapbox.natural-earth-2/{z}/{x}/{y}.png'),

grump = L.tileLayer('http://localhost:8000/v2/gdal_urban_rural_pop/{z}/{x}/{y}.png');

var baselayer = {
    'Natural Earth' : naturalEarth,
    'GRUMP Population' : grump
};

Then I do the same for some overlays:

// overlay layers
osm = L.tileLayer('http://{s}.tiles.mapbox.com/v3/mapbox.world-bank-borders-en/{z}/{x}/{y}.png'),

df_admin0 = L.tileLayer('http://localhost:8000/v2/dwelling_fractions/{z}/{x}/{y}.png', {
        wax: 'http://localhost:8000/v2/dwelling_fractions.json'
    }),

df_port = L.tileLayer('http://localhost:8000/v2/PRT_dwelling_fractions/{z}/{x}/{y}.png', {
        wax: 'http://localhost:8000/v2/PRT_dwelling_fractions.json'
    });

//add layer group
hazard_curve_group = L.layerGroup();

hazard_curve1 = L.tileLayer('http://localhost:8000/v2/hazard_curve/{z}/{x}/{y}.png', {
        wax: 'http://localhost:8000/v2/hazard_curve.json'
    }).addTo(hazard_curve_group);

hazard_curve2 = L.tileLayer('http://localhost:8000/v2/hazard_curve2/{z}/{x}/{y}.png', {
        wax: 'http://localhost:8000/v2/hazard_curve.json2'
    }).addTo(hazard_curve_group);

var overlays = {
    "OpenStreetMap" : osm,
    "Dwelling Fractions Admin 0" : df_admin0,
    "Dwelling Fractions Portugal" : df_port,
    "Hazard Curves test" : hazard_curve_group,
};

// create the map
var map = L.map('map', {
    ...
});

//add it all to the map
L.control.layers(baselayer, overlays, {collapsed: true, position: 'topleft'}).addTo(map);


texts: Raster Layer Statistics CSV Output using QGIS Graphic Modeller

texts:
I am currently working on a project where I am using historical cloud coverage from MODIS Cloud Fraction data to try and estimate future cloud coverage of any particular month. I currently have 10-years worth of data for January through to December and have generated global mean, minimum and maximum raster layers (Step 1) for each individual month for that 10-year period. I have detailed the methodology in the flow diagram below to show the steps taken. I currently have 3 layers of global statistics (mean, minimum and maximum) for each month of the 10-year period (Raster 2).

As I am doing this is in Graphic Modeller, is there any way that the final output in Step 4 be written to a CSV file? Conversely, is there a way to script this process (I don't code well) that will make things a bit easier/ more specific if I want to specify output statistics such as minimum, maximum, mean, +1 standard deviation and -1 standard deviation? This is all so that I can display (at least) the minimum, maximum and mean cloud cover from January to December for any given AOI on a graph.


texts: Overlaying a mbtiles tile layer produced with tilemill over another layer maintaining transperency

texts:
So using leaflet, I'm overlaying my base tile layer as google satellite imagery. I want to overlay another layer (Mbtiles layer) produced with tilemill over this satellite layer. I was successful in doing this. I set the background of the mbtiles file as transparent on tilemill. Im using tilestache to serve the tiles and leaflet for a JS API.

Now, I want this mbtiles layer to only show the coloured layers (roads, polygons etc.) and the background to be transparent (to show the google satellite layer wherever the mbtiles layer is transparent). But this is not happening. The background of the MBtiles layer is white and there is no transparency. Where am i going wrong?


texts: Approximating GPS coordinates for a point from n expected locations

texts:
I'm trying to find the coordinates of my point of attention (point X,marked blue). I used the GPS device of my car to collect the coordinates according to where I parked my vehicle each time I visited point x. So after following this exercise for 16 days, I managed to obtained 16 sets of coordinates, spread round about my point of attention.

After plotting these coordinates on the map, I observed the following:
 Twice or thrice out of ten times, my GPS device gave a wrong set of  coordinates which turned out to be quiet far from point X. Also due to traffic, occasionally I'm unable to park close to point x and hence in this case too, the coordinates obtained are far from point X.

Problem : From the 16 sets of coordinates obtained, what process do I use to narrow down to one set of coordinates which is in close proximity to my point of attention (point X)?


texts: Landsat 7 - 8 surface reflectance bands from ESPA (LEDAPS)

texts:
I have downloaded some Level 1 and higher level imagery data including surface reflectance as well.

I should process these data in Idrisi (Selva) package.
The Level 1 images' pixel values are between 0-255. The surface reflectances' values are between 0-10000 (valid range, according to the XML metadata documentation). Unfortunately Idrisi can not interpret these values correctly. The Level 1 (DN) composits look just fine, whilest the composits (from same bands) of surface reflectance look very weak (pale). I aware the 0.0001 scale of SR, but its application did not help.

Please advise how could I overcome on this issue (might be by image enhancements, or reclass, or some transformation?).


texts: Is it possible to get satellite image tile sets from something like openmaptiles tiles offers?

texts:
I am able to get an offline tile set for the USA from here. However, I was wondering if it is possible to get satellite imagery of those tile sets. I have been checking online and most require to be connected to the internet and cache for 30 days then reconnected.

I would like to download satellite imagery tile set that I could just keep for like a year or so then download later.

I am currently using tileserver-gl to host the mbtiles tile set for offline mode.

Restrictions involved:

Javascript (Leaflet)
tileserver-gl requires .mbtiles


texts: Wavelet vs Discrete compression for imagery

texts:
What is the main reasons to choose wavelet (jp2, sid, ecw) instead of discrete (JPEG) compression methods for raster data?

Lets limit this question to aerial imagery and lossy compression.

Why should we give preference to wavelet compression instead of "classic" JPEG? It will give smaller size but may be there any drawbacks, what do you think?  Will it lead to more CPU processing for decomression (during viewing)?


texts: 'Locking' intervals in interval determination

texts:
Working with graduated maps in Qgis 2.8 I'd like to lock a certain interval value while letting Qgis determine the rest - in relation to eg. missing values.

I know how to work with queries, but this makes Qgis not display certain values at all, where as i would like for it to display these values but disregard them in interval determination. Query filtering furthermore requires manual addition of the empty interval class whereas the method I'm looking for automates this by "always" having the class interval locked.


texts: Automating ArcMap with ArcObjects for unit testing?

texts:
I need to run unit tests against ArcMap extension which closely coupled with specific edit workspace and map state. When I try to initialize ArcMap 10 application with code like this:

ESRI.ArcGIS.RuntimeManager.Bind(ESRI.ArcGIS.ProductCode.Desktop);
IDocument pdoc = new MxDocumentClass();
_application = pdoc.Parent;

I receive error:

  The runtime application type must be specified before license initialization

Is there any workaround or better way for unit testing?


texts: What WMS server would be suited for serving many aerial photos?

texts:
Given a big set of aerial photos taken. 7000x5000 pixels.

What is the best way to serve these in a WMS. When reading up about different WMSs, I get confused about the complexity and have a hard time finding the simple answer to how we can serve large about of images.

One of the confusing things are raster layers. I assume we need to serve raster layers since our original data is aerial photos from a plane. My readings seeem to lead me on the path that a raster layer is one big picture which the WMS then slices up. But its really not possible to stich all the aerial photos together to one big picture. The other path is to create many layers with one aerial photo in each and therefore have a raster layer per image, which also seems stupid.

Is there any WMS implementations out there where we can hook our images into such when a request comes in for a BBOX, then we tell it what image it should use and then it slices it up and not having to have one big image in memory with the hole world. (An answer to all this could be, how do openstreetmap create its internal map?)


texts: Is there a way to edit field name, type & length after importing from Excel?

texts: Is there a way to edit the name, type and length given to each column on the attribute table once the information has been imported from Excel?


texts: Saving vector layer with joins converts joined attributes from double to integer in QGIS?

texts: I'm using QGIS 2.16.3. I'm trying to use the "Dissolve with Stats" plugin, which requires me to save the layer with the joined attributes. However, saving the layer converts the joined fields from double to integer, and I need to preserve that precision.


texts:
Does anyone have any tutorial, link, video etc on how I can use Open Street Map for routes?
I know google maps and saw that I can easily do that. But this application will have an internal module and we can not afford the license.


texts:
I don't understand why I should choose readOGR() over readShapePoly.

What are all the different ways of doing this, and why should I choose one over another?


texts: FUSION's PolyClipData creates corrupt LAS file?

texts:
I am trying to learn LASPY and am following a tutorial. I got an error when starting up using:

import laspy
infile = laspy.file.File("./laspytest/data/simple.las", mode="r")

where I got a storage error:

LaspyException: Error mapping file: [Error 8] Not enough storage is available to process this command

To try and get around this I tried FUSION's polyclipdata routine to parse the LAS file into 4 parcels:

command = 'polyclipdata /multifile ' + temp + '\\pointfile2.shp ' + temp + '\\pointfile3.shp {}'.format(in_las_file)
process = subprocess.Popen(command, stdout=subprocess.PIPE, stderr=subprocess.PIPE, shell=True)
(stdoutdata, stderrdata) = process.communicate()
if stderrdata:
   print stderrdata

When I try and work with LASPY again, I get the following error:

ValueError: buffer size must be a multiple of element size

I thought that this was a LASPY error and created a question to resolve it, however, after testing "clean" data, LASPY was determined to not be the problem (laspy ValueError: buffer size must be a multiple of element size). I am using Python 2.7.

Why is FUSION's polyclipdata not producing a usable product?

I started looking into R's lidr package and ran the readLAS command. I got the following messages:

ERROR: 'end-of-file' after 8196120 of 8196121 points
Warning messages:
1: Invalid file: header states the file contains 8196121 points but 8196120 were found. Header has been updated.
2: Invalid file: header states the file contains 7570432 1st returns but 1028359 were found. Header has been updated.
3: Invalid file: header states the file contains 482368 2nd returns but 979894 were found. Header has been updated.
4: Invalid file: header states the file contains 114533 3rd returns but 936803 were found. Header has been updated.
5: Invalid file: header states the file contains 25014 4th returns but 997835 were found. Header has been updated.
6: Invalid file: header states the file contains 3477 5th returns but 1053020 were found. Header has been updated.
7: 7369979 points below 0 found.
8: Dataset is invalid: 1049976 points with a return number of 0 found.

It looks as though the datasets created by FUSION messed up the header?


texts: Different band file format when downloading Landsat 8 from earthexplorer

texts:
When downloading Landsat 8 OLI/TIRS Land Surface Reflectance from the USGS earthexplorer service, the usual archive has the following files in it :

LC82330672013150LGN00_B10.TIF
LC82330672013150LGN00_B11.TIF
LC82330672013150LGN00_B1.TIF
...
LC82330672013150LGN00_B9.TIF
LC82330672013150LGN00_BQA.TIF
LC82330672013150LGN00_MTL.txt

For some scenes (e.g. LC82270652013236LGN00) though, we get the following

LC82270652013236LGN00_cfmask_conf.tif
LC82270652013236LGN00_cfmask.tif
LC82270652013236LGN00_sr_band1.tif
...
LC82270652013236LGN00_sr_band6.tif
LC82270652013236LGN00_sr_band7.tif
LC82270652013236LGN00_sr_cloud.tif
LC82270652013236LGN00.xml

The data stored in the .tif files is very different. When the _MTL.txt file is present, we have band 3 values from 6000 to 10000 for a scene in the amazonian forest. But when the .xml file is present, the values range from 0 to about 1500 (also for an amazonian forest scene which should have similar histogram).

All of the Landsat 8 information I could find assumes we get the _MTL.txt file and I was unable to find any information about the other one. There doesn't seem to be any option in earthexplorer to control which type of data you want to download.

Does someone know how to convert between the two ?


texts: FeatureLayer selectFeatures hides the feature only on one layer

texts:
I am using arcgis 3.23. I have defined the following SimpleLineSymbol to use to highlight selected features in different FeatureLayers.

              app.selectionLineSymbol = new SimpleLineSymbol();
              app.selectionLineSymbol.setWidth(15);
              app.selectionLineSymbol.setColor(new Color([255, 255, 255, 0.60]));

The layers in question are defined as follows

      var blazeColorsTa = new FeatureLayer(app.blazeColorsTaUrl, {
          id: "blazeColorsTa",
          outFields: ["*"]
      });
      blazeColorsTa.hide();
      layers.push(blazeColorsTa);

      var longTrails = new FeatureLayer(app.longTrailsUrl, {
          id: "longTrails",
          outFields: ["*"]
      });
      layers.push(longTrails);

I use the following code to highlight selected features from the blazeColorsTa layer and the selected features show through the selection symbol as desired.

                var blazeColorLayer = app.map.getLayer("blazeColorsTa");
                blazeColorLayer.show();
                blazeColorLayer.setSelectionSymbol(app.selectionLineSymbol);
                blazeColorLayer.selectFeatures(queryParams,FeatureLayer.SELECTION_NEW);

However, when I use identical code to highlight selected features from the longTrails layer the selected features are hidden by the selection symbol.

                var longTrailLayer = app.map.getLayer("longTrails");
                longTrailLayer.show();
                longTrailLayer.setSelectionSymbol(app.selectionLineSymbol);
                longTrailLayer.selectFeatures(queryParams,FeatureLayer.SELECTION_NEW);

I assume the problem has to do with a difference in how the layers are created in the GIS service but the only difference I can see is the longTrails layer is has both the Min Scale and Max Scale set to 0.  Where should I look to see what is causing the issue?

EDIT:

Here is what the trail looks like before selection:

Here is what it should look like after selection.  I changed the color to light green so you could see it better.  The selection symbol works when it is called by the popup.  It also looks like this when I use my selection code on the blazeColorsTa layer.

Here is what it looks like when I call it with the code I listed above on the longTrails layer where the feature is hidden by the selection symbol.

Just in case it is helpful, here is the code which defines the popup mentioned earlier.

            popup = new Popup({
                  lineSymbol: app.selectionLineSymbol,
                  visibleWhenEmpty: false,
                  hideDelay: 0,
                  zoomFactor: 1,
                  keepHighlightOnHide: true
              }, domConstruct.create("div"));

I had to make a custom zoom for the popup and I use the selection code to highlight the feature I have zoomed to.


texts: Doubles are rounded to integers when leaving editing mode in a shp

texts:
I'm very new to QGIS. I have been trying to calculate areas of interception and overlapping percentages of different shp polygon layers. I've watched tutorials and have managed to do it well. However, there's an issue that I havn't been able to solve:

For calculating areas in a polygon layer, I go to the field calculator and input "$area", selecting the "real" type. Clicking OK enables edit mode and propperly calculates areas as doubles. However, when I exit edit mode by clicking on the pencil and saving, all areas are rounded to the nearest integer. The datatype of the attribute is also changed to Integer64

My question is why does this happen and how do I stop it?

I'm working with QGIS 2.16(also tested it with 2.14) and for the following example I used the "continents.shp" layer in the world data from the QGIS official tutorial.


texts:
I've received a CSV file exported from an unknown application containing polygon coordinates in the following format:

Title | Latitude 1 | Longitude 1 | Latitude 2 | Longitude 2 | ... | Latitude 50 | Longitude 50

How can I generate polygons from this data in Mapinfo Professional (12.0)?


texts: Keeping Properties When Converting JSON to SHP

texts:
I'm  trying to convert the Nielsen map of DMAs from a JSON to a shapefile. The properties of the json is as formatted like this:

  "id": 546,
          "properties": {
            "name": "dma:",
            "latitude": 33.896366,
            "tvperc": 88.6,
            "dma": 546,
            "dma1": "Columbia, SC",
            "cableperc": 50.5,
            "adsperc": 38.8,
            "longitude": -80.9102195

However, when I drag the json to QGIS, the ID and the name column are the only thing that are retained. Similarly, when I try to convert this to a shapefile through ogr2ogr, the same thing occurs. Is there a way I can keep the properties? I know I can use the -select parameter but there are no examples on how to use this anywhere in the documentation. My code is below.

ogr2ogr -f "ESRI Shapefile" nielsentopo.json "usa-dma.shp"


texts: different zoom level for same shapefile arcmap

texts: I have a road Shapefile in ArcMap which includes three different type of roads highway, motorway, streets. I want to display all these roads at different scale. One way to do this is to create three new shapefiles and set zoom level for each of these but this is not what i want. I heared that ArcMap modelbuilder can do this but can't find any help.


texts: Upgrading to ArcGIS 10.3 geodatabase on windows R2 SP2 and Oracle 10g (64 bit) 10.2.0.5?

texts:
We have Oracle 10g (64 bit) 10.2.0.5 on Windows Server 2008 Enterprise SP2, and we are planning to migrate our geodatabase from 9.3.1 to 10.3.

As per Esri "Oracle database requirements for ArcGIS 10.3.x" the maximim OS is SP1.

Does anyone know what is limitation/risk/impacts if we keep our OS SP2 not SP1 and migrate our geodatabase?

And what if we migrate our ArcGIS Server/Desktop to 10.3 and kept using geodatabase 9.3.1?


texts: Extract values from raster by using polygons

texts:
I am working with multiple raster files that cover all of united states. The file format of the rasters - GeoTIFF. I have around 100 geotiffs that I need to process. Then I have another shapefile with polygons. What I need to do is extract information (SUM of all values) from the raster files by using the polygon shapes and get data in a table format. I have tried ArcMap 10.3 Spatial Analyst extension's Zonal Statistics (Table) for the SUM, but I was not satisfied with the results - there are way too many 0s in places where there should be some values even if they are small.

I am asking for a good way to extract values in ArcMap, QGIS, possible R or other software. The GeoTIFFS dimensions are 3600 pixels to 859 pixels and they are around 25mb in size (each).

If there any model builder suggestions with functions information that would be great!


texts: Use Google Earth Engine to stack NDVI timeseries layers into one image and export

texts:
I would like to export the NDVI with the date in the filename in one image. Now I can only export them respectively.
And when I export the NDVI image, it shows the error "Error: Exported bands must have compatible data types; found inconsistent types: Int16 and Byte." This is the code:

//study area
var roi = /* color: #d63000 */ee.Geometry.Polygon(
        [[[108.04782666015626, 22.81419793085548],
          [108.70837963867189, 22.829387337342464],
          [108.69876660156251, 23.183323796154976],
          [108.05194653320314, 23.16438649072554]]]);

var geometry = /* color: #d63000 */ee.Geometry.Point([108.50995736836717, 23.048234117697987]);

//cloud mask
function maskL8sr(image) {
// Bits 3 and 5 are cloud shadow and cloud, respectively.
 var cloudShadowBitMask = (1 << 3);
 var cloudsBitMask = (1 << 5);
// Get the pixel QA band.
 var qa = image.select('pixel_qa');
// Both flags should be set to zero, indicating clear conditions.
 var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
                 .and(qa.bitwiseAnd(cloudsBitMask).eq(0));
return image.updateMask(mask);
}

//dataset
var L8 = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
.map(maskL8sr);

var rgb_vis = {min: 0, max: 3000, bands:['B4', 'B3', 'B2']};
var falsecolor = {min: 0, max: 3000, bands:['B5', 'B4', 'B3']};
var filtered = L8.filterDate('2018-01-01', '2018-12-31')
.filterBounds(geometry);
//Map.addLayer(filtered, falsecolor, 'falsecolor');

//add ndwi
var addNDVI = function(image) {
  return image.addBands(image.normalizedDifference(['B5','B4']).rename('NDVI'));
};
var imagen2=L8.map(addNDVI);

var NDVI = imagen2
print(ee.Image(NDVI.first()));

var ndvifilter = NDVI.filterDate('2018-01-01', '2018-12-31')
.filterBounds(geometry);

// Visualize NDVI
var ndviPalette = ['FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718',
  '74A901', '66A000', '529400', '3E8601', '207401', '056201',
  '004C00', '023B01', '012E01', '011D01', '011301'];
Map.addLayer(ndvifilter.select('NDVI'),
            {min:0, max: 1, palette: ndviPalette}, 'ndvi');

//List NDVI
var list=ndvifilter.toList(200);
    for (var i=0;i<200;i++){
         var ndvibatch=ee.Image(list.get(i));
         var date = ndvibatch.date().format('yyyy-MM-dd').getInfo()
         var name= 'NDVI_'+date
         print(name);
    // Export NDVI
         Export.image.toDrive({
               image: ndvibatch,
               description: name,
               fileNamePrefix: name,
               scale: 30,
               region:roi,
               maxPixels: 1e13
               });
    }

The Code link: https://code.earthengine.google.com/b6f6cf802ac8a6667ec20ca6dedffce7


texts: I have a series of points that represent oil and gas wells (over 1,000), and an emissions value for each well. I created 3 concentric buffers that represent increasing distance from the well, and assigned a fraction of the emission value to each buffer. There are many instances where wells are close together, and their buffers overlap. In these overlapping areas, I would like the emission values of the buffers to add. I tried using the raster to polygon tool, but the overlapping emissions did not sum.


texts: How to exclude the area of intersecting polygons from the total sample?

texts:
machine translation

I am interested in solving such a problem. There is a table

id | coordsPolygon
 1 | 0103000020E6100 ...
 2 | 01030000301F300 ...
...

I have to make a query that returns the total surface area covered by polygons. However, there are polygons that intersect with each other, or enter into each other, which adds to the usual query type
sum (ST_Area (coordsPolygon)) excess area.

Question: how can I get the total surface area covered by polygons?

интересует решение такой проблемы. Есть таблица

id|coordsPolygon
 1 |0103000020E6100...
 2 |01030000301F300...
...

Я должен произвести запрос, который возвращает общую площадь поверхности,  покрываемой полигонами. Однако есть полигоны, которые пересекаются друг с другом, либо входят друг в друга, что добавляет к обычному запросу типа
sum(ST_Area(coordsPolygon))
лишнюю площадь.

Вопрос: как я могу получить общую площадь поверхности,  покрываемой полигонами?


texts: Converting a PyQGIS script that works for one input to work with multiple inputs

texts:
A script I have written to perform a number of tasks on a single input needs amending to work with multiple inputs given an input folder. The working script is this and takes a input raster through some processes to output a shp:

##in_rast=raster
##out_shp=output vector

#Take grid and create extent
out1=processing.runalg('gdalogr:rastercalculator', in_rast,'1',None,'1',None,'1',None,'1',None,'1',None,'1','A>0',None,5,None,None)
#Ratser extent to polygon
out2=processing.runalg('gdalogr:polygonize', out1['OUTPUT'],'DN',None)
#Calculate area m^2 and add to new field 'Area_m2'
out3=processing.runalg('qgis:fieldcalculator', out2['OUTPUT'],'Area_m2',0,10.0,3.0,True,'$area',None)
#Select polygons >5m^2 and extract to new shp
out4=processing.runalg('qgis:extractbyattribute', out3['OUTPUT_LAYER'],'Area_m2',3,'5',None)
#Fill holes less than 200m^2
out5=processing.runalg('qgis:fillholes', out4['OUTPUT'],200.0,out_shp)

I have made a start on updating this to work on multiple inputs but being relativity new to python my understanding of the correct methods to use has reached its limit. The updated script looks as follows, only working and outputting for the first process:

##Input_raster=multiple raster
##result_ras=folder

import os

layers = Input_raster.split(';')
for raster in layers:
    result_name = os.path.splitext(os.path.basename(raster)) [0] + "_Extent.tif"
    out1=processing.runalg('gdalogr:rastercalculator', raster,'1',None,'1',None,'1',None,'1',None,'1',None,'1','A>0',None,5,None,result_ras + "/" + result_name)

The question is; how can I continue to write the rest of the processes into my script? I am only interested in the output for the final step so all intermediate outputs can be saved as temp. Ultimately I'd like the output to take the name of the original input raster and append _Outline.shp.


texts: How to query a JTS QuadTree using a query envelope with holes

texts:
I have a Java Topology Suite QuadTree containing Points, now I'm trying to query the QuadTree with a polygon that contains a hole.  Here's how I create the polygon.

//Method will create a searchable Polygon in the shape of a donut.
public Polygon searchAreaPoly(double outterBoundry, double innerBoundry) {

    Coordinate center = MyIndex.getWGSCoord(longitude, latitude);
    GeometricShapeFactory gsf = new GeometricShapeFactory();

    gsf.setCentre(center);
    gsf.setNumPoints(20);
    LinearRing[] holes = null;
    GeometryFactory fact = new GeometryFactory();

    //create the hole in the shape if innerBoundry is > 0
    if (innerBoundry > 0) {
        gsf.setSize(innerBoundry);
        Polygon innerPoly = gsf.createCircle();
        Coordinate[] innerCoords = innerPoly.getCoordinates();
        LinearRing hole = fact.createLinearRing(innerCoords);
        holes = new LinearRing[]{hole};
    }

    //create the outter boundry
    gsf.setSize(outterBoundry);
    Polygon outterPoly = gsf.createCircle();
    Coordinate[] outterCoords = outterPoly.getCoordinates();
    LinearRing shell = fact.createLinearRing(outterCoords);

    Polygon searchPoly = fact.createPolygon(shell, holes);

    return searchPoly;
}

Here's how I am querying the QuadTree;

//create polygon to search for points the intersect with it
Polygon polygon = data_bounds.searchAreaPoly(shell, hole);

List<MyQuadNode> items = quadTree.query(polygon.getEnvelopeInternal());

However the QuadTree returns a list that doesn't remove the Points the fall within the hole, no matter how big the hole.  Does anyone know how to do this, or even if it is possible?


texts: Changing distance units in joined layer in QGIS 3.10

texts:
I am performing a Vector analysis using the 'Join attributes by nearest' tool. I am trying to evaluate the distance between a central highway and each vector point in a layer (which each represent buildings).
The Joined layer shows each of these distances but the unit for distance used is currently degrees. How can I get distances to be shown in metres or kilometres?


texts: Converting .pbf to PostrgreSQL table using osm2pgsql?

texts:
I try to import a .pbf file into a PostgreSQL table using osm2pgsql in python using the subprocess module.

The command works perfectly in the command-line, but when I try to use it in subprocess it errors out with:

  osm2pgsql version 0.92.0 (64 bit id space)

  Osm2pgsql failed due to ERROR: Couldn't open style file
  'C:/libs/share/osm2pgsql/default.style': No such file or directory

The Command Prompt is:

osm2pgsql -c -l -r pbf -p singapore-d osm -U postgres -H localhost -S D:/default.style D:/malaysia-singapore-brunei-latest.osm.pbf

The Python Code is:

import subprocess

subprocess.call(['osm2pgsql', '-c', '-l', '-r', 'pbf', '-p', 'singapore', '-d', 'osm', '-U', 'postgres', '-H',
             'localhost', '-S', 'D:/default.style','D:/malaysia-singapore-brunei-latest.osm.pbf'],
            shell=True, env={'PATH': 'C:\osm2pgsql'})


texts:
In QGIS 1.8 I have two rasters derived from the same multi-band tif image

ndvi - each pixel contains NDVI values. Using
/Plugin/Analysis Point Sampling tool
a new vector layer is created and the NDVI values of selected pixels are saved in the attribute table of this new layer. All is OK here.
geo
a georeferenced raster of the same multi-band photo

How do I get the coordinates of the same pixels from raster (2) geo?

I would like to store them into the same attribute table where the NDVI values are.


texts: How can i set Vector layer Fill as Transparent Fill in QGIS with Python

texts:
I want to fill vector layer as transparent fill. I am using following code but it fills it with white colour and inner boundaries for vector layer get hid

Code i am using :

   from PyQt4 import QtGui
   layer = iface.addVectorLayer("path_to_vector_layer", "layer_name", "ogr")
   myVectorLayer = iface.activeLayer()
   myRenderer  = myVectorLayer.rendererV2()

   if myVectorLayer.geometryType() == QGis.Polygon:

     mySymbol1 = QgsFillSymbolV2.createSimple({'color':'#ffffff',
                                      'color_border':'#000000',
                                      'width_border':'0.6'})
   myRenderer.setSymbol(mySymbol1)
   myVectorLayer.triggerRepaint()


texts: Adding OpenLayers basemap to a QGIS project changes the CRS to WGS 84 / Pseudo Mercator

texts:
I have a QGIS (2.10.1) project that contains several shapefiles, all of which are EPSG:3338 (NAD83 / Alaska Albers).  The project property CRS is also EPSG:3338.  All of the shapefiles display correctly, and the project is to my liking.

However, when I subsequently add a basemap via the Web > OpenLayers plugin, the project properties are immediately forced to EPSG:3857 (WGS 84 / Pseudo Mercator), and the map display changes accordingly.

I do not want the map to display in Pseudo Mercator.  How do I maintain the project's original NAD83 / Alaska Albers CRS while displaying an OpenLayers basemap?

Note that the project's "on the fly CRS transformation" setting has no effect.   Regardless of its initial setting, it always will be checked on after adding a basemap.


texts:
I'm using this ArcGIS map service as the basis for an application I'm developing.

However, this base map is in English, and for my application, it's necessary that the labels be displayed in Brazilian Portuguese.

Is it possible to change the language of this base map?

If this is not possible, is there an alternative to the ArcGIS map service where you can change the language?


texts: How do I get raster attribute 'labels' from a WMS GetFeatureInfo' request?

texts: I have a raster with 95 different integer values (1-95). When I do a GetFeatureInfo request, GeoServer sends back the raster cell integer value. That makes sense, but what if I want a text value (a descriptive label). I have the descriptive label in my style (.sld) and the legend is generated using this label. I'd like to get that label back through GetFeatureInfo. ArcGIS Server solves this problem by associating an attribute table to the raster (a row for each unique raster value). Is there a way to accomplish this sort of thing in GeoServer/WMS?


texts: counting values in an update statement postgres

texts:
I have a repetitive method I want to execute and I expected underneath script to work. I found out the count() function is not allowed in postgres update statements. What is the most logical and simple way to update my column with the count of occurences I encounter, without creating a new table.

UPDATE
     table1 as tbl1
SET
     new_column = count(tbl2.*)
FROM
     table2 as tbl2
WHERE
     ST_Intersects(tbl1.geom, tbl2.geopunt)


texts: Copy features from one layer to another with ogr and python

texts:
I have two layers with the same feature definitions named lyr and lyr_add.  These layers were populated from a JSON based mapserver.  I would like to append all of the features from lyr_add to lyr, but, when I use lyr.CreateFeature(feature) for a single feature, it is not appended.  I would like to append the features to the lyr prior to writing it to a file.  Running my python script does not result in any errors or exceptions.  What should I modify in my code to successfully append the features?

An abbreviated version of my code that I have attempted is:

import ogr
import os
import sys
import urllib, json

os.chdir(r'C:\Users\pipi\Documents\Rogaine')

xmin = 16501408.594
ymin = -4372525.989
xmax = 16520000
ymax = -4320000
extent = str(xmin) + ',' + str(ymin) + ',' + str(xmax) + ','+ str(ymax)

dict = {}
dict['Creek_Unnamed']=49

def count(url):
    """Returns a count of all of the features within the provided boundary"""
    url_count = '&where=1%3D1&returnCountOnly=true'
    url = url + url_count

    response = urllib.urlopen(url)
    data = json.loads(response.read())
    count = data[u'count']
    return count

def object_list(url):
    """Returns a list of the OBJECTID'S for all of the features within the provided boundary"""
    url_id = '&returnIdsOnly=true'
    url = url + url_id
    response = urllib.urlopen(url)
    data = json.loads(response.read())
    list = data['objectIds']
    return list

url_server = "http://maps.six.nsw.gov.au/arcgis/rest/services/sixmaps/LPIMap/MapServer/"
url_2 = "/query?geometry="
url_3 = "&geometryType=esriGeometryEnvelope"
url_4 = "&f=pjson"

for key in dict:
    url = url_server + str(dict[key]) + url_2 + extent + url_3 + url_4
    count = count(url)
    if count == 0:
        pass
    else:
        ds = ogr.Open(url)
        lyr = ds.GetLayer()
        if count > 1000:
            object_list = object_list(url)
            low_bound = object_list[x*1000]
            if ((x+1)*1000) > count:
                up_bound = object_list[-1]
            else:
                up_bound = object_list[((x+1)*1000-1)]

            url_loop = "&where=OBJECTID+>%3D+" + str(low_bound) + "+AND+OBJECTID+<%3D+" +str(up_bound)
            url_2 = url + url_loop
            ds_add = ogr.Open(url_2)
            lyr_add = ds_add.GetLayer()
            for i in range(0, lyr_add.GetFeatureCount()):
                feature = lyr_add.GetFeature(i)
                lyr.CreateFeature(feature)
            del ds_add

        del ds, lyr


texts: Creating vector MbTile file from source of lat / long points?

texts:
I need to create a vector mbTile, not a raster one, to be compatible with the application that uses it.

I have a postgres database that has many rows, each one representing an item with a lat and long position. Of course I can export that to CSV if required.

What should I use to create the vector mbtile? I have tried with tilemill, but it creates a raster file.


texts: Multiple overriding rules for styling in QGIS 3.4.0

texts:
I am a total beginner to QGIS and still struggling with the basics. After a quick research on this forum I could not find the answer to my current issue.
What I am trying to do is to colour the regions of a country shapefile where using the information stored in a series of categorical variables. So for each variable if that category is "yes" I want the region to be coloured in one colour (also I would later like to have a legend for the colours, but I'll worry about that later..). I have started doing this using rules-based styling in the attribute table, but I am unsure how to do it when there is conflicting rules, when a region has multiple entries of "yes". If I use single rules, then colours are overwritten by the rule hierarchy. Ideally I would like to have this region colours in stripes.
Does anybody have an idea how to do this?
Here is a description of my data structure.


texts: Landsat 8 cloud band: how to understand their (bit) coding?

texts:
I am trying to understand how to use the Lansdat 8 cloud quality band. The user guide L8_Surface-Reflectance-Code-LASRC-Product-Guide.pdf does not explain much, but just gives tables (Table 7-1, 7-2 and 7-3) I don't really understand how to read those, and how to use them.

Could someone help walk through this table? My understanding is I need to convert the results of the pixel_qa band into bits. In R, I use intToBits() and get a sequence of bits, say for pixel value 388:

intToBits(388)
#>  [1] 00 00 01 00 00 00 00 01 01 00 00 00 00 00 00 00 00 00 00 00 00 00 00
#> [24] 00 00 00 00 00 00 00 00 00

How do I interpret this? What about the Cloud confidence Attribute/ Do I have to do bite addition of columns 6 and 7? I imagine also that some attribute contradict each other (cannot have clear and cloud and high cloud confidence in the same time?)


texts:
How to delete GCP set in GEOTIFF?

I have some GEOTIFF files.
They have been set for GCP, but I made a mistake.
Therefore, I want to take the following procedure.

Delete the set GCP.
Grant new GCP.

2 is I know.

gdal_translate - gcp [col plxel] [row plxel] [lon] [lat] [input file name] [output file name]

However, I do not know the method of 1.


texts: Splitting lines once at specified distance from each end?

texts:
This is related to an old problem of mine, regarding arrowheads on lines, which doesn't look very nice unless the line is very thin. It's even an ArcGIS idea with quite a lot of votes (although no apparent attention from ESRI).

My current workaround consists of splitting the lines at a specified distance from the end, distance depending on the intended scale. Then change the symbology on the end part to contain arrows. Hopefully they will connect nicely to the rest of the line, which then can be a lot thicker.

However, this time I need arrows on more lines than I care to edit manually, and I need arrows in both ends. So I figured out I need to split my lines at about 110 m from each end. How can I accomplish that? There are no vertices, points or anything at that location which I can use. All my lines are of different lengths. Every function I have tried (native ArcGIS functions, ET Geowizards) can only split the lines multiple times in equal segments, and only count from one side. I only need one segment at each end, and it has to be 110 m.

I guess I could use the split function twice on each line, select the segments I need and be happy with that. It seems highly ineffecient though.

How should I go about fixing this, or are there any other workarounds to get my arrows look nice? Python solutions could work as well.

I have a Basic license, Python 2.7, Spatial analyst and 3D analyst to work with.

I have a number of lines that looks like this:


texts:
I'm using GeoServer and I need to put the value of NDVI from four different wms into a chart. I'm using also OpenLayer for my map.
map.on('singleclick', function(evt) {
  var viewResolution = view.getResolution();

  var ndvi_array = [], ndvi_20200706, ...

  var url_20200706 = sourceNDVI_20200706.getFeatureInfoUrl(
    evt.coordinate, viewResolution, 'EPSG:3857',
    {
      'INFO_FORMAT': 'application/json',
      'propertyName': 'GRAY_INDEX',
    });
  console.log(url_20200706);
    fetch(url_20200706)
      .then(function (response) {
        return response.text(); })
      .then(function (data) {

        var json = JSON.parse(data).features[0];

        var ndvi_20200706 = json.properties.GRAY_INDEX;

      }).catch((error) => {
        console.warn(error)
      });
  console.log("NDVI 6 luglio 2020: "+ndvi_20200706);
  ndvi_array.push(ndvi_20200706);

   var url_20190706 = sourceNDVI_20190706.getFeatureInfoUrl(
    evt.coordinate, viewResolution, 'EPSG:3857',
    {
      'INFO_FORMAT': 'application/json',
      'propertyName': 'GRAY_INDEX',
    });
  console.log(url_20190706);
    fetch(url_20190706)
      .then(function (response) {
        return response.text(); })
      .then(function (data) {

        var json = JSON.parse(data).features[0];

        var ndvi_20190706 = json.properties.GRAY_INDEX;

      }).catch((error) => {
        console.warn(error)
      });

.
.
.

console.log("Array: "+ndvi_array);

});

I'm be able to see on the map all wms.
EDIT 1
Probably my question isn't clear and I've edit the code with all my purpouse.
What I want to do is put the NDVI index of all WMS into an array because I need to use an array for a chart.
EDIT 2(after @Mike indications)
I've edited all code, now in console log I see this when I click on a point:

NDVI 6 luglio 2020: undefined
Array:


texts:
I have created this shaded relief image in arcmap (first image). Is there a way to paint part of it according to geology and keep the relief texture (like the green color on the second image)?
Creating a different layer in order to paint it on another software (e.g. photoshop) will make the relief texture faint and not very visible, so I was wondering if I can do it through arcmap.


texts: Using polygons to create mask on 2d grid with Python?

texts:
I have some polygons (Canadian provinces), read in with GeoPandas, and want to use these to create a mask to apply to gridded data on a 2-d latitude-longitude grid (read from a netcdf file using iris). An end goal would be to only have data for a given province remaining, with the rest of the data masked out. So the mask would be 1's for grid boxes within the province, and 0's or NaN's for grid boxes outside the province.

The polygons can be obtained from the shapefile here:
https://www.dropbox.com/s/o5elu01fetwnobx/CAN_adm1.shp?dl=0

The netcdf file I am using can be downloaded here:
https://www.dropbox.com/s/kxb2v2rq17m7lp7/t2m.20090815.nc?dl=0

I imagine there are two approaches here but I am struggling with both:

1) Use the polygon to create a mask on the latitude-longitude grid so that this can be applied to lots of datafiles outside of python (preferred)

2) Use the polygon to mask the data that have been read in and extract only the data inside the province of interest, to work with interactively.

My code so far:

import iris
import geopandas as gpd

#read the shapefile and extract the polygon for a single province
#(province names stored as variable 'NAME_1')
Canada=gpd.read_file('CAN_adm1.shp')
BritishColumbia=Canada[Canada['NAME_1'] == 'British Columbia']

#get the latitude-longitude grid from netcdf file
cubelist=iris.load('t2m.20090815.nc')
cube=cubelist[0]
lats=cube.coord('latitude').points
lons=cube.coord('longitude').points

#create 2d grid from lats and lons (may not be necessary?)
[lon2d,lat2d]=np.meshgrid(lons,lats)


texts: UTM to WGS84 relation between longitude and northing

texts:
Problably I'm doing something wrong or perhaps I have a misknowledge.

I tought that Latitude ~ Northing, so if I do a conversion of two position with the same latitude I will get same northing, but not.

Converting:

40.97ºN 5.7ºW --> 272811E 4538938N (30N)

40.97ºN 5.65ºW --> 277018.6E 4538809.1N (30N)

I was wrong? I'm using a wrong tool?
I know that it's no noob question, sorry.

If I do the same in lat 0 obviously I get 0 Northing so I'm lost, I don't understand the northing difference.

Can you explain it or give a link where I can understand it.


texts: How to run a spatial query with a raster and vector layer in QGIS?

texts:
I have a soil map (vector-polygon layer) and a map of the crops farmed in the same area (raster layer, from the US Cropland Data Layer).

I would like to calculate for a given soil type, the surface farmed per crop. How do I do that?

The Zonal Statistics and ZonalStats plugins didn't solve my problem.

Hope I am clear enough.
Thanks


texts: Get value, made with expression builder, with python script

texts:
I am trying to get a value created with the expression builder. The expression is in a virtual field. The expression works but when I try to get the value from my python script, it always returns None. Because my python script is a custom form I automatically get the right layer object. This is how I am trying to get the calculated value.

for f in layer.getFeatures():
        print(f["exprAge"])

Maybe I should say that the features are created with the add feature tool and after I identify a feature with the identify tool the script will run.

I figured out that, because my layer is a shapefile, I can't access the value. Only QGIS can see virtual fields, for other layers  those fields don't exist.

How can I get the value of a virtual field or can I directly "send" the value to another non-virtual field in QGIS?


texts: Looping multiple feature classes in python for ArcGIS 10

texts:
I currently have a large list of feature datasets and dbf tables in a geodatabase. I would like to join each shapefile with its matching DBF file using the "Make Query Table" tool and then "Copy Features" to produce a new shapefile with specified fields.

I am trying to automate this process so I don't have to manually select field names for each set of data (output table would use the same fields). By using the Iterate Datasets in ModelBuilder, I was hoping to do this. However, the Make Query Table will not connect to the geodatabase and still requires me to plug in a specific table.

Does anyone know I can use looping to run this process?

Are there any similar examples using python?

Example of Geodatabase:
AAA1_part1.shp
AAA1_part1.dbf
AAA2_part2.shp
AAA2_part2.shp
.....
.....


texts: Graphical modeler in QGIS 3.0.3 - Access Vector Field in Calculator

texts:
I want to create an algorithm with the graphical modeler in QGIS 3.0.3 where i define one Vector Layer and a Vector Field from this Layer as my inputvariables. Then i want to access the field-values in the field calculator. For example:
formula in fieldcalculator:

($area/1000000) * INPUT_FIELD

But this won't work. The fieldcalculator doesn't know or can't access the field values.

When i use an algorithm and define the INPUT_FIELD as an input-variable it works, but i can't access the field-values in the calculator.

How can i acess the field-values from my input Vector Field in the fieldcalculator?


texts: Should I consider any other geocoding services?

texts:
I need to do geocode lookups so that I can show members of my website other members who are in their vicinity.  My approach will be to obtain their latitude and longitude from their mobile device and then do a proximity calculation using the Haversine formula.  If they won't allow me to use their location or if they're on a computer, I'll do a geocode lookup using either the postal code or city/region/country they provided when they signed up to obtain their latitude and longitude.

I want to obtain global data, i.e. not just North America and Europe, but Asia and South and Central America if possible.  My platform is Python and Django.  I can't spend a lot of money.  I'd also like to be able to obtain some type of support (phone, email, chat).

I've done a brief assessment of the following geocode services and am wondering if there's anyone else I should consider?

Cloudmade:  Yes.  They appear to meet my requirements.
MapQuest: Yes.  They appear to meet my requirements.
Google: No.  Their terms state that I have to display a map which I don't want to do.
Yahoo:  No.  Their business doesn't depend on this service.  Who knows how long it will be around?
Microsoft:  No.  Bing may get sold.  Again, Microsoft doesn't rely on this service so I consider them too risky.
ArcGIS: No.  They meet my requirements but they cost at least $1,500 for their Desktop product.
deCarta: No.  Their service looks good but their site says they don't have geocode data for Africa, Asia, and most of South and Central America.
Cartodb: Yes.  They appear to meet my requirements.
Texas A&M Geocoding Services:  No.  Their data appears to be U.S.-centric.
Data Science Toolkit:  No.  Support seems questionable.
NetAcuity:  No.  This is overkill for my purposes and would likely be far too expensive.
MaxMind GeoIP:  No.  This requires IP addresses which aren't always accurate.
TwoFishes:  No.  This isn't a commercial product.  No support.
Boundless:  No.  Another expensive enterprise-grade product.
Pitney Bowes:  No.  Enterprise-grade

Would Cloudmade and MapQuest be my best bets?  Is there anyone else I should be looking at?


texts: Setting up service compatible with GeoPy?

texts:
GeoPy is a great tool, but the providers lack data for growing parts of the City that students are moving into, and we can't geocode them accurately because these providers can't update their data as frequently as we need.

The most up-to-date data that we need to geocode with is served up daily by the City and County of Denver in the form of address points, parcels, and street centerlines.

Is there any documentation or examples of creating a service that can talk to GeoPy that is based on the datasets I've described above?

On our end, we can update data from the City on a daily basis into a PostGIS database and publish it via GeoServer through the various formats available to an internal (or external?) service.

From there our automated GeoPy scripts can keep our student (and more) address points accurate to the City's data.

Other local options include the PostGIS Geocoder (TIGER data won't work), MMQGIS on street centerlines (can't automate, and would prefer to use parcel or point datasets), and of course we're looking for something non-proprietary as even that solution wasn't working for us.


texts: How to reproject all the vector files in a folder to a particular projection using ogr2ogr?

texts:
We have a folder with all sorts of vector files like tab,shp,dxf. Many of the files are in different projection. I would like to convert all the files to EPSG: 4326 and eventually  the non-shapefiles to shapefiles.

I've used the following command

For  %%f in ("*.shp",".tab",".dxf") do ogr2ogr -t_SRS EPSG:4326 "ESRI Shapefile" "%~dpnf.shp" "%%f"

It say's.."Failure COuldn't fetch requested layer".

And how to save the output files in a separate folder, not in the same folder itself

Edit:

I've attached the image of the error