texts:
I try to get a list of all fields - only they one which are not hidden - from QqsVectorLayer in QGIS 3 using python API. In QGIS 2 I used to do it like that:

conf = layer.editFormConfig()
for field in layer.fields():
    if conf.widgetType(field.name()) != 'Hidden':
        print(field.alias())

What's the proper way to do that in QGIS 3?


texts: Estimating biomass in mangroves in google earth engine

texts:
I am working on a project estimating biomass production and carbon storage in mangrove forests.

I looked upon examples estimating biomass production using NPP and GPP such as https://mygeoblog.com/2016/10/19/biomass-production-from-npp/ and https://mygeoblog.com/2016/10/18/measuring-carbon-storage-from-space/ .

These examples are using Modis MOD17A3H.006: Net Primary Production Yearly L4 Global 500 m and MYD17A2H.006: Gross Primary Productivity 8-Day L4 Global 500 m  data.

But these products are lacking data in mangroves sites (I am focussing on Sundarbans and some other mangrove sites in India). Of which I think the reason is that wetlands are assigned fill values https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table/mod17a2h_v006 and https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table/mod17a3h_v006 .

Can anyone guide me what dataset would be best for mangroves? Can you guide with examples and algorithm to use?


texts: get latitude & longitude as per address given for leaflet

texts:
In my project I am using leaflet to display points.
I am getting address from my database & want to show point for that particular address.
I dident found anything which will give me latitude & longitude as per address provided.
Can any one help me on this.


texts: Usage of ELSE Conditional Statement in QGIS Field Calculator

texts:
I have an issue about how to change some data in a field not affecting unrelated data.

I want to replace '(space)' between words in label_geo_ field data by this '_' sign when (if) this sign appears in another field - label_eng data.

So I tried to write this WHEN CASE (see below).

CASE
WHEN "label_eng" LIKE '%_%' THEN replace("label_geo_", ' ', '_')
ELSE ''
END

Unfortunately it changes all data in field label_geo_ and replaces '(space)' by '_' everywhere in the field.

What should I write in ELSE to exclude all those data that do not have '_' sign in label_eng field and part of label_geo_ data will stay untouched?

(Update for image)

I want to put '_' sign only there when label_eng data has it.


texts: Is it possible to replace the SQLite library in QGIS?

texts:
On Linux I was able to use the JSON1 extension by default on a layer from a SpatiaLite database connection, but on Windows (QGIS 3.0 via OSGeo4W) it isn't enabled. Is it possible to replace the SQLite/SpatiaLite library with one where JSON1 is enabled, and if yes, how?

Update: JSON1 is enabled as of osgeo4w's sqlite3 3.17.0-4.


texts: Dynamic raster layers served via MapServer

texts:
I have some processing that runs on my machine producing between five to ten rasters layer a day. I want to serve these layers using MapServer.

I initially start with a folder that has 300 raster layers all defined in the Mapfile. In the end I could end up with + - 1500 images. How can I serve these layers and incrementally add each new file generated to the Mapfile?

I notice that the more images I have the larger the Mapfile. Any reference to documentation and examples is welcome. Can I use variables to define my Mapfile layers. All these raster layer have a simple raster style.


texts: QGIS Layer Attribute Table Field Calculator Default Values

texts:
QGIS 2.18.9
In a large project, I have more than 50 layers to which I want to add a new real number field.  The default values which are displayed in Field Calculator for field type, size and precision are Integer, 10 and 0 respectively.  To change these for each individual layer, I seem to have to enter Decimal Number, 6 and 2 every time which is very time-consuming. Is there a way I can change the default values?


texts:
Tried to add "basemap.at" WMTS (Austrian Basemap grey) in GeoServer 2.12.1 as Datastore/Layer.

The preview just shows the image attached.

In QGIS and Leaflet the WMTS works well - the getCapabilities-request shows that the service CRS is epsg:3857. Tried also with setting "use defined CRS" because geoserver shows as native one epsg:4326.


texts: How Do I Snap a CSV to a SHP Using ArcPy?

texts:
I'm a veteran Python programmer very new to ArcGIS.

I'm trying to do what seems like it would be a simple operation in arcpy, but I don't think I know enough about the ArcGIS terminology to do more than fumble around.

Here's what I'd like to do:

Load a .shp file containing terrain features (lines).
Load a .csv file containing X-Y coordinates, setting the projection.
Snap the coordinates in the csv to the .shp file.
Save the new csv data to disk.

I can do this in the ArcMap UI without any trouble, but automating this with arcpy is a confusing mess to me. To load the data files do I use CreateFeatureDataset_management, CreateTable_management, CopyFeatures_management, AddXY_management, Project_management, CopyRows_management, or something else entirely?

I would very much welcome a tutorial, overview, or clue hammer. It seems that all of the ESRI documentation starts "in the middle" and I need to start at the beginning.

I don't suppose there's something like a "record what actions I perform in the GUI and export it as an arcpy script" function?


texts: Saving to LAS file using laspy, correctly

texts:
I have a text file, each line in the text file has X,Y and Z coordinates eg:

551124.5841823021 470778.1953796738 318.1167521614109

I want to read the data and export it to LAS format, by using laspy.
The code I use, inspired from :
header = laspy.header.Header()

with open(input_file, "r") as file:
     lines = list(line for line in (l.strip() for l in file) if line)
output_X = []
output_Y = []
output_Z = []
for line in lines:
    splited = line.split()
    output_X.append(float((splited[0])))
    output_Y.append(float((splited[1])))
    output_Z.append(float((splited[2])))
output_file = input_file.rsplit(".", 1)[0]+"_changed.las"
outfile = laspy.file.File(r"{}".format(output_file),
                          mode="w", header=header)

outfile.X = numpy.array(output_X)
outfile.Y = numpy.array(output_Y)
outfile.Z = numpy.array(output_Z)

outfile.close()

The values I expect to find should be of the same type as those from txt file, however, if I read the file I get:

[((551124, 470778, 318, 0, 0, 0, 0, 0, 0),)
((551124, 470779, 318, 0, 0, 0, 0, 0, 0),)
((551124, 470780, 318, 0, 0, 0, 0, 0, 0),) ...
((551438, 470707, 315, 0, 0, 0, 0, 0, 0),)
((551438, 470708, 315, 0, 0, 0, 0, 0, 0),)
((551438, 470709, 315, 0, 0, 0, 0, 0, 0),)]

I want the entries to be of type float, with decimal. I know it should be possible, since I have LAS files received from a third party, that have decimals.
Should I be setting the scale for the LAS file on export?
Also, if I open the output file in a software that reads LAS files, ArcGIS, points are not represented.
What am I doing wrong?


texts: Find closest point on a MulitLineString for each point of a MultiPoint

texts:
I'm working on a project to connect buildings to electric power plant via a network. To this end I'm beginning by trying to connect every buildings to the network.

I've importing the shp files to a database via shp2pgsql and now I've got a table of MultiPoints (buildings, with p the number of rows) and another one of MultiLineString (network, with n the number of rows).

I used ST_ClosestPoint and ST_ShortestLine on these tables but I'm facing the same problem with both of them : it returns a table of n x p rows because it computes the closest point (or the shortest line) from each point to each line of the network. My goal is to find the closest point on any line of the network for each building.

Right now the best I've got is this query :

SELECT min(ST_Length(ST_ShortestLine(network.geom, buildings.geom))) as length, id FROM network, building group by id

It returns a table of the length of the shortest line for each building but I can't find any way to get the geometry of the concern line (in order to know the location of the point on the network).

I could do a sub-query to do it point by point, order the n rows by calculing the length of each line and keep the one with the shortest length but the final goal is to use this query on files with over 200 000 buildings and a network of over 50 000 lines so I think I need to find the most efficient way.

Moreover, I did not find any way to find the "closest closest point" using the result of the query ST_ClosestPoints(network, buildings).

Does anyone think of a way of doing this efficiently ?


texts: Convert WGS-1984-Pseudo Mercator to OSGB1936

texts: I've imported UK wide coverage of ESC wind exposure data as a raster to ArcMap 10.2. The data is in WGS 1984 Pseudo Mercator which puts it way out against the UK OS standard grid. I can't get it to transform to OSGB1936 in ArcMap. In QGIS it seems to work fine and correlates with other UK vector data. Any suggestions?


texts: Copy attributes from existing shapefile to new shapefile (Python 2.7.13, GDAL)

texts:
I am trying to copy a selection of features from a shapefile to a separate shapefile, but I am stuck with trying to set the attributes of the separate shapefile to be the same as that of the original shapefile. I first thought about creating a list of attributes with the ' Name', 'Type', 'Width' and 'Precision' of each attribute and than use that list to define the attributes of the new shapefile. But I am unsure how to go about doing this?

My code so far is just constructing a list of existing attributes and creating a new shapefile:

attributes = []
for i in range(layerDef.GetFieldCount()):
    fieldName = layerDef.GetFieldDefn(i).GetName()
    fieldTypeCode = layerDef.GetFieldDefn(i).GetType()
    fieldType = layerDef.GetFieldDefn(i).GetFieldTypeName(fieldTypeCode)
    fieldWidth = layerDef.GetFieldDefn(i).GetWidth()
    GetPrecision = layerDef.GetFieldDefn(i).GetPrecision()
    attributesList = fieldName, fieldType, fieldWidth, GetPrecision
    attributes.append(attributesList)

newShapefile = driver.CreateDataSource("BuildingsFilter.shp")
srs = osr.SpatialReference()
srs.ImportFromWkt(DataSrsWkt)
newLayer = newShapefile.CreateLayer("BuildingsFilter, srs, geom_type = layer.GetLayerDefn().GetGeomType())


texts: How to access members of GeometryCollection - Python

texts:
I am using 'intersection' function from Shapely to find junctions of a road network. As i found, intersection returns GeometryCollection. In my case it could be Point, LineString or perhaps MultiLineString. I want to segregate these objects based on their geometry. I only want Point objects for further use. I found similar post here

It suggests to use as.GeometryCollection however, in my case it gives:
p = l1.intersection(l2).asGeometryCollection()
AttributeError: 'GeometryCollection' object has no attribute 'asGeometryCollection'

here is the sample code: (I expect to get points in junctions[])

with open('./lines.geojson') as f:
routes = json.load(f)

junctions =[]
for i in range(len(routes['features'])-1):
    j =i+1
    while (j < len(routes['features'])):

        l1 = LineString(routes['features'][i]['geometry']['coordinates'])
        l2 = LineString(routes['features'][j]['geometry']['coordinates'])

        p = l1.intersection(l2).asGeometryCollection()
        junctions.append(p[0])

        j = j+1


texts: Seeking tools to write smartphone application for bus map?

texts:
I noticed that few cities provide a bus map for smartphones. As a result, tourists usually resort to just taking the subway instead because it's easy to figure out, although riding a bus is a better way to discover a city.

I was wondering what tools are available to write a smartphone application (Android and iPhone) that would do this:

Show the whole bus network
Allow the user to point their current
location (or get it from GPS, if available), then point to where
they'd like to go, and display the different bus lines they should
take to get there.

Does this require extensive computer knowledge to write?


texts:
Is it possible to show the delaunay circumcircles in QGIS? I know how to make the Voronoi polygons and Delaunay triangles.

I found a web page that shows the circumcircles at:

http://www.scribd.com/doc/147332295/A-Gentle-Introduction-to-Gis

It appears that the author did this with QGIS, but he didn't explain how.

Image from Wikipedia:


texts:
How do I convert a non-closed line (think of it as a C) to a polygon filling the area in ArcMap?

I think I can convert a ring to a polygon with the Construct Polygon function on the Advanced Editing Toolbar. However I get the following error: "No polygons were created. Ensure the selected features form continuous rings and overlap the polygons."

Obviously, I need to make my C into an O. It's just a matter of adding a line segment between the two endpoints. How do I add that segment and make the line into a ring?


texts: Equidistant Cartesian projection for NYC

texts:
I'm working on a project in which I cluster locations in New York City based on their distance. My dataset contains latitude and longitude and thus I'd like to project it to a 2D Cartesian space in which distances are preserved.

I'm not really familiar with all those concepts, but I chose to use the geotools java library to project my coordinates. However I don't know to which CRS I should project in order to be able to do so, could you help me?

I've seen some cylindrical equidistant CRS but I assume the coordinates will then be in a cylindrical space. Should I then project this space on a flat surface? Is there any equidistant Cartesian CRS?
How could I find a corresponding CRS in the online databases?


texts: Calculating raster statistics for each band in QGIS?

texts:
i need to calculate raster statistics for RGBN (4 channels, multiband image) rasters in QGIS.

In Arcgis if you use "calculate statistics" in ArcCatalog, statistics are calculated for every image band.

In Qgis if you use "Raster layer statistics" you get a log file with a summary of image stats, so you cannot extract statistics for each band.

Do you know how this is achieved in Qgis?


texts:
I'm using Ubuntu 19.10 and I've just installed QGIS 3.10 using this lines in console:

#add lines
sudo nano /etc/apt/sources.list
deb https://qgis.org/ubuntu-ltr eoan main
deb-src https://qgis.org/ubuntu-ltr eoan main

#public key
sudo apt-key adv --keyserver keyserver.ubuntu.com --recv-key 51F523511C7028C3
#install
sudo apt-get update && sudo apt-get install qgis qgis-plugin-grass

Everything seems to work just fine but, I found that I can't manage GPKG files. As far as I know, they should be in the gdal library. I've installed 2.4.2 version and I cant find GPKG format within it.

> ogr2ogr --version
GDAL 2.4.2, released 2019/06/28
> ogr2ogr --formats | grep GPKG
...

> gdalinfo --formats | grep ECW
ECW -raster- (rw+): ERDAS Compressed Wavelets (SDK 5.4)
JP2ECW -raster,vector- (rw+v): ERDAS JPEG2000 (SDK 5.4)
> gdalinfo --formats | grep GPKG
....

Should I upgrade gdal? It would be enough with some extra libraries?


texts: Compare two forest loss classifications and fill missing years: accuracy assesment or correlation?

texts:
I have a time series of forest loss rasters originating from two different classifications with some missing yearly data. Rasters represents forest loss (1) or no disturbance (0) in specific year. I think that one classification underestimates the other.

Example of time series (t1, t2.,,,) and missing years from one classification:

I am curious if there is a way to compare the accuracy between overlapping years, and based on this value calculate the "underestimation" between two classifications, and forest loss values for missing years?

I think I can calculate accuracy assessment for overlapping years, but maybe some sort of raster correlation is better to predict forest loss values on missing years from class1, based on class2?

My dummy example illustrates rasters from two classifications in one year, calculates the accuracy and the "underestimation" between both.

Dummy data (one year)

library(raster)

val1<-c(0,1,1,1,0,
           1,1,1,1,1,
           1,1,1,0,0,
           0,1,1,0,0,
           0,1,1,0,1)

val2<-c(0,1,1,1,0,
        1,1,1,1,1,
        1,1,0,0,0,
        0,1,0,0,0,
        0,1,0,0,1)

r1<- raster(ncol = 5, nrow = 5)
r2<- raster(ncol = 5, nrow = 5)

values(r1)<-val1
values(r2)<-val2

par(mfrow = c(1, 2))
plot(r1, main = "r1:", legend = FALSE)
plot(r2, main = "r2:\nUnderestimate")

# Calculate error matrix
out<-table(val1, val2)

#       val2
# val1  0  1
#    0  9  0
#    1  3 13

# User accuracy
diag(out) /rowSums(out)

# Producer accuracy
diag(out)/colSums(out)

sums <- vector()

for(i in 1:dim(out)[1]){
  sums[i] <- out[i,i]

}

o.accuracy <- sum(sums)/sum(out)

What are your suggestions to predict forest loss data for missing years for classification1 based on classification2? Should I go by accuracy assessemnt on overlapping years or should I apply a correlation between rasters?


texts: Use arcpy to connect to database without connection file

texts:
It seems like arcpy can only work with connection files. This creates one of two problems:

Deployment requires manual creation of connection files.
We have to write a bunch of code to recreate the connection file on each run. (Creating it only the first time is not feasible since the connection information may change.)

As a software developer who does not have access to the production environment, I want to minimize manual steps in the deployment, and I want to keep my code simple. I also want this information to be very easy to change and not require opening up ArcMap (which would be manual), so I want to store it in a simple text config file and establish connections on the fly in my script.

Is there a way to create a connection to a database without creating a connection file? An "in memory only" connection (like is possible with pretty much any other database technology, even Oracle, or a number of other libraries, like GDAL)?

I'm using ArcGIS 10.2, and I will only be establishing query layers.


texts: How to create ArcSDE Geodatabase in PostgreSQL without ArcGIS Desktop?

texts:
I am programming a small application using Arcobjects, and I need to create an ArcSDE geodatabase on PostgreSQL.

In ArcSDE at 10.2 there is no Post Installation Wizard.

How do I create a Geodatabase in PostgreSQL using ArcSDE line Command? i.e. without needing to use ArcGIS Desktop (Create Entreprise Geodatabase toolbox).


texts: Add tiles created with Qtiles as a layer into Qgis

texts:
How do I add tiles that have been created using the Qtiles plugin as a layer into QGis and show in my project?

Ideally I would like to be able to export this as part of a QGis2Web project. I have seen the tutorial here which shows how to add the layer as a basemap after exporting through QGis2web, but I would like to be able to have in in my project as a working layer rather than a base map.


texts: Creating shapefile layer with OSGeo4W shell

texts:
I want to create an new empty shapefile with OSGeo4W shell.

Is there a command line like:

gdal_create -f "ESRI Shapefile" -wkt "point" -srs "EPSG:31256" -o output.shp


texts: Rules for scale specific layer style in GeoServer

texts: I've got a very large polygon land cover map (millions of polygons) over a large area. Our first try at serving the data through GeoServer WMS was not very successful due to the large amount of small polygons. So it was suggested to use scale specific styling (along with spatial indexes) to show only the larger polygons at a broader scale and refine as the user zooms in. This process of setting which polygons appears at each scale, based on polygon size, is a trial and error approach? Or are there rules that one could follow in order to identify which features to display depending on the scale?


texts: How to export a layout to GeoPDF format using QGIS on Ubuntu 20.04?

texts:
I'd like to export a layout in GeoPDF format using QGIS 3.12.3-București (1:3.12.3+32focal) installed on Ubuntu 20.04 Focal Fossa from https://qgis.org/ubuntu repository.
QGIS is compiled and runs against GDAL/OGR 3.0.4 (3.0.4+dfsg-1build3).
The poppler library (libpoppler97 0.86.1-0ubuntu1) is installed with the following packages:

libpoppler-cpp0v5/focal,now 0.86.1-0ubuntu1 amd64
libpoppler-dev/focal,now 0.86.1-0ubuntu1 amd64
libpoppler-glib8/focal,now 0.86.1-0ubuntu1 amd64
libpoppler-private-dev/focal,now 0.86.1-0ubuntu1 amd64
libpoppler97/focal,now 0.86.1-0ubuntu1 amd64
poppler-data/focal,focal,now 0.4.9-2 all
poppler-utils/focal,now 0.86.1-0ubuntu1 amd64

Both
gdalinfo --formats | grep -i pdf

and
ogrinfo --formats | grep -i pdf

output:
PDF -raster,vector- (rw+vs): Geospatial PDF

However, the "Create Geospatial PDF (GeoPDF)" checkbox is disabled in the PDF Export Options window with the following explanation:
"GDAL PDF driver was not built with PDF read support. A build with PDF read support is required for GeoPDF creation."
What else should I do and what other package should I install in order to enable the GeoPDF exporting feature?


texts: Same SQL expression works with search cursor but not make feature layer

texts:
I need to get features in a File GeoDatabase feature class whose "fieldName" (string field) value matches the value of "strValue" (string variable). While I know there are features satisfying this condition, Make Feature Layer seems to return an empty layer (the "Describe" object of it has an empty "FIDSet" list).  On the other hand search cursor with the same SQL expression successfully returns those features. I'm scratching my head what I am missing.

This is the where clause:

wc = arcpy.AddFieldDelimiters(FGDB_FC, fieldName) + ' = ' + '\'' + strValue + '\''

This worked as expected:

rows = arcpy.SearchCursor(FGDB_FC, wc)

But this returned a seemingly empty feature layer (dd.FID has no element):

arcpy.MakeFeatureLayer_management(FGDB_FC, "fl", wc)
dd = arcpy.Describe("fl")
print dd.FIDSet.split(';')


texts: 2D canvas library for OpenLayers. Drawing custom interactive objects

texts:
I am looking for the 2D canvas library like Konva.js that I can use above OpenLayers 5.x.
I need to draw shapes with font inside that are interactive - can be dragged over the map to the appropriate location and rotated.
On the image below there is how it looks like with Konva.js over the image background.
Can you advise proper libs or techniques so the same can be implemented over the OpenLayers map?


texts: COM Exception when setting the DisplayField property

texts:
I have a  C# .NET application referencing ArcGIS Engine 10.2.

The last line of this fragment of code causes a COM Exception:

gridLayer = new FeatureLayerClass();
gridLayer.FeatureClass = layersFeatureWorkspace.OpenFeatureClass(layerName);
gridLayer.Name = "Density Grid";
gridLayer.Cached = true;
gridLayer.Selectable = true;
gridLayer.ShowTips = true;
gridLayer.DisplayField = "NumEvents"; // Causes COM exception immediately!

The COM error code is -2147467259, which doesn't appear in this list of error ranges: http://resources.arcgis.com/en/help/arcobjects-net/conceptualhelp/index.html#//0001000002zz000000

What is going on here and how can I resolve the issue?


texts: Creating donut polygon in existing polygon using ArcGIS Desktop?

texts:
I am tracing polygons in ArcGIS with the help of editor tool but I need to exclude some areas from the existing polygons.

How can I make it?

I went through google, but instructions provided are either not clear or i am not able to understand what they say.


texts: Matching points to grid-cells - using join attributes by location

texts:
I have point data in a shapefile covering most of Africa. I've made 0.1x0.1 grids using the make grid option from QGIS and now I want to link the grid_id to my point data. Below you will find a picture of all my data points an the grid and a zoom in where the grids can be distinguished from each other - I want every point to get assigned a grid_id number base don which grid they are located in.
I've tried using join attributes by location, where my input layer was the point file and my join layer was the grids I've just made. I've experimented with other geometric predicates, but I would think I need to use the within option. However, when I run the program, he does not match any of the grid ids to my point data.
I've also tried using the MMQGIS spatial join option, but this won't allow me to select a spatial operation once I've selected the layers I want to combine.
Does anyone have any suggestions on how to link a grid_id to point data?


texts: How to change the version requested from QGIS Plugin Repository?

texts:
I built QGIS from GitHub master branch (at about version 2.14.3) and am running from the build directory but it seems to be hardcoded to connect the the plugin repository sending a version of 1.19:

You don't seem to be able to edit this parameter.  Where does this value come from and how can I change it to fetch version 2 plugins?


texts: .net arcobjects delete point feature by one of the attribute value

texts:
What I want is to delete a feature by objectID,  currently I have got the objectID of the related feature, but I have no idea how to delete the feature, can anybody give me some advice?

The following is what I have now.

Public Sub RemovePoint(ByVal lpCreatWS As IWorkspace, ByVal lpFeatureLayer As IFeatureLayer,  ByVal lsLayerName As String, ByVal lpFeatCursor As IFeatureCursor)

Dim lpCreatVer As IVersion
Dim lpVersionEdit As IVersionEdit
Dim lpWorkSpaceEdit As IWorkspaceEdit
Dim lpVerFC As IFeatureClass

lpCreatVer = lpCreatWS
lpCreatVer.RefreshVersion()

lpVersionEdit = lpCreatWS
lpWorkSpaceEdit = lpCreatWS

Try

      lpFeature = lpFeatCursor.NextFeature

      lpWorkSpaceEdit.StartEditing(True)

      Do Until lpFeature Is Nothing

            lpVerFC = lpObjSDEUtils.GetFeatureClass(lpCreatWS, lsLayerName)
            lpVerFeat = lpVerFC.GetFeature(lpFeature.OID)

Catch ex As Exception

End Try


texts: Using Google Earth Web as an application development platform

texts:
I'm looking for an option to use Google Earth Web as a platform to develop a custom application, which would involve reviewing/editing a large list of map locations and connections between them.
The 3D content of GEW would be great for this. However, I could not find any documentation of using it for this purpose (similar to how we could use Leaflet or OpenLayers for 2D content).  Is using GEW for building custom apps possible?


texts: How do I use gdal command-line utility to extract raster attributes as input for gdalwarp?

texts:
I need to reproject a raster into the same extent and resolution as an existing raster, and want to use the gdal command-line utilities to do it.

How do I extract these values in a format I can use as input for gdalwarp?

Here's what I've tried:

function gdal_extent_gdalwarp_te() {
    if [ -z "$1" ]; then
        echo "Missing arguments. Syntax:"
        echo "  gdal_extent <input_raster>"
        return
    fi
    EXTENT=$(gdalinfo "$1" |\
        grep "Lower Left\|Upper Right" |\
        sed "s/Lower Left  //g;s/Upper Right //g;s/).*//g" |\
        #tr "\n" " " |\
        sed 's/ *$//g' |\
        tr -d "\n[(,]'")
    echo -n "$EXTENT"
}

function gdal_pixelsize_gdalwarp_tr() {
    if [ -z "$1" ]; then
        echo "Missing arguments. Syntax:"
        echo "  gdal_pixelsize_gdalwarp_tr <input_raster>"
        return
    fi
    EXTENT=$(gdalinfo "$1" |\
        grep "Pixel Size" |\
        sed "s/Pixel Size =//g; s/,/ /g" |\
        tr "\n" " " |\
        tr -d "[(,])-")
    echo -n "$EXTENT"
}

ExtractedSRS=$(gdalsrsinfo -o proj4 "SRSraster.tif" | tr -d "'")
echo "d"
ExtractedExtent=$(gdal_extent_gdalwarp_te "SRSraster.tif")
echo "e"
ExtractedResolution=$(gdal_pixelsize_gdalwarp_tr "SRSraster.tif")

#in this case I'm trying to tile the input rasters into a new output raster
gdalwarp -t_srs "$ExtractedSRS" -te "$ExtractedExtent" -tr "$ExtractedResolution" inputrasters*.tif outputraster.tif

This gives the error:

Creating output file that is -124997P x 0L.
ERROR 1: Attempt to create -124997x0 dataset is illegal,sizes must be larger than zero.

If I try leaving out the -te option, like this,

gdalwarp -t_srs "$ExtractedSRS" -tr "$ExtractedResolution" inputrasters*.tif outputraster.tif

I get:

Usage: gdalwarp [--help-general] [--formats]
    [-s_srs srs_def] [-t_srs srs_def] [-to "NAME=VALUE"]
    [-order n | -tps | -rpc | -geoloc] [-et err_threshold]
    [-refine_gcps tolerance [minimum_gcps]]
    [-te xmin ymin xmax ymax] [-tr xres yres] [-tap] [-ts width height]
    [-wo "NAME=VALUE"] [-ot Byte/Int16/...] [-wt Byte/Int16]
    [-srcnodata "value [value...]"] [-dstnodata "value [value...]"] -dstalpha
    [-r resampling_method] [-wm memory_in_mb] [-multi] [-q]
    [-cutline datasource] [-cl layer] [-cwhere expression]
    [-csql statement] [-cblend dist_in_pixels] [-crop_to_cutline]
    [-of format] [-co "NAME=VALUE"]* [-overwrite]
    [-nomd] [-cvmd meta_conflict_value]
    srcfile* dstfile

Available resampling methods:
    near (default), bilinear, cubic, cubicspline, lanczos, average, mode.

FAILURE: Wrong value for -tr parameters.


texts:
There are two shapefiles. One is a compilation of municipalities on the coast with borders extending into ocean. The second shapefile is a coastal border.

I'm trying to clip the municipality shapefile with the coastal border shapefile to remove the municipality sea territory from my image, but I've had little success. The municipality boundaries are contiguous, so the coastal border on the second shapefile needs to cut and remove the excess shape.

I'm using QGIS 2.14.3-Essen. I've googled around and tried several strategies. I've used a lot of the GeoProcessing Tools (difference, intersect, clip, union) and what returned was an invisible shapefile.

What should I try to make it work?

Below is a copy of the screen settings for the clip. The results was yet another invisible layer.

Below this is a screencap of municipality boundaries to be cut (CITY) and the coastline in red to cut them (SMC_BORDER).


texts:
I am working on a project which runs on ArcGIS platform using VB.NET. It is to model a sub-surface fault plane (currently a plane rectangular model with some dip) and then to find required distances from various points on surface of the earth. These distances are not necessarily nearest distances.

One idea I got and decided to use is TIN network.

I want to know if this is productive given that the number of points I know on surface would be only 5.

Is there any better option to model this?

And one more point I want to remind is that we can not model perfectly vertical faults using that.


texts:
I am new in PostGIS and want to learn. I faced the query where I have to find one point in two different layers or table. I created a query like this but it is not giving the desired result:

SELECT b.name as "Boroughs_Name", r.name as "Restaurant_Name", l.lsoa11nm as "Lsoa_Name"
FROM italian_restaurants as r, .london_boroughs as b, .london_lsoa as l
WHERE r.name= 'Pasta House' and st_intersects(b.geom, r.geom) or st_intersects(b.geom, l.geom)


texts:
I'm working to build a dynamic map, which pulls data from GeoServer to get the list of available layers, and then uses Leaflet to allow the user to determine which specific layers to display.  Leaflet and OpenLayers both support the context of a base layer, in which the user selects one base, with as many overlays as they then need.

I'm successfully pulling the list of layers from GeoServer's GetCapabilities call.  However, I've not determined a way to suggest which layers are good candidates for base, versus which should be considered overlays.  I've cross-checked the responses returned for WMS and WFS, as well as version 1.1 versus 1.3.  Suspect it's a concept issue on my part.  Can someone clarify how to distinguish between base and overlay layers for layer lists returned from GeoServer?


texts: Fixes to achieve topology in QGis not working

texts:
I imported a .shp file to QGIS with polygons from a non-topological database. Polygons were adjacent and shared borders, but most nodes didn't coincide. I then proceeded to fix all nodes and intersections manually (more than 3000) enabling snapping (with a 10 pixel tolerance) and topological editing. Sometimes topology worked and moving the node moved both borders, most of the time it did not. I finished correcting and checking all the nodes, topology checker indicates that there are no invalid geometries but a whole lot of gaps and intersections... except they DO NOT exist (I zoomed all the way to 50:1 and still can't see them). And all nodes seem to be corectly snapped to one another. Tried using v.clean in all possible ways but nothing happens because I suppose there is no cleaning needed. My guess is that since the shapefile is from a non-topological origin QGIS is not able to create topology for all the nodes even if they have been snapped together.

Does anyone have any idea how to fix this?

Or are there still elements missing from QGIS to achieve this?


texts: Restart ArcGIS Server Map Service From Command Line

texts: Is there a way to restart an ArcGIS Server v10 Map Service from the command line in Windows 7?


texts: Subtract MultiPolygon table from LineString table

texts:
I've two tables, one with a roadnetwork and another one with some polygon. The roadnetwork table has a geometry column defined as

COLUMN geom geometry(LineString,4326);

while the polygon table has a geometry column defined as

COLUMN geom geometry(MultiPolygon,4326);

I want to perform a Boolean operation, by subtracting the polygon geometry from the roadnetwork geometry, in order to make an hole inside my roadnetwork.

I've tried the following query, that I've found in stackoverflow:

UPDATE public.roadnetwork
SET geom = ST_Multi(ST_Difference(public.roadnetwork.geom, public.sub_polygon.geom))
FROM public.sub_polygon
WHERE ST_Intersects(public.roadnetwork.geom, public.sub_polygon.geom);

But I got the following error:

subtractpolygon.sql:4: ERROR:  Geometry type (MultiLineString) does not match column type (LineString)

How can I subtract the polygons defined in polygon table from the roadnetwork table?

EDIT:

If I run this query:

SELECT ST_AsText(ST_Difference(public.roadnetwork.geom, public.sub_polygon.geom))
FROM public.roadnetwork, public.sub_polygon
WHERE ST_Intersects(public.roadnetwork.geom, public.sub_polygon.geom)

I obtain some results, in which there are both LINESTRING and GEOMETRYCOLLECTION EMPTY:

 GEOMETRYCOLLECTION EMPTY
 GEOMETRYCOLLECTION EMPTY
 LINESTRING(51.5118811558203 25.3205662647387,51.512059 25.3202912,51.512134 25.3201413,51.5121684 25.320028,51.5121684 25.319968,51.5121396 25.3198288)
 GEOMETRYCOLLECTION EMPTY
 LINESTRING(51.5128903 25.3196673,51.5129743 25.3198392,51.5130642 25.3199292,51.5132841 25.3200641,51.5136939 25.320329,51.5139688 25.3205439,51.5143176 25.3208249,51.5145508 25.321047,51.5146999055034 25.3211917269465)
 GEOMETRYCOLLECTION EMPTY
 GEOMETRYCOLLECTION EMPTY
 GEOMETRYCOLLECTION EMPTY
 LINESTRING(51.5170740202354 25.331309846686,51.5170939 25.3329918,51.5171189 25.3333467,51.5171813 25.3336378,51.5172913 25.3339226,51.5174162 25.3341775,51.5175362 25.3343724,51.5176661 25.3345324,51.5178627 25.3346364)

I think that when I update the geom column I got the error because I've empty GEOMETRYCOLLECTION in query results, so it cannot be inserted in geom table that wants LineString. How can I change the query in order to filter unwanted data?

EDIT2

I've updated my query in order to remove empty geometry collections that don't affect my results:

WITH RESULTS as (
SELECT ST_Difference(public.roadnetwork.geom, public.sub_polygon.geom)
FROM public.roadnetwork, public.sub_polygon
),

 FILTERED_RESULTS as (
  SELECT * FROM RESULTS WHERE RESULTS.st_difference <> '0107000020E610000000000000'
)

 UPDATE public.roadnetwork
 SET geom = FILTERED_RESULTS.st_difference FROM FILTERED_RESULTS;

The problem is that now I've the following error

Geometry type (MultiLineString) does not match column type (LineString)

In this case I can't filter MultiLineString because I need these data. Probably they are results of ST_Difference because I've not them in my roadnetrowk table. So I need to find a way to process FILTERED_RESULTS in my query in order to replace MultiLinestring with corresponding LineString. I've tried ST_LineMerge but with no change in my data.

Problem in graphic way

I've this road network and this polygon:

I want to obtain this:


texts: Can't open QGIS 3.8 plugin manager on OSX

texts:
I have not been able to open the plugin manager in QGIS 3.8.  I am running Sierra OSX.  This is the message that I get when I try to open the plugin manager:

  An error occurred during execution of following code:
  pyplugin_installer.instance().showPluginManagerWhenReady()

  Traceback (most recent call last):
    File "", line 1, in
  NameError: name 'pyplugin_installer' is not defined

  Python version:
  3.6.8 (default, Jun 25 2019, 16:00:16)
  [GCC 4.2.1 Compatible Apple LLVM 9.0.0 (clang-900.0.39.2)]

  QGIS version:
  3.8.0-Zanzibar 'Zanzibar', exported

Also (this might be a related problem), when I open QGIS, I get the following dialogs:

  Couldn't load plugin 'DigitizingTools' due to an error when calling its classFactory() method
  ModuleNotFoundError: No module named 'math'

  Couldn't load plugin 'ImageServerConnector' due to an error when calling its classFactory() method
  ModuleNotFoundError: No module named 'math'

  Couldn't load plugin 'latlontools' due to an error when calling its classFactory() method
  ModuleNotFoundError: No module named 'math'

  Couldn't load plugin 'qgis2web' due to an error when calling its classFactory() method
  ModuleNotFoundError: No module named 'math'
  Couldn't load plugin 'quick_map_services' due to an error when calling its classFactory() method
  ModuleNotFoundError: No module named 'math'

  Couldn't load plugin 'MetaSearch' due to an error when calling its classFactory() method
  ModuleNotFoundError: No module named 'math'

  Couldn't load plugin 'processing'

  Couldn't load plugin 'processing'

  An error occurred during execution of following code:
  pyplugin_installer.initPluginInstaller()

I am running Python 3.7.3 using PyEnv


texts: Connecting to Microsoft SQL Server from QGIS

texts:
I can't connect my QGIS to SQL Server.

This is my database:

This is my connection configuration to QGIS

If I press test connection it says "Driver not loaded". If i load database, there is nothing happening.

What is wrong here?


texts: Creating a Heatmap of Delivery Probability in QGIS3

texts:
Im using QGIS 3.4.2 to create a heatmap to plot the source place where a network message is generated, in a vehicular network.
I separated the messages in two groups, one in each layer: delivered and not delivered.
So far is working 100%.

Now I want to combine both layers in order to generate a single heatmap, which must show the % rate of delivery combined.
For instance, if in a determined region 4 messages were generated, and only one of them could be delivered, the % rate should indicate 25%...

Currently I'm using two separated and very basic CSV files, one for delivered messages and another for not_delivered messages:

Longitude,Latitude
-5483420,-2938848
-5490938,-2938382
-5480800,-2953976

But I can easily merge them and create a single CSV similar to this:

Longitude,Latitude,Delivered
-5483420,-2938848,Y
-5490938,-2938382,N
-5480800,-2953976,N

I cannot find a single example online.


texts: Date column changes when I upload a shapefile with this column type to Google Earth Engine

texts:
I have a FeatureCollection in my Assets in Google Earth Engine, which is the result of uploading some shapefile containing some productive areas from the company I work for.
When I import this collection to my script and print it, I can see the columns which contain date values in the previous shapefile are now changed to type Long in the FeatureCollection.
How can I import this FeatureCollection without having this problem?


texts: Finding closest segment to a point using PostGIS

texts:
I am trying to find the closest road segment to a point that I have and return an attribute from the closest segment to the given points.

I know in PostGIS I need to use ST_ClosestPoint, St_distance and min functions, but I keep having syntax errors.

So far my code looks like

WITH closest_segment As
)
 SELECT
  ST_ClosestPoint(ST_GeogFromText('POINT(-122.206111 47.983056)')::GEOMETRY,
    shape::GEOMETRY)AS segment_location  --shape is the geography WKB

 FROM
  public.road_segments
)
 RETURNING closest_segment
WITH new_segment AS
)
 SELECT
  MiN(ST_Distance(closest_segment::GEOM,ST_GeogFromText('POINT(-122.206111
   47.983056)')::GEOMETRY)
 FROM
  closest_segment
)
 RETURNING closest_segment

My knowledge of PostGIS and coding in it is very limited.


texts: Projections for solar eclipse umbral paths

texts:
Most projections of the path of totality for the 2017 solar eclipse look like

 (from http://xjubier.free.fr/en/site_pages/solar_eclipses/TSE_2017_GoogleMapFull.html).

But there are also maps online using different projections such that the path is a straight line or mostly straight line, as in

 (from http://www.eclipse-maps.com)

I cannot find any CRS for the latter.  The ideal CRS for me would produce a straight path (just has to be mostly straight for the part of the path across the US) and in addition the path would be straight up and down with "west" side on the top.  I am using QGIS but I am amenable to switching to another open source package if necessary.

@MichaelZeiler said (and I have confirmed) that for this particular eclipse and this part of the run (the US) a standard projection happens to work well enough to produce a satisfyingly straight enough umbral path: EPSG:102003 = USA_Contiguous_Albers_Equal_Area_Conic

+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs

Now I am trying to figure out how to transform that into a CRS that produces maps that are rotated 67.560 degrees clockwise relative to the above.

In other words I want to change the CRS so that I get

Notice how all the labels are horizontal.  This is an advantage of changing the CRS as opposed to the Canvas and other solutions - labels just work normally for all new layers.  I still don't understand how to change the CRS after reading a number of pages.

Xavier Jubier said that a good approximation of a straight line can also be obtained with EPSG:102009 = North_America_Lambert_Conformal_Conic

+proj=lcc +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs

I have learned after many searches that there are obscure projections that will give more exact and general astronomy related solutions and anyone wanting to take this more difficult approach will want to know the keyword quincuncial.


texts: Change default field properties for number format in ArcMap?

texts: Anyone who has worked with numbers that have more than 14 digits has run into ArcMap's wonderful habit of converting any number of such length to e-notation. While changing this is a simple matter of going into the Field Properties and setting the number format, doing this repeatedly for every field in every table grows tiresome. Is there some secret environment setting or method to make number formats default to something besides "Numeric"?


texts:
I have several layers from OpenStreetMap that I wanted to convert to WGS 84 4326 due to my using several base layers that are assigned to that CRS.

No matter what I try, whenever I check the layer they all say "this layer does not have a specified CRS"

Things I have tried:

Set project to WGS 84 4326.Turned OTF on and off.
Set CRS in layer options and saved it and reloaded it.
QGIS settings to prompt me when I add a layer missing a CRS which it does not prompt.

I have used a raster image to compare against my OSM layers and it would change accordingly to my projection change. Although it would never match the location of these OSM layers regardless of projection.


texts: Subset moving window when using focal on raster in R

texts:
I have a classified Landsat image, and I'd like to pass a filter over it converting certain cells to 0.  For some categories, I'd like to set cells to 0 unless an entire 7x7 window contains only that cell value.  For other categories, I'd like to set cells to 0 unless an entire 5x5 window contains only that cell value. And for others, set to 0 unless at 3x3 window, etc.  I've tried two ways to do this:

One is to pass multiple filters over the same raster, and then combine their outputs.  I tried this and it takes an incredibly long time, so I think it would be more efficient to pass one 7x7 filter over the raster, and have a conditional function that subsets the window to 5x5 or 3x3 based on the value of the category it is passing over.

However, it doesn't seem to be working, I think because I am having trouble referencing the 3x3 sub window within the 7x7 moving window.  How can I write a function to be used within the focal() function that references only certain cells within the moving window?

The code I am using:

r <- raster('myraster.tif')

FocalSelect <- function(cells){
  if (cells[c(25)]==4){
    return(4)}
  if (cells[c(25)]==9){
    return(9)}
  if (cells[c(17,18,19,24,25,26,31,32,33)]==rep(2,9)){
    return(2)}
  if (cells[c(17,18,19,24,25,26,31,32,33)]==rep(5,9)){
    return(5)}
  if (cells[c(17,18,19,24,25,26,31,32,33)]==rep(8,9)){
    return(8)}
  if (cells[c(seq(9,13),seq(16,20),seq(23,27),seq(30,34),seq(37,41))]==rep(1,25)){
    return(1)}
  if (cells[c(seq(9,13),seq(16,20),seq(23,27),seq(30,34),seq(37,41))]==rep(6,25)){
    return(6)}
  if (cells[c(seq(9,13),seq(16,20),seq(23,27),seq(30,34),seq(37,41))]==rep(7,25)){
    return(7)}
  if (cells==rep(3,49)){
    return(3)}
  if (cells==rep(10,49)){
    return(10)}
  if (cells==rep(11,49)){
    return(11)}
  else{
    return(0)
  }
}

out <- focal(r, matrix(rep(1,49),nrow=7),fun=FocalSelect)

and this give me these warnings:

1: In if (cells[c(17, 18, 19, 24, 25, 26, 31, 32, 33)] ==  ... :
  the condition has length > 1 and only the first element will be used
2: In if (cells[c(17, 18, 19, 24, 25, 26, 31, 32, 33)] ==  ... :
  the condition has length > 1 and only the first element will be used
3: In if (cells[c(17, 18, 19, 24, 25, 26, 31, 32, 33)] ==  ... :
  the condition has length > 1 and only the first element will be used
4: In if (cells[c(seq(9, 13), seq(16, 20), seq(23, 27), seq(30,  ... :
  the condition has length > 1 and only the first element will be used
5: In if (cells[c(seq(9, 13), seq(16, 20), seq(23, 27), seq(30,  ... :
  the condition has length > 1 and only the first element will be used
6: In if (cells[c(seq(9, 13), seq(16, 20), seq(23, 27), seq(30,  ... :
  the condition has length > 1 and only the first element will be used
7: In if (cells == rep(3, 49)) { ... :
  the condition has length > 1 and only the first element will be used
8: In if (cells[c(17, 18, 19, 24, 25, 26, 31, 32, 33)] ==  ... :
  the condition has length > 1 and only the first element will be used
9: In if (cells[c(17, 18, 19, 24, 25, 26, 31, 32, 33)] ==  ... :
  the condition has length > 1 and only the first element will be used
10: In if (cells[c(17, 18, 19, 24, 25, 26, 31, 32, 33)] ==  ... :
  the condition has length > 1 and only the first element will be used

EDIT: It looks like I was subsetting the moving window correctly, but my if statement was phrased wrong.

if (cells[c(17,18,19,24,25,26,31,32,33)]==rep(2,9)){

should be

if (identical(cells[c(17,18,19,24,25,26,31,32,33)],rep(2,9))){


texts:
Sorry for my English

I have leaflet.js (image of map 5913x7863) and there L.CRS.Simple

I had

var imageBounds = [[1468.9700,-11.7050], [1968.2705,363.7705]];
var file0 = new L.imageOverlay(imageUrl, imageBounds).addTo(map);

and everything was good, but I understood, that L.imageOverlay - it's wrong way and I need tiles. So now I have

L.tileLayer('images/tiles/{z}/{x}/{y}.png', {
continuousWorld: true,
noWrap: true,
maxZoom: 6,
minZoom: 3
}).addTo(map);

But I don't know how to do L.tileLayer with this bounds (coordinates) [1468.9700,-11.7050], [1968.2705,363.7705]

Now bounds near this: [0,0], [-123,92], but I need [1468.9700,-11.7050], [1968.2705,363.7705]
And I don't understand why is so hard to do bounds with tiles as easy with imageOverlay?

I have tried this:

var southWest = L.latLng(1468.9700,-11.7050),
northEast = L.latLng(1968.2705,363.7705),
mybounds = L.latLngBounds(southWest, northEast);
...
bounds: mybounds

but it have not helped.

Update:

I trying to resolve my problem, but it don't help:

var w = 5913;
var h = 7863;
var mapMinZoom = 0;
var mapMaxZoom = 5;
var _map = L.map('map', {
    maxZoom: mapMaxZoom,
    minZoom: mapMinZoom,
    crs: L.CRS.Simple,
    zoomControl: false,
    attributionControl: false,
  });

var _mapBounds = new L.LatLngBounds(
_map.unproject([0, h], mapMaxZoom),
_map.unproject([w, 0], mapMaxZoom));
_map.setMaxBounds(_mapBounds);

var _mapCenter = _map.unproject([w/2, h/2], mapMaxZoom);
_map.setView(_mapCenter, 3);

var _tileLayer = L.tileLayer('images/tiles/{z}/{x}/{y}.png', {
    minZoom: mapMinZoom, maxZoom: mapMaxZoom,
    bounds: _mapBounds,
    continuousWorld: true,
    noWrap:true,
    tileSize:256,
    crs: L.CRS.Simple,
    detectRetina:false
  }).addTo(_map);

Now coordinates is near -247,184
and unproject coordinates is near 1480,1970

If you have a solution for L.CRS.Simple and XY coordinates - please help!


texts: Creating bookmarked PDF from Mapbook (Python) export script

texts:
I regularly use an export script to create multi-page PDFs from Data-Driven-Pages based Map-books.  I am able to export and append each individual page to the PDF, and save the PDF in the single-page layout with the "Page Thumbnails" pane opened.  I need to be able to (Pythonically) set the Bookmark text for each page within the PDF.  I can do this manually from Adobe, but it would save a lot of time and headache if it could be automated.

Here is the script I use, at its most basic level without any bells or whistles.  At the bottom of the script is where I set the page-layout (pdf_open_view="USE_THUMBS").  I can use "USE_BOOKMARKS" but this doesn't actually apply any bookmark text to the individual pages.

import arcpy, os
mxdPath = r"C:\temp\MyMapBook.mxd"
tempMap = arcpy.mapping.MapDocument(mxdPath)
tempDDP = tempMap.dataDrivenPages

outDir      = r"C:\temp"   #PDF Folder where the individual pages and Combo are created
comboPDF    = r"\MultiPage_PDF.pdf"        #Name of Combo PDF
finalpdf_filename = os.path.join(outDir + comboPDF)
finalPdf = arcpy.mapping.PDFDocumentCreate(finalpdf_filename)

for pgIndex in range(1, tempDDP.pageCount + 1):
    tempDDP.currentPageID = pgIndex
    figName = str(tempDDP.pageRow.sdsFeatureName)
    figNumber = str(tempDDP.pageRow.Fig_Num_2015FSP)
    pageName = "Page-" + figNumber + "_" + figName + ".pdf"
    individualPagePDF = os.path.join(outDir + "\\" + pageName)

    #do something to the page before exporting. e.g. add inset map, add table etc.

    ###################################################################
    tempDDP.exportToPDF(individualPagePDF, "CURRENT") # export to signle page PDF
    finalPdf.appendPages(individualPagePDF) # append the page to the "combo" PDF
    print pageName + " added to final mapbook pdf"
    ###################################################################

    #undo whatever page modifications you just did before the export

##### Delete Objects
del tempMap, tempDDP

finalPdf.updateDocProperties(pdf_open_view="USE_THUMBS",pdf_layout="SINGLE_PAGE")
finalPdf.saveAndClose()


texts: Duplicate labels in ArcGIS Server REST map service OpenLayers 5

texts:
I'm trying to export PLSS grid from ArcGIS map server and I have a duplicate labels in PLSS map service. Using OpenLayers 5

My code
    mapPLSSTownshipLayer = new ol.layer.Tile({
    extent: [-13884991, 2870341, -7455066, 6338219],
    source: new ol.source.TileArcGISRest({
        url: 'https://gis.blm.gov/arcgis/rest/services/Cadastral/BLM_Natl_PLSS_CadNSDI/MapServer/export?dynamicLayers=[{"id":1,"source":{"type":"mapLayer","mapLayerId":1},"minScale":4000000.000000047,"maxScale":0,"drawingInfo":{"renderer":{"type":"simple","symbol":{"type":"esriSFS","color":[0,255,0,0],"outline":{"type":"esriSLS","color":[38,38,38,255],"width":2}}},"labelingInfo":[{"labelPlacement":"esriServerPolygonPlacementAlwaysHorizontal","labelExpression":"[TWNSHPLAB]","minScale":550000,"maxScale":0,"symbol":{"type":"esriTS","color":[0,0,0,255],"haloColor":[255,255,255,255],"haloSize":2,"font":{"family":"Arial","size":13}}}]}}]&TRANSPARENT=true&SIZE=700%2C700',
        params: {
            'DPI':210,
            'FORMAT':'PNG24'
        },
        ratio: 1
    }),
    maxResolution: 300
});

I tried change image size but my params seems overwriting by default

Any suggestions?


texts: Restrict the scale of the search widget in ArcGIS Flex API

texts:
I wonder if I could restrict the scale scope only for the searching widget (and not the extent of the map)? We want to narrow down the scale for searching in order to avoid too many results which may cause intranet traffic jam.

I am using Adobe Flash Builder and ArcGIS Server 10.2.2.


texts: How do I import AutoCAD [DWG] drawings to QGIS [SHP] using the right CRS settings?

texts:
We did some topographical surveys in Indonesia which are drawn in AutoCAD. It's our aim to publish them with QGIS and its webclient.

I thought the best way is to georeference the data in AutoCAD 2014. Eventually, I opened the DWG-file in AutoCAD Map 3D 2014 (see AutoCAD-Map-3D-Screenshot) because I know of it's possibility to export as SHP-file. Finally, I imported the SHP-file in QGIS (nightly built) but the result is not the expected one, see QGIS-screenshot :

The white dot(s) represents GPS-data which was correctly imported and rendered in the southern part of Nias. The blue dot(s) represents the imported SHP-files which should actually be located at the same location but are -- as you can see -- in the middle of nowhere.

Do you have any clue about a mistake I made?

While exporting/ importing the data I payed attention to the CRS and set it to "DGN95.UTM-47N" in AutoCAD Map 3D and to "EPSG:23867 - DGN95 / UTM zone 47N" in QGIS -- which, in my opinion, should be the same one and lead to a right result. Unfortunately, it doesn't. Or is there something I missed?

I already thought about the "Data source encoding" which I can set in QGIS. Right now it's set to "System". I guess Debian Wheezy (which I use for QGIS) and Windows 7 (which I use for AutoCAD [Map 3D]) use a different one -- but I'd also suppose that I would not see anything at all if I import data with a wrong encoding. Some ideas?

Any help is very appreciated. Thank you.


texts: I am looking to create a raster filter which reclassifies single or contiguous pixels less than, for example, 30 pixels (shown in white below) into another class (shown in black). The filter, majority filter, aggregate and other spatial analyst generalization tools seem to be too over reaching in this case. Any Python or Arc solutions and ideas would be appreciated.


texts: Finding Two locations with conditional search on DEM?

texts:
I want to find two locations H difference in elevation between two locations on a DEM.

Ex: point A and B are my output after search on DEM.  conditions are:

buffer R distance around point A and
difference of H elevation to point A

I am very poor in programing and looking for a tool or program/script to achieve this.

No...Q is storage capacity of reservoir to generate min power as power capacity is propotional to Qh (Q storage capacity, h head availability).

Yes..I did my analysis as you mentioned above procedure by using archydro.
But it can used if you know the location and general power project which has one reservoir.

My case require two reservoirs with difference of elevation of two outlets H.
I want to find two locations at a time by search(condition of H) on DEM.


texts: How to create dashboard combining maps and charts open source way?

texts:
As far as I am aware there is no open source user interface (framework) which would allow easily design website as combination of interactive maps and charts like ArcGIS dashboard tool does. And you can't go without at least basic coding skills and combination of tools. What is currently the closest it gets to create interactive web GIS dashboard in open source world?

In particular, what would it take to create open source alternative dashboard to coronavirus dashboard? We talk client side only solution, which must be capable to connect to GIS servers (consume GIS services).

My understanding is currently you must have some coding skills, and use some combination of following tools:

Bootstrap - main framework to design page / grid
Leaflet or OpenLayers to incorporate interactive maps
D3.js or Chart.js to incorporate interactive charts

By initial search combination of Bootstrap, Leaflet and D3.js seems to be the most used way these days to incorporate interactive maps and charts into dashboard like website. Is that correct?


texts: Moving Feature Class to Different Feature Dataset within Same File Geodatabase using ArcPy?

texts:
I want to progammatically move a feature class into a different feature dataset within the same file geodatabase. I'm wondering if there is a way of doing this without including code to rename and delete data.

My current code looks something like this:

for fc in fc_list:
        out_fc = fc+'_1'
        arcpy.FeatureClassToFeatureClass_conversion(fc, archive, out_fc)
        arcpy.Delete_management(fc)
        arcpy.Rename_management(out_fc, fc)


texts:
I'm trying to georeference the .grd files provided here http://open.canada.ca/data/en/dataset/3d4b68a5-13bc-48bb-ad10-801128aa6604.

The documentation states "The CANGRD grid is in polar stereographic projection with a 50 km spatial resolution. The grid is a 125 (columns) by 95 (rows) matrix, where the SW corner (0,0) is at 40.0451°N latitude and 129.8530°W longitude. The projection is true at 60.0°N and centered on 110.0°W."

The header of each .grd file contains this:

DSAA
125 95
0.0 124
0.0 94
-3.051277 3.268792

I tried with the following Proj4 string, but the whole dataset appears tilted.

+proj=stere +lat_0=40.0451 +lat_ts=60 +lon_0=-129.8530 +x_0=0 +y_0=0 +datum=WGS84 +k=1 +to_meter=50000

Any idea how to properly georeference these files?


texts: Restrict Map Services available to a custom ArcGIS Server Print Service

texts:
I am setting up a custom print service in ArcGIS Server to use our own templates, but as this print service will run as the ArcGIS Server account it will potentially have access to all secured services.

How can I restrict/white list the map services an ArcGIS Server Print Service will have access to?

I do not wish to restrict access to the print service, it is to be publicly available, however I only want it to be available to print/export other public map services.

I can't see anything in the documentation, however could possibly setup some form of request filtering in IIS.


texts: Kartograph - AttributeError: _Shape instance has no attribute 'parts'

texts:
I am new to GIS and open source tools so bear with me.

I am attempting to convert a shapefile to svg by running kartograph as a python module.  My source file (provided by the local GIS Administration) contains nearly 50,000 parcels of land and I only want to carry over a single attribute for each parcel (which I will use for colorizing later). Here is my code:

from kartograph import Kartograph
K = Kartograph()

config = {
    "layers": {
        "county": {
            "src": "parcels.shp",
            "attributes": {
                    "id": "XXDSDS"
            }
        }
    },
}

K.generate(config, outfile='test-file.svg')

The following error results:

  Traceback (most recent call last):
    File "test.py", line 19, in
      K.generate(config, outfile='test-file.svg')
    File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/kartograph.py-0.6.8-py2.7.egg/kartograph/kartograph.py", line 46, in generate
      _map = Map(opts, self.layerCache, format=format)
    File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/kartograph.py-0.6.8-py2.7.egg/kartograph/map.py", line 48, in init
      me.proj = me._init_projection()
    File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/kartograph.py-0.6.8-py2.7.egg/kartograph/map.py", line 88, in _init_projection
      map_center = self.__get_map_center()
    File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/kartograph.py-0.6.8-py2.7.egg/kartograph/map.py", line 140, in __get_map_center
      features = self._get_bounding_geometry()
    File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/kartograph.py-0.6.8-py2.7.egg/kartograph/map.py", line 257, in _get_bounding_geometry
      charset=layer.options['charset']
    File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/kartograph.py-0.6.8-py2.7.egg/kartograph/layersource/shplayer.py", line 121, in get_features
      geom = shape2geometry(shp, ignore_holes=ignore_holes, min_area=min_area, bbox=bbox, proj=self.proj)
    File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/kartograph.py-0.6.8-py2.7.egg/kartograph/layersource/shplayer.py", line 153, in shape2geometry
      geom = shape2polygon(shp, ignore_holes=ignore_holes, min_area=min_area, proj=proj)
    File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/kartograph.py-0.6.8-py2.7.egg/kartograph/layersource/shplayer.py", line 170, in shape2polygon
      parts = shp.parts[:]
  AttributeError: _Shape instance has no attribute 'parts'

I have run the code using a different source file (U.S. county map) and the conversion worked, so I don't think there are any code or version problems.  I assume kartograph is looking for something in my source file that it can't find or doesn't exist.

Can anyone help me troubleshoot the AttributeError?

If it matters, I'm working on a Mac 10.9.5

UPDATE: following gene's guidance in the comments, I attempted to read my shapefile using pyshp (shapefile).  shapefile appears to read the source file just fine and appears to recognize the 'parts' attribute. Here is an example of my work:

>>> import shapefile
>>> sf = shapefile.Reader("TaxParcels.shp")
>>> shapes = sf.shapes()
>>> len(shapes)
70073
>>> for name in dir(shapes[45873]):
...     if not name.startswith('__'):
...             name
...
'bbox'
'parts'
'points'
'shapeType'
>>> shapes[876].shapeType
5
>>> shapes[876].parts
[0]
>>> fields = sf.fields
>>> len(fields)
206
>>> len(sf.records())
70073
>>>

I don't have enough experience yet to know what the next step is in troubleshooting.  I even tried stripping the "attributes" parameter from my script and I came up with the same error.  Could the problem be with how Kartograph is handling the shapefile?


texts: How to add one attribute column to another with a space between the two in QGIS

texts: I am wondering if it possible to add one attribute column to another with a space between the two in QGIS. I figured out you can add them together but I need a space between the texts


texts: Repeat polygon marker line by interval: neighboring polygons

texts:
As you can see in the screenshot below, I use a polygon marker line with interval. If there are 2 neighboring polygons, the markers will be doubled - of course.

But I need to suppress this effect. Any ideas?

EDIT: In my case, the polygons are in the SAME layer.


texts: Change coordinate system of a WMS service

texts:
I have a WMS service that is projected to UTM/SAD-69.

My workspace is in WGS-84 with no projection.

I want to change the coordinate system of my WMS service.

Is it possible in ArcMap 10.0?

The wms coordinate system is:

More information from the Layer Properties (using ArcGIS 10.2) that show it to have a Map Service Type of Cached.


texts:
I want to create a web application that takes in 5 images does some calculation and shows the result on a basemap. all these images would be on geotif format single band and with values on pixels, I want to show the output of all images on a basemap and download option  also.
Users will connect to my website and upload these files hit submit and would be redirected to next page where output would be displayed.

I have installed django on my server and developed front-end how to connect it to backend, for example a gdal operation in backend that will give output.


texts: How to remove the key control so users won't be able to use - and + to zoom?

texts:
I am using Openlayers on Drupal 7 and I'm trying to remove the key control so the person won't be able to use the key - and + to zoom.

for (map in Drupal.settings.openlayers.maps) {
    data = jQuery('#' + map).data('openlayers').openlayers;
    jQuery(data.controls).each(function(){
        data.removeControl((this));
    });
}

So with this code I delete all the controls, but in the end it's still working !

On the file openlayers_behavior_keyboarddefaults.js if I comment the line Drupal.openlayers.addControl(data.openlayers, 'KeyboardDefaults');

The keyboard zoom is not working !

In drupal there is an option witch can do that, but that not what I'm looking for.

If I don't find any solution I still can change the value in the DataBase but I would like not using this method, (I'd prefere to learn OL instead :))

Do you have any ideas ?

Thanks !


texts: Connecting to ArcGIS ImageServer in QGIS with custom renderingRule

texts:
How to connect ArcGIS ImageServer service in QGIS with custom 'renderingRule'?
Example service URL: https://ahn.arcgisonline.nl/arcgis/rest/services/Hoogtebestand/AHN3_i/ImageServer
renderingRule:

"AHN2 - Hillshade"


texts:
I made a QGIS plugin in QGIS 1.8 doing a lot of XSLT transformation and other XML handling with the library lxml. After upgrading to QGIS 2.0 this plugin does not work anymore:

"ImportError: No module named lxml.etree"

What libraries do one use to handle XML and XSLT in Python in QGIS 2.0?
Would it be a mess to try and get lxml working, in regard to distributing the plugin to other computers?


texts:
I am working with the v2.3.3 NuGet of GDAL by Tamas Szekeres for the .NET environment, and can successfully open the etopo1_ice_g_f4.flt raster to cut basic tiles.

The problem I am having is that even when I apply the -r bilinear option to my translation options, the generated image remains pixelated.

Specifically, I make the following calls in code to get the initial raster (simplifying the file names for this example):

GdalConfiguration.ConfigureGdal();

var etopoPath = @"F:\testing\etopo1_ice_g_f4.flt";
var colorPath = @"F:\testing\color.txt";
var pngPath = @"F:\testing\tile.png";
var tifPath = @"F:\testing\tile.tif";

var raster = Gdal.Open(etopoPath, Access.GA_ReadOnly);

I then set the appropriate projection extents, which will produce a 256x256 tile (in OpenMaps, it would be tile x = 306, tile y = 389, zoom = 10):

// shows us tile 306 389 z10
const double lonWest = -72.42;
const double lonEast = -72.02;
const double latNorth = 39.64;
const double latSouth = 39.24;

I set up the gdal_translate options to include bilinear interpolation:

var translateOptions = new[] {
    "-of", "GTiff",
    "-projwin", $"{lonWest}", $"{latNorth}", $"{lonEast}", $"{latSouth}",
    "-outsize", "256 256",
    "-r", "bilinear",
    "-q"
};

Gdal.wrapper_GDALTranslate(tifPath, raster, new GDALTranslateOptions(translateOptions), null, null);

And finally, process the GeoTIFF as a color-relief tile:

var demOptions = new[] {
    "-of", "PNG",
    "-q"
};

var tifRaster = Gdal.Open(tifPath, Access.GA_ReadOnly);
Gdal.wrapper_GDALDEMProcessing(pngPath, tifRaster, "color-relief", colorPath, new GDALDEMProcessingOptions(demOptions), null, null);

And it works, producing the following PNG image:

Now, the problem is that this does not seem to have been interpolated. The same image is produced even when -r bilinear is removed from the options.

Using the same elevation data (etopo1), and doing my own bilinear interpolation processing, I get this image:

I would like to not have to get the raw elevation data and do my own bilinear interpolation... I'd rather have the image properly interpolated coming out of gdal_translate.

So, am I missing something?


texts: Projecting latitude/longitude coordinates to Cartesian?

texts:
Let's say I have a square whose vertices are defined by coordinates in latitude/longitude format. I want to represent this square (and which ever shapes it contains, also defined as latitude/longitude) in a Cartesian graphic, such that the left-bottom vertex corresponds to point (0,0).

If this is possible, then how?


texts: Clipping a Geoserver raster layer with a polygon

texts:
I have an image mosaic raster layer with an extent that covers the contiguous United States. I also have a wfs layer with various polygons some covering subsets of the US and some outside of the US.

Essentially I'd like to have a request that returns a clipped version of the raster using one of these polygon.

As an attempt to do this I've installed the cross layer filtering extension

The following query either returns the entire raster or none of it based on if any part of the entire geometry of the raster intersects the polygon.

http://geoserverinstance/geoserver/wms?&layers=gdd:agdd&styles=gdd:agdd_web,&FORMAT=image/png&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&EXCEPTIONS=application/vnd.ogc.se_inimage&SRS=EPSG:4269&CQL_FILTER=INTERSECTS(the_geom, querySingle('gdd:redmaple', 'geom', 'gid=1'))&bbox=-125.020833333333,24.0625,-66.479166666662,49.937500000002&width=1600&height=800

In the query I just pass "the_geom" for the raster layer. Is there a way for it to do the intersection on a pixel by pixel basis masking out the pixels that don't intersect?

Is it not possible to do this clipping activity using this extension? I know there is also a wps extension that might do it, but I've had issues installing it. So if this won't work I can go back to trying to figure those issues out.


texts: Openlayers zoomToExtent pause until image is retrieved

texts:
Problem:

layer transition effect: resize
Do programmatic map.zoomToExtent() to somewhere OUTSIDE of current wms map bounds
OpenLayers setCenters to target bounds and zooms in.
While new wms image is being downloaded, OpenLayers displays 'enlarged' map
Once wms image is loaded, it's transitioned into map viewport

Is there a way to 'prevent' OpenLayers from set-centering at #3 until wms image is downloaded ?

What I want:

as above
as above
as above
OpenLayers doesn't setCenter nor zoom-in (or just doesn't update map)
if new WMS image is loaded, then update map.


texts:
I create an object derived from BaseCustomLayer which draws a polyline.
I call Map.addLayer(myBaseCustomLayer) and the function

public override void Draw(esriDrawPhase drawPhase, IDisplay Display, ITrackCancel trackCancel)

is hit.  BUT if I set DynamicMapEnabled = true the draw function is not hit !

This is driving me crazy.  Thanks for any hints.


texts:
I was looking on was NMEA sentences. From which I got a series of know how doubts.

I got to know about the terms like $GPGGA, GPHDT, GPGSV, GPGLL from GPS data. Also I could see many new terms like PASCD, PAGCD, PADCH.

What forms the source for these sentences, and what does these refer to ?
I couldn't get any information on these sentences, kindly help me out in this


texts: Add all the shapefile in a map from a folder

texts:
I'd like to know how can I add, in the map, all the shapefile which are in a specific folder, in one shot; using an algorithm of the Graphical Modeler of QGIS (if it exists). I noticed a similar algorithm, but it can only add a single shapefile.
I need this function using the graphical modeler and it should be a part of a complex model.


texts: How to render a raster multiple times across the date line?

texts:
Using EPSG:4326 and WGS84, I have a GeoTIFF with longitude extents (167, 208) - so it crosses the date line. I'd like MapServer to display the image regardless of whether the user's requesting bbox comes from the west or the east, however:

[...]bbox=150,0,250,90[...]

displays the image, while

[...]bbox=-210,0,-110,90[...]

displays a blank image.

Everything else unchanged, I'd like to see the same image returned in both cases. Some research brought up this 3+ year old mail from FrankW:

http://lists.osgeo.org/pipermail/mapserver-dev/2010-February/009713.html

Points #3 & 4 imply that adding PROCESSING "RESAMPLE=AVERAGE" should solve this problem, but I've added that to my mapfile with no luck. Has this changed in the past few years? Is there another way I can get the desired behavior? Is there a setting I could check for that could be preventing the image from mirroring every +/- 360 degrees longitude?


texts:
I have a large polygon (A) with a specific area.
I also have another polygon (B) with area less than polygon A.
I want to buffer the polygon B so that its area becomes equal to the area of A.
from shapely.geometry import Polygon
ulx,uly= -105.645292, 28.094511
lrx, lry = -101.985052, 24.885167
poly_A = Polygon([(ulx,uly),(ulx,lry),(lrx,lry),(lrx,uly)])
print (poly_A.area)

ulx,uly= -110.398408,  58.136267
lrx, lry = -108.689382,  57.137692
poly_B = Polygon([(ulx,uly),(ulx,lry),(lrx,lry),(lrx,uly)])
print (poly_B.area)

buffer_size = (poly_A.area - poly_B.area) / poly_B.length
print (buffer_size)

poly_C = poly_B.buffer(buffer_size)
print (poly_C.area)

After knowing buffer_size, I can get bufferred polygon C  with area equals to polygon A.
Note. All polygons dealt are rectangles.
Both polygons are in latlon coordinates


texts: How to merge feature class with the same name and save in other gdb?

texts:
I have tree feature class with the same name in different geodatabases for example:

one.gdb\fc1,fc2,fc3 | two.gdb\fc1,fc2,fc3  | tree.gdb\fc1,fc2,fc3

what I want to build a single merge in another geodatabase, for example, result.gdb\fc1_merge, fc2_merge. fc3_merge. ¿How could I build a script python to merge the fcs in one single fc in another geodatabase?

So far, this is my script:

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

arcpy.env.workspace = 'C:\FOLDER_GDB' #<--- All my gdb inside
inWorkspace = arcpy.env.workspace
OutPutWorkspaces = 'C:\FOLDER_RESULT\Result.gdb'
fc = 'fc1'
Matches = []
for root, dirs, files in os.walk(inWorkspace):
    for shps in files:
        if shps == 'fc1':
           match =(os.path.join(root, shps))
           matches.append (match)

arcpy.Merge_management(Matches,OutPutWorkspaces+fc+'_Merge')

print 'Done!!'

When I run the script I get de next error:

Traceback (most recent call last):
  File "<string>", line 254, in run_nodebug
  File "D:\SedesolVQB\Merge_prueba.py", line 18, in <module>
    arcpy.Merge_management(Matches,OutPutWorkspaces+fc+'_Merge.shp')
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\management.py", line 3816, in Merge
    raise e
arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000735: Input Datasets: Value is required
Failed to execute (Merge).

Does anybody can help me to fix my script


texts: Batch renaming of shapefiles using ArcPy and folder name where shapefile is stored?

texts:
I am very new to Python and wondered if anyone had any idea how to do this?

I have ~100+ shapefiles named 'land boundary' within the folder of their site name.

Is this possible using Python in ArcGIS?


texts: Calculating flood volume using LiDAR DEM and HEC-RAS floodplain polygon?

texts:
I’ve got a vexing GIS analysis problem which is
conceptually fairly straightforward, but in reality more challenging to figure out….

Basically, I need to calculate the VOLUME of floodwater within a
100 year floodplai, using ArcGIS 10.0 (all extensions available).

My input datasets are:

LiDAR DEM (1m pixel or 3m pixel) of bare ground (raster)
100 year floodplain (output of USACE HEC-RAS analysis),
Polygon Shapefile, but attribute table has no useful content.

So far I have overlaid the floodplain polygons over the LiDAR raster,
[using Extract by Mask tool from Spatial Analyst], so that I have the ground elevations are along the outside edge of the floodplain (think of the top of the bathtub).
Key concept is that the top of the floodplain gradually falls in elevation
as you move from the headwaters down through the watershed.

Then, (this seems to be the impossible part) I need to figure out
how I can take the elevations along the outside edge pixels of the
raster floodplain (output of previous Extract by Mask) and extend them to the
equivalent elevation on the other side of the floodplain, and
“raise” all the pixels in between the 2 outside edges of the floodplain
to the “top level” of this 100-year floodplain.

Once I have this “top elevation” of the floodplain water, I should
be able to use the CUT AND FILL tools in ArcGIS to simply subtract
the LiDAR ground elevations from the “top floodplain elevation”
and compute a volume between these two rasters.

I’ve been researching this for days and have yet to find a solution.

Perhaps I’m thinking of this completely wrong?

I would have expected to find others who have dealt with this problem,
but while there are similar problems, none appear to offer a solution.


texts: Creating drop down of derived fields from CSV using ArcPy?

texts:
I'm trying to create a drop down menu in a Python script in ArcMap where one input is a CSV and the other input is a drop down of the headers in that CSV. Is this possible at all?

In the image above, "Candidate List" is the supplied CSV and "Candidate Field" is the field header that the user selects for the column that the data is located in. Would I have to put something in the script its self or in the update parameters?

with open(Candidate List, "rb") as f:
    reader = csv.reader()
    i = reader.next()

Wherein I want I to be the fields for that drop down.


texts: Using different criterion (CO2 emissions, fuel consumption) on network analyst

texts:
i just started to use the network analysis tool from QGIS.

The problem is that the only criterion available to make the "shortest" path are time and length.

I have a layer with information like CO2 emissions, fuel consumption etc, how do I route the path that has less fuel consumption, less CO2 emissions?

I just want to add as a criterion the CO2 emissions and fuel consumption fields if possible.


texts:
I'm trying to learn how to use ps.map in GRASS. I thought that running it in interactive mode from the command line would be a good place to start, but I'm getting the following error:

> ps.map                                                                                ~
Traceback (most recent call last):
  File "/Applications/GRASS-6.4.app/Contents/MacOS/etc/wxpython/gui_modules/menuform.py", line 59, in <module>
    import globalvar
  File "/Applications/GRASS-6.4.app/Contents/MacOS/etc/wxpython/gui_modules/globalvar.py", line 76, in <module>
    import wx
  File "/Applications/GRASS-6.4.app/Contents/MacOS/etc/python/wx/__init__.py", line 45, in <module>
    from wx._core import *
  File "/Applications/GRASS-6.4.app/Contents/MacOS/etc/python/wx/_core.py", line 4, in <module>
    import _core_
ImportError: dlopen(/Applications/GRASS-6.4.app/Contents/MacOS/etc/python/wx/_core_.so, 2): no suitable image found.  Did find:
    /Applications/GRASS-6.4.app/Contents/MacOS/etc/python/wx/_core_.so: mach-o, but wrong architecture

UPDATE: In response to Markus's comment that terminal messages on GRASS startup might help to diagnose the problem, the following two chunks provides all terminal output from GRASS startup.

Chunk 1: from starting /Applications/Grass-6.4.app to GUI Startup window.

Last login: Sat Mar 23 10:08:40 on ttys002
gregory@Gregorys-MacBook> '/Applications/GRASS-6.4.app/Contents/MacOS/grass.sh'; exit
Rebuilding Addon HTML manual pages index...
Rebuilding Addon menu...
dyld: DYLD_ environment variables being ignored because main executable (/usr/bin/osascript) is code signed with entitlements
Python 2.6.7 found.
Cleaning up temporary files ...
Starting GRASS ...

Chunk 2: from selecting a location/mapset in the GUI startup window and clicking "Start GRASS".

          __________  ___   __________    _______________
         / ____/ __ \/   | / ___/ ___/   / ____/  _/ ___/
        / / __/ /_/ / /| | \__ \\_  \   / / __ / / \__ \
       / /_/ / _, _/ ___ |___/ /__/ /  / /_/ // / ___/ /
       \____/_/ |_/_/  |_/____/____/   \____/___//____/

Welcome to GRASS 6.4.2 (2012)
GRASS homepage:                          http://grass.osgeo.org/
This version running thru:               shell (/bin/zsh)
Help is available with the command:      g.manual -i
See the licence terms with:              g.version -c
If required, restart the GUI with:       g.gui wxpython
When ready to quit enter:                exit

gregory@Gregorys-MacBook> Mar 23 11:31:05 Gregorys-MacBook.local pythonw2.6[2520] <Error>: CGContextRestoreGState: invalid context 0x0
Mar 23 11:31:05 Gregorys-MacBook.local pythonw2.6[2520] <Error>: CGContextRestoreGState: invalid context 0x0
Mar 23 11:31:05 Gregorys-MacBook.local pythonw2.6[2520] <Error>: CGContextRestoreGState: invalid context 0x0
Mar 23 11:31:05 Gregorys-MacBook.local pythonw2.6[2520] <Error>: CGContextRestoreGState: invalid context 0x0
Mar 23 11:31:05 Gregorys-MacBook.local pythonw2.6[2520] <Error>: CGContextRestoreGState: invalid context 0x0


texts: Using ArcGIS JavaScript API Directions with embedded coordinates?

texts:
I want to use ArcGis Javascript API for directions with the users current loaction as the start location, and a predetermined location as the stop destination. All the examples I can find use a widget for start and stop locations.

How can I set a predetermined location for the destination, and a current location for the start rather than using the widgets?


texts: Removing identical items from grouped layers in legend of ArcMap?

texts:
I am using ArcMap 10.5 and in the process of creating a map template for people to use in my team. Part of the template is a map layout with a legend. The problem I am having is that I have some layers which are grouped together (Population points, as seen below) and all have the same symbology. In each map created by users of the template, they may have none, one or more than one of the features in this group activated.

The problem I am having is that if all of them are switched on, they all show up in the legend. I have tried adding only one of the layers to the Legend, but this means that if the end user doesn't add that specific layer to their layout then the legend will not display the key for these groups of items.

I am trying to find a way to label the legend so that if one or more of these grouped layers is switched on, then only one label will appear in the legend (i.e. no duplicates).

Essentially this is what I would like the end product to look like:


texts:
I need to create an Image with a Geometry Feature Property which is a String so that I can map an ImageCollection with that property.
I tried reduceToVectors but the output vector doesn't have the geometry property I need....


texts:
I have a couple of QGIS applications (one 3.2.2-Bonn and one 2.18.9 Las Palmas) running on Ubuntu boxes (both 18.04.1 LTS).

Files are loaded from a SAMBA file share and load fine. But when I try to save the projects, QGIS gives an error:

  couldn't save to the location or there was a zip error

If I try a save-as it says the same thing and creates a file of zero bytes.

Is this a QGIS bug or am I doing wrong when I mount the file share?

I mount them through the 'other locations' tab in the Nautilus.

I should have noted that the SAMBA is a windows server and the files were all created in the windows version of QGIS - not sure if that makes a difference.

The terminal command produces this:

And the error has changed but now looks like this:


texts: What is the function of the ESRI.ArcGIS.AddIns.targets file?

texts:
I am currently upgrading from 10.2 to 10.3 and as expected there are always some issues with recompiling addins and upgrading the the visual studio version.

There is now a known issue for which there is a workaround described here.

My question is what is the function of the ESRI.ArcGIS.AddIns.targets file?  What does it do and why does renaming the new ESRI.ArcGIS.AddIns.11.targets to the former ESRI.ArcGIS.AddIns.targets work?  Is it what packages the actual dll for use with ArcGIS applications?


texts: Iteratively export raster layers along with basemap in ArcMap

texts:
My MXD has ~50 raster layers and basemap layer (from Esri basemap). The objective is to export each raster layer + basemap layer into PNG files one by one.

In Export to PNG documentation I am not seeing any parameter to select layers I want to export so a simple loop isn't a solution.

Thus I am looking for any guidance on iteratively exporting layers along with basemap.


texts: Automate exporting one image by year for several years GEE

texts:
I am using the code below to export one image with 12 bands for year the year 2015, if want to generate the images for each year from 2000 on, is there any way to automate this to generate one image per year, without having to change the dates of the code each time?

// Export Burned Area Australia
//

var col =  ee.ImageCollection("MODIS/006/MCD64A1")
             .filterDate('2015-01-01', '2016-01-01')
             .select('BurnDate');

// Define the regional bounds
var region = ee.Geometry.Polygon(
        [[[139.81147500927887, -22.584034738321385],
          [139.81147500927887, -39.52984428265981],
          [155.19233438427887, -39.52984428265981],
          [155.19233438427887, -22.584034738321385]]], null, false);

Map.addLayer(region)

// Clip and add a date band
var clipToRegion = function(img) {
  var dateString = ee.Date(img.get('system:time_start')).format('yyyy-MM-dd');
  var clipped = img.clip(region)
  return clipped.rename(dateString)
};

var burned_area = col.map(clipToRegion);
var scal = burned_area.first().projection().nominalScale();
var check = ee.Image(burned_area.first());
Map.addLayer(check, {palette: ['000000', '00FFFF'], max: 366}, 'check');
print('Scale', scal )

// Stack One layer by year
//
var stackCollection = function(collection) {
  // Create an initial image.
  var first = ee.Image(collection.first()).select([]);

  // Write a function that appends a band to an image.
  var appendBands = function(image, previous) {
      var dateString = ee.Date(image.get('system:time_start')).format('yyyy-MM-dd');
      return ee.Image(previous).addBands(image);
  };
  return ee.Image(collection.iterate(appendBands, first));
};

var evi_img = stackCollection(burned_area);
print("EVI image stack",evi_img);

Export.image.toDrive({
  image: evi_img,
  description: 'BurnedAreaAustralia',
  scale: 463.3127165275,
  region: region,
  fileFormat: 'GeoTIFF',
  formatOptions: {
    cloudOptimized: true
  }
});


texts: Raster turns into 4 separate rasters after clipping

texts:
I have a raster with a .sid extension.  After I clip my raster, it produces 4 separate new rasters (3 of which are black and white).  When I delete any of the Black and White rasters, the clipped raster that's in color no longer works.  Is there a way to get around this, or at least combine them all together into a single raster?

I have tried the Composite Bands tool, and that comes up with an error.


texts:
I run my exe on cmd by passing argument and get this error:

Cannot mix incompatible Qt library (version 0x50902) with this library (version 0x50901)

How can I fix the problem?

It is an standalone application and I give the path on cmd window at the beginning like that (my exe file name is "path2Shape"):

path2Shape.py C:\Users\Mustafa Uçar\Desktop\Tutorial\Mississippi\mississippi.shp

Here is my code of the exe file:

import sys
from qgis.core import *
from qgis.gui import *
from PyQt5.QtCore import *
from PyQt5 import QtGui

app = QgsApplication([], True)
path = "C:/OSGeo4W64/apps/qgis-dev"
app.setPrefixPath(path, True)
app.initQgis()

canvas = QgsMapCanvas()
title = "PyQGIS Standalone Application Example"
canvas.setWindowTitle(title)
canvas.setCanvasColor(Qt.white)

filepath = sys.argv[0]

layer = QgsVectorLayer(filepath, "testlayer", "ogr")
QgsProject.instance().addMapLayer(layer)
canvas.setExtent(layer.extent())
canvas.setLayers([layer])
canvas.zoomToFullExtent()

canvas.show()
sys.exit(app.exec_())


texts: Centralizing shapefile data into database

texts:
I've got hundreds of shapefiles from various different GIS projects that I want to start consolidating into a single database platform, currently attempting this with Postgres  / PostGIS.

Hardly any of the data is standardized - meaning that it's a lot of the same data types, but the particular attribute names / types don't match.

Where should I start tackling this?  Should I develop a standard model to migrate each shapefile into first (e.g. Hydro_line, transport_line, Hydro_poly standards, etc)?

An alternative is to just import each shapefile into Postgres individually, so each shp becomes a table in the database, but I'm not sure about this in terms of performance and organization.  Feels sort of like delaying the inevitable...

Any advice on dealing with this daunting task?


texts: GDAL error - no create method for format - QGIS

texts:
I am trying to clip a raster layer. I am getting the error..

ERROR 6: GDALDriver::Create() ... no create method implemented for this format.

The documentation (link below) says to use GeoTIFF (gtiff). The format GeoTIFF exists, but the extensions are .tif .tiff .TIF & .TIFF.

What am I missing?

gdal.org/gdal_translate

Edit:
Raster > Extraction > Clipper... Select a region. Give it a name. Hit OK. Note the tiff file types available.


texts: Google Earth Engine: calculate number of weeks where SST over a threshold in past 12 weeks

texts:
I've been using the MODIS-Aqua/L3SMI dataset to investigate sea surface temperatures (SST) for a specific region.

What I'd like to do is for a specific time epriod (could be day, week or month), find out how many of the previous weeks have had a SST anomaly of > 1 degree over the average.

To start with I've calculated the mean SST for the region over the previous ten years.

var collection = ee.ImageCollection('NASA/OCEANDATA/MODIS-Aqua/L3SMI')
  .select('sst')
  .map(function(image){return image.clip(EEZ)});

// Define reference conditions from the first 10 years of data.
 var reference = collection.filterDate('2002-01-01', '2012-12-31')

// Sort chronologically in descending order.
  .sort('system:time_start', false);

print('reference', reference)

// Compute the mean of the first 10 years.
 var mean = reference.mean();

 print('ref mean', mean)

Then subtracted than mean from the mean SST from another series of images from the region to get the SST anomaly.

var series = collection.filterDate('2016-01-01', '2018-12 31').map(function(image) {
    return image.subtract(mean).set('system:time_start',
 image.get('system:time_start'));
 });

 print('SSTanomaly', series)

From here however, I'm not sure how to calculate how many consecutive weeks this anomaly is 1 or more before a specific date (could be a day, week or month).


texts: Adding nodes automatically without clicking the mouse

texts:
Is there a way in QGIS (or maybe in another GIS (ArcMap, for example))) to automatically create nodes at a given interval when digitizing a raster?
For example, I decided to create a Polyline object. I created a starting point and then just move the mouse in the desired direction. At this time QGIS is listening for the line length. If the length of the drawn line reaches 10m, a node is created automatically in the place of the mouse pointer. At the end of procedure, as usual, right click to complete the polyline.

I attach this pseudocode:
my_length = 10
start_point = create_point_on_left_click_event(get_coords_from_mouse_pointer())
while !right_click_event:
    #line_length_listetning...
    on_event(line_length == my_length):
        node = create_node(get_coords_from_mouse_pointer())
        new_start_point = node

end_point = create_point_on_right_click_event(get_coords_from_mouse_pointer())
#completion Polyline


texts: ArcMap remote desktop connection error with Windows10

texts:
Issue with ArcMap 10.5 (and above) with Windows10. For some reason I have been getting the following error when trying to start an instance during a remote desktop session (most probably related to the graphics card).

  LoadLibrary failed with error 87:The parameter is incorrect.


texts: How to force QueryTask to always POST requests?

texts:
On ArcGIS API for JavaScript, how can we force QueryTask to always perform POST requests?

I know that I can force POST for esriRequest via the usePost options parameter. Perhaps a similar, undocumented, method exists for the rest of the API?

My application specifications require all requests to be performed via POST. So, a solution applicable to all ArcGIS API for JavaScript requests would be even better.


texts: Find nearest polygon (from GeoSeries) to point (from GeoSeries)

texts:
I have two GeoSeries, consisting of points and polygons.  I want to find the polygon in dataframe B that is closest to each point in dataframe A.  The polygons are rooftops from https://github.com/Microsoft/USBuildingFootprints, which I already geocoded using https://github.com/Bonsanto/polygon-geohasher.

I'm currently computing the 7 digit geohash of each point, and merging on buildings in neighboring 7 digit geohashes using geotools.expand.  This is better than doing a full outer merge, but relies on explode.  My general approach was to minimize calls of distance, since computing the distance from a point to a polygon is expensive.

The code is a bit slow (~20 minutes to match 100k rows), and so I'm trying to make it faster.  My searching points to r-trees, but the sklearn implementation seems to be geared towards identifying the nearest point, rather than the nearest polygon. I'm interested in the left join rather than the right join.

Code below:

import pandas as pd
import numpy as np
import geopandas
import geohash
from shapely.geometry import Point

def match_func(df):
    point = Point(df.iloc[0,:][['lat', 'long']])
    df.loc[:, 'dist'] = geopandas.GeoSeries(df.geometry).distance(point)
    df = df.sort_values('dist')
    return(df.head(1))

def main(file):

    x           = import_points()
    rooftop_df  = import_rooftops()

    x['id'] = range(1, len(x) + 1)

    def neighbor_fun(lat,long):
         return(geohash.encode(lat,long,precision=7))

    func1 = np.vectorize(neighbor_fun)

    x['g7_neighbor'] = func1(x['lat'], x['long'])
    x = x.explode('g7_neighbor')
    x = x.merge(rooftop_df, left_on='g7_neighbor', right_on='geo7')

    xg = x.groupby('id')
    xout = pd.concat([match_fun2(df_group) for group_name, df_group in xg])
    return(xout)


texts: Adding buffer to tracklog in real-time in ArcPad?

texts:
Hardware: Trimble Juno T41/5
Software: ArcPad 10.2

I am looking for a method or plugin for ArcPad that will help with spraying for noxious weed control, specifically a way to add a "boom width" buffer to a track log and display it in real time. This plugin should allow the user to identify gaps and overlap in chemical application. I know there are several stand alone systems that do similar functions but I already own the Trimble handheld and don't really want to spend more for a very basic function.


texts: MultiLineString to separate individual lines using Python with GDAL/OGR, Fiona, Shapely

texts:
There are lots of MultiLineString features in a polyline shapefile. The problem is to split these lines into separate individual lines. There is a similar question with PostGIS, but I need to code in Python.

I have been planning to extract the geometry from each MultiLineString feature and create separate features and append a shpfile with those new features. Later, in the second pass, I would delete the MultiLineString object which is no more required as it is now split into its components. It will require reading data from original shpfile and appending the newly created features to a copy of the shpfile. Then the appended (copy) shpfile needs to be rid of the MultiLineString features. It seems cumbersome as I find no way to read and write the shapefile simultaneously.

There should be better ideas to do that. How to split these lines into separate individual lines using Python with OGR, Shapely, Fiona or Pyshp? A snippet of code to give an idea is below.

import fiona
from shapely.geometry import LineString

original_shpfile = fiona.open(path_to_road_network_shapefile)
with fiona.open(path_to_road_network_shapefile, 'a') as copy_shpfile:
    rec = original_shpfile.next()
    while rec:
        try:
            vertices_on_link = rec['geometry']['coordinates']
            LineString(vertices_on_link) #reconstruct the line from vertices
            #if multipart,
            #has multiple lists within the list,
            #one for each of the Linestring component. also throws AssertionError

        except AssertionError:
            number_of_parts = len(vertices_on_link)
            for x in range(0,number_of_parts):
                part = rec
                part['geometry']['coordinates'] = vertices_on_link[x]
                copy_shpfile.write(part)
        rec = original_shpfile.next()

#in the next phase, delete each of the the MultiLineString features from the copy of the shapefile


texts: QGIS 3.16+ transparency layer export to pdf not working

texts:
I've got problem with transparency when exporting my layout to PDF file. I have layer with transparency defined in Symbology -> Layer Rendering -> Opacity.

In map view and in layout i can see transparency in my layout:

However, adter exporting to PDF transparency is gone and edge of polygon looks weird:

My PDF export settings are:

What solve this issue, but I can't use it because of other reasons:

Checking "Print as raster" in Layout Settings

Define transparency right in the color of symbol, not in "Opacity" of layer.
Exporting as image
using older varsion of QGIS (works in 3.2)

As I've mentioned, i can't use any of this solution.
What I've tried:

Enabling/Disabling "simplify geometry"
Enabling/Disabling "Always export as vector"
Enabling/Disabling "Export RDF metadata"
messing up with other layer-, layout- and export settings. All things that works I've mentioned above.

Someone has this exact problem earlier, but no satisfying answer was provided
Do you have any ideas how to fix that issue? Am I doing something wrong?